Skip to content

Restructure gammatone implementation using Eigen #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
terpste wants to merge 1 commit into
base: eigen-replace-arma
Choose a base branch
from
Open

Restructure gammatone implementation using Eigen #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
144 changes: 88 additions & 56 deletions src/gammatone_filterbank.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,81 +16,113 @@

#include <utility>
#include <valarray>
#include <cstring>

#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<double>& 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<double> &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<double> GammatoneFilterBank::ApplyFilter(
const std::valarray<double>& signal) {
AMatrix<double> output(num_bands_, signal.size());
std::valarray<double> 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
43 changes: 24 additions & 19 deletions src/gammatone_spectrogram_builder.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -63,32 +67,33 @@ absl::StatusOr<Spectrogram> GammatoneSpectrogramBuilder::Build(
" required minimum)."));
}
size_t num_cols = 1 + floor((sig.NumRows() - window.size) / hop_size);
AMatrix<double> 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<double> 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<double> 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<double> 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<double> out_matrix(out_eigen);

// Order the center freq bands from lowest to highest.
std::vector<double> ordered_cfb;
ordered_cfb.reserve(erb_rslt.centerFreqs.size());
Expand Down
17 changes: 17 additions & 0 deletions src/include/equivalent_rectangular_bandwidth.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*/
Expand Down
33 changes: 18 additions & 15 deletions src/include/gammatone_filterbank.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@

#include "amatrix.h"

#include "Eigen/Dense"

namespace Visqol {

/**
Expand Down Expand Up @@ -66,7 +68,7 @@ class GammatoneFilterBank {
*
* @return The filtered output.
*/
AMatrix<double> ApplyFilter(const std::valarray<double>& signal);
Eigen::MatrixXd ApplyFilter(const Eigen::MatrixXd &signal);

/**
* Set the equivalent rectangular bandwidth filter coefficients that are to
Expand Down Expand Up @@ -95,20 +97,21 @@ class GammatoneFilterBank {
*/
double min_freq_;

std::valarray<std::valarray<double>> fltr_cond_1_;
std::valarray<std::valarray<double>> fltr_cond_2_;
std::valarray<std::valarray<double>> fltr_cond_3_;
std::valarray<std::valarray<double>> fltr_cond_4_;
std::valarray<double> fltr_coeff_A0_;
std::valarray<double> fltr_coeff_A11_;
std::valarray<double> fltr_coeff_A12_;
std::valarray<double> fltr_coeff_A13_;
std::valarray<double> fltr_coeff_A14_;
std::valarray<double> fltr_coeff_A2_;
std::valarray<double> fltr_coeff_B0_;
std::valarray<double> fltr_coeff_B1_;
std::valarray<double> fltr_coeff_B2_;
std::valarray<double> 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

Expand Down
14 changes: 7 additions & 7 deletions tests/gammatone_filterbank_test.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,24 +25,24 @@ 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<double> k10Samples{
std::valarray<double>{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);
AMatrix<double> filter_coeffs = AMatrix<double>(erb.filterCoeffs);
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
Expand Down