From a475fa6faff47c2342548967c124daca8605bbdf Mon Sep 17 00:00:00 2001 From: Stergios Terpinas Date: Thu, 15 Jun 2023 13:02:51 +0300 Subject: [PATCH] Restructure gammatone implementation using Eigen. --- src/gammatone_filterbank.cc | 144 +++++++++++------- src/gammatone_spectrogram_builder.cc | 43 +++--- .../equivalent_rectangular_bandwidth.h | 17 +++ src/include/gammatone_filterbank.h | 33 ++-- tests/gammatone_filterbank_test.cc | 14 +- 5 files changed, 154 insertions(+), 97 deletions(-) diff --git a/src/gammatone_filterbank.cc b/src/gammatone_filterbank.cc index 86832266..44fff6c1 100644 --- a/src/gammatone_filterbank.cc +++ b/src/gammatone_filterbank.cc @@ -16,81 +16,113 @@ #include #include +#include #include "amatrix.h" +#include "equivalent_rectangular_bandwidth.h" #include "signal_filter.h" +using namespace Eigen; + namespace Visqol { GammatoneFilterBank::GammatoneFilterBank(const size_t num_bands, const double min_freq) - : num_bands_(num_bands), - min_freq_(min_freq), - fltr_cond_1_({0.0, 0.0}, num_bands), - fltr_cond_2_({0.0, 0.0}, num_bands), - fltr_cond_3_({0.0, 0.0}, num_bands), - fltr_cond_4_({0.0, 0.0}, num_bands) {} + : num_bands_(num_bands), min_freq_(min_freq) + { + fltr_cond_1_ = MatrixXd::Zero(num_bands, kFilterLength); + fltr_cond_2_ = MatrixXd::Zero(num_bands, kFilterLength); + fltr_cond_3_ = MatrixXd::Zero(num_bands, kFilterLength); + fltr_cond_4_ = MatrixXd::Zero(num_bands, kFilterLength); + a1 = MatrixXd::Zero(num_bands, kFilterLength); + a2 = MatrixXd::Zero(num_bands, kFilterLength); + a3 = MatrixXd::Zero(num_bands, kFilterLength); + a4 = MatrixXd::Zero(num_bands, kFilterLength); + b = MatrixXd::Zero(num_bands, kFilterLength); + } size_t GammatoneFilterBank::GetNumBands() const { return num_bands_; } double GammatoneFilterBank::GetMinFreq() const { return min_freq_; } void GammatoneFilterBank::ResetFilterConditions() { - fltr_cond_1_ = {{0.0, 0.0}, num_bands_}; - fltr_cond_2_ = {{0.0, 0.0}, num_bands_}; - fltr_cond_3_ = {{0.0, 0.0}, num_bands_}; - fltr_cond_4_ = {{0.0, 0.0}, num_bands_}; + fltr_cond_1_.fill(0.0); + fltr_cond_2_.fill(0.0); + fltr_cond_3_.fill(0.0); + fltr_cond_4_.fill(0.0); } void GammatoneFilterBank::SetFilterCoefficients( - const AMatrix& filter_coeffs) { - fltr_coeff_A0_ = filter_coeffs.GetColumn(0).ToValArray(); - fltr_coeff_A11_ = filter_coeffs.GetColumn(1).ToValArray(); - fltr_coeff_A12_ = filter_coeffs.GetColumn(2).ToValArray(); - fltr_coeff_A13_ = filter_coeffs.GetColumn(3).ToValArray(); - fltr_coeff_A14_ = filter_coeffs.GetColumn(4).ToValArray(); - fltr_coeff_A2_ = filter_coeffs.GetColumn(5).ToValArray(); - fltr_coeff_B0_ = filter_coeffs.GetColumn(6).ToValArray(); - fltr_coeff_B1_ = filter_coeffs.GetColumn(7).ToValArray(); - fltr_coeff_B2_ = filter_coeffs.GetColumn(8).ToValArray(); - fltr_coeff_gain_ = filter_coeffs.GetColumn(9).ToValArray(); + const AMatrix &filter_coeffs) { + for (size_t i = 0; i < num_bands_; ++i) { + Eigen::RowVector3d a1_i = {filter_coeffs(i, ErbFiltersResult::A0) / filter_coeffs(i, ErbFiltersResult::Gain), + filter_coeffs(i, ErbFiltersResult::A11) / filter_coeffs(i, ErbFiltersResult::Gain), + filter_coeffs(i, ErbFiltersResult::A2) / filter_coeffs(i, ErbFiltersResult::Gain)}; + a1.row(i) = a1_i; + + Eigen::RowVector3d a2_i = {filter_coeffs(i, ErbFiltersResult::A0), filter_coeffs(i, ErbFiltersResult::A12), filter_coeffs(i, ErbFiltersResult::A2)}; + a2.row(i) = a2_i; + + Eigen::RowVector3d a3_i = {filter_coeffs(i, ErbFiltersResult::A0), filter_coeffs(i, ErbFiltersResult::A13), filter_coeffs(i, ErbFiltersResult::A2)}; + a3.row(i) = a3_i; + + Eigen::RowVector3d a4_i = {filter_coeffs(i, ErbFiltersResult::A0), filter_coeffs(i, ErbFiltersResult::A14), filter_coeffs(i, ErbFiltersResult::A2)}; + a4.row(i) = a4_i; + + Eigen::RowVector3d b_i = {filter_coeffs(i, ErbFiltersResult::B0), filter_coeffs(i, ErbFiltersResult::B1), filter_coeffs(i, ErbFiltersResult::B2)}; + b.row(i) = b_i; + } } -AMatrix GammatoneFilterBank::ApplyFilter( - const std::valarray& signal) { - AMatrix output(num_bands_, signal.size()); - std::valarray a1, a2, a3, a4, b; +void BatchFilterSignal(const MatrixXd &numer_coeffs, + const MatrixXd &denom_coeffs, + const MatrixXd &input_signal, + MatrixXd &output_signal, + MatrixXd &init_conditions) { + size_t i, n = denom_coeffs.cols(); + for (size_t m = 0; m < input_signal.rows(); m++) { + output_signal.col(m) = numer_coeffs.col(0) * input_signal(m, 0) + init_conditions.col(0); + for (i = 1; i < n; i++) { + init_conditions.col(i - 1) = numer_coeffs.col(i) * input_signal(m, 0) + + init_conditions.col(i) - denom_coeffs.col(i).cwiseProduct(output_signal.col(m)); + } + } +} - // Loop over each filter coefficient now to produce a filtered column. - for (size_t chan = 0; chan < num_bands_; chan++) { - a1 = {fltr_coeff_A0_[chan] / fltr_coeff_gain_[chan], - fltr_coeff_A11_[chan] / fltr_coeff_gain_[chan], - fltr_coeff_A2_[chan] / fltr_coeff_gain_[chan]}; - a2 = {fltr_coeff_A0_[chan], fltr_coeff_A12_[chan], fltr_coeff_A2_[chan]}; - a3 = {fltr_coeff_A0_[chan], fltr_coeff_A13_[chan], fltr_coeff_A2_[chan]}; - a4 = {fltr_coeff_A0_[chan], fltr_coeff_A14_[chan], fltr_coeff_A2_[chan]}; - b = {fltr_coeff_B0_[chan], fltr_coeff_B1_[chan], fltr_coeff_B2_[chan]}; - - // First filter - auto fltr_rslt = SignalFilter::Filter(a1, b, signal, fltr_cond_1_[chan]); - fltr_cond_1_[chan] = fltr_rslt.finalConditions; - - // Second filter - fltr_rslt = SignalFilter::Filter(a2, b, std::move(fltr_rslt.filteredSignal), - fltr_cond_2_[chan]); - fltr_cond_2_[chan] = fltr_rslt.finalConditions; - - // Third filter - fltr_rslt = SignalFilter::Filter(a3, b, std::move(fltr_rslt.filteredSignal), - fltr_cond_3_[chan]); - fltr_cond_3_[chan] = fltr_rslt.finalConditions; - - // Fourth filter - fltr_rslt = SignalFilter::Filter(a4, b, std::move(fltr_rslt.filteredSignal), - fltr_cond_4_[chan]); - fltr_cond_4_[chan] = fltr_rslt.finalConditions; - - output.SetRow(chan, std::move(fltr_rslt.filteredSignal)); +void BatchFilterBands(const MatrixXd &numer_coeffs, + const MatrixXd &denom_coeffs, + const MatrixXd &input_signal, + MatrixXd &output_signal, + MatrixXd &init_conditions) { + size_t i, n = denom_coeffs.cols(); + for (size_t m = 0; m < input_signal.cols(); m++) { + output_signal.col(m) = numer_coeffs.col(0).cwiseProduct(input_signal.col(m)) + init_conditions.col(0); + for (i = 1; i < n; i++) { + init_conditions.col(i - 1) = numer_coeffs.col(i).cwiseProduct(input_signal.col(m)) + + init_conditions.col(i) - denom_coeffs.col(i).cwiseProduct(output_signal.col(m)); + } } - return output; +} + +MatrixXd GammatoneFilterBank::ApplyFilter(const MatrixXd &signal) { + + MatrixXd helper1 = Eigen::MatrixXd::Zero(num_bands_, signal.rows()); + MatrixXd helper2 = Eigen::MatrixXd::Zero(num_bands_, signal.rows()); + + // Loop over each filter coefficient now to produce a filtered column. + // First filter + Visqol::BatchFilterSignal(a1, b, signal, helper1, fltr_cond_1_); + + // Second filter + Visqol::BatchFilterBands(a2, b, helper1, helper2, fltr_cond_2_); + + // Third filter + helper1.fill(0.0); + Visqol::BatchFilterBands(a3, b, helper2, helper1, fltr_cond_3_); + + // Fourth filter + helper2.fill(0.0); + Visqol::BatchFilterBands(a4, b, helper1, helper2, fltr_cond_4_); + + return helper2; } } // namespace Visqol diff --git a/src/gammatone_spectrogram_builder.cc b/src/gammatone_spectrogram_builder.cc index 6314f8eb..c539c81a 100644 --- a/src/gammatone_spectrogram_builder.cc +++ b/src/gammatone_spectrogram_builder.cc @@ -27,6 +27,10 @@ #include "signal_filter.h" #include "spectrogram.h" +#include "Eigen/Dense" + +using namespace Eigen; + namespace Visqol { const double GammatoneSpectrogramBuilder::kSpeechModeMaxFreq = 8000.0; @@ -63,32 +67,33 @@ absl::StatusOr GammatoneSpectrogramBuilder::Build( " required minimum).")); } size_t num_cols = 1 + floor((sig.NumRows() - window.size) / hop_size); - AMatrix out_matrix(filter_bank_.GetNumBands(), num_cols); - auto sig_val_arr = sig.GetColumn(0).ToValArray(); - for (size_t i = 0; i < out_matrix.NumCols(); i++) { - const size_t start_col = i * hop_size; - // Select the next frame from the input signal to filter. - const std::slice_array frame = - sig_val_arr[std::slice(start_col, window.size, 1)]; + MatrixXd out_eigen (filter_bank_.GetNumBands(), num_cols); + MatrixXd frame_eigen (window.size, 1); + MatrixXd hann_window (window.size, 1); + for (int i = 0; i < window.size; ++i) { + hann_window(i, 0) = 0.5 - (0.5 * cos(2.0 * M_PI * i / (window.size - 1))); + } + + // run the windowing + for (size_t i = 0; i < out_eigen.cols(); i++) { + const size_t start_row = i * hop_size; + // select the next frame from the input signal to filter. + frame_eigen = sig.GetBackendMat().col(0).middleRows(start_row, window.size); // Apply a Hann window to reduce artifacts. - const std::valarray windowed_frame = window.ApplyHannWindow(frame); + frame_eigen.array() *= hann_window.array(); - // Apply the filter. + // apply the filter filter_bank_.ResetFilterConditions(); - auto filtered_signal = filter_bank_.ApplyFilter(windowed_frame); - // Calculate the mean of each row. - std::transform(filtered_signal.begin(), filtered_signal.end(), - filtered_signal.begin(), - [](decltype(*filtered_signal.begin())& d) { return d * d; }); - AMatrix row_means = filtered_signal.Mean(kDimension::ROW); - std::transform(row_means.begin(), row_means.end(), row_means.begin(), - [](decltype(*row_means.begin())& d) { return sqrt(d); }); - // Set this filtered frame as a column in the spectrogram. - out_matrix.SetColumn(i, std::move(row_means)); + MatrixXd filtered_signal = filter_bank_.ApplyFilter(frame_eigen); + + // calculate the mean of each row + out_eigen.col(i) = filtered_signal.array().square().rowwise().mean().sqrt(); } + AMatrix out_matrix(out_eigen); + // Order the center freq bands from lowest to highest. std::vector ordered_cfb; ordered_cfb.reserve(erb_rslt.centerFreqs.size()); diff --git a/src/include/equivalent_rectangular_bandwidth.h b/src/include/equivalent_rectangular_bandwidth.h index fcf01c8b..41c3a051 100644 --- a/src/include/equivalent_rectangular_bandwidth.h +++ b/src/include/equivalent_rectangular_bandwidth.h @@ -27,6 +27,23 @@ namespace Visqol { * creation. */ struct ErbFiltersResult { + /** + * An enumeration that assigns the ERB filters coefficient names to the corresponding + * column indices of the ErbFiltersResult::filterCoeffs matrix. + */ + enum CoeffsMap { + A0 = 0, + A11 = 1, + A12 = 2, + A13 = 3, + A14 = 4, + A2 = 5, + B0 = 6, + B1 = 7, + B2 = 8, + Gain = 9 + }; + /** * The filter coefficients for the ERB filter. */ diff --git a/src/include/gammatone_filterbank.h b/src/include/gammatone_filterbank.h index 4102bff9..4ac75d1d 100644 --- a/src/include/gammatone_filterbank.h +++ b/src/include/gammatone_filterbank.h @@ -22,6 +22,8 @@ #include "amatrix.h" +#include "Eigen/Dense" + namespace Visqol { /** @@ -66,7 +68,7 @@ class GammatoneFilterBank { * * @return The filtered output. */ - AMatrix ApplyFilter(const std::valarray& signal); + Eigen::MatrixXd ApplyFilter(const Eigen::MatrixXd &signal); /** * Set the equivalent rectangular bandwidth filter coefficients that are to @@ -95,20 +97,21 @@ class GammatoneFilterBank { */ double min_freq_; - std::valarray> fltr_cond_1_; - std::valarray> fltr_cond_2_; - std::valarray> fltr_cond_3_; - std::valarray> fltr_cond_4_; - std::valarray fltr_coeff_A0_; - std::valarray fltr_coeff_A11_; - std::valarray fltr_coeff_A12_; - std::valarray fltr_coeff_A13_; - std::valarray fltr_coeff_A14_; - std::valarray fltr_coeff_A2_; - std::valarray fltr_coeff_B0_; - std::valarray fltr_coeff_B1_; - std::valarray fltr_coeff_B2_; - std::valarray fltr_coeff_gain_; + Eigen::MatrixXd fltr_cond_1_; + Eigen::MatrixXd fltr_cond_2_; + Eigen::MatrixXd fltr_cond_3_; + Eigen::MatrixXd fltr_cond_4_; + Eigen::MatrixXd a1; + Eigen::MatrixXd a2; + Eigen::MatrixXd a3; + Eigen::MatrixXd a4; + Eigen::MatrixXd b; + + /** + * This is dependent to the filters order. For a 2nd order filter + * we need 3 coeffs for numerator and denominator. + */ + const int kFilterLength = 3; }; } // namespace Visqol diff --git a/tests/gammatone_filterbank_test.cc b/tests/gammatone_filterbank_test.cc index 1ab8ae0e..0940f4b2 100644 --- a/tests/gammatone_filterbank_test.cc +++ b/tests/gammatone_filterbank_test.cc @@ -25,13 +25,13 @@ const size_t kSampleRate = 48000; const size_t kNumBands = 32; const double kMinFreq = 50; -// The contents of this input signal are random figures. -const AMatrix k10Samples{ - std::valarray{0.2, 0.4, 0.6, 0.8, 0.9, 0.1, 0.3, 0.5, 0.7, 0.9}}; - // Ensure that the output matrix from the signal filter has the correct // dimensions. TEST(ApplyFilterTest, basic_positive_flow) { + // The contents of this input signal are random figures. + Eigen::MatrixXd k10Samples (10, 1); + k10Samples << 0.2, 0.4, 0.6, 0.8, 0.9, 0.1, 0.3, 0.5, 0.7, 0.9; + auto filter_bank = GammatoneFilterBank{kNumBands, kMinFreq}; auto erb = EquivalentRectangularBandwidth::MakeFilters( kSampleRate, kNumBands, kMinFreq, kSampleRate / 2); @@ -39,10 +39,10 @@ TEST(ApplyFilterTest, basic_positive_flow) { filter_coeffs = filter_coeffs.FlipUpDown(); filter_bank.SetFilterCoefficients(filter_coeffs); auto filtered_signal = - filter_bank.ApplyFilter(k10Samples.GetColumn(0).ToValArray()); + filter_bank.ApplyFilter(k10Samples); - ASSERT_EQ(k10Samples.NumElements(), filtered_signal.NumCols()); - ASSERT_EQ(kNumBands, filtered_signal.NumRows()); + ASSERT_EQ(k10Samples.size(), filtered_signal.cols()); + ASSERT_EQ(kNumBands, filtered_signal.rows()); } } // namespace