From 80d6f85a18e8d35250a67340ecc414b62267391e Mon Sep 17 00:00:00 2001 From: Stergios Terpinas Date: Thu, 15 Jun 2023 13:01:57 +0300 Subject: [PATCH 1/2] Replace Armadillo with Eigen. --- BUILD | 2 +- WORKSPACE | 19 ++-- src/amatrix.cc | 230 ++++++++++++++++++++++++------------------ src/include/amatrix.h | 19 ++-- 4 files changed, 152 insertions(+), 118 deletions(-) diff --git a/BUILD b/BUILD index 4b1f358..76d660e 100644 --- a/BUILD +++ b/BUILD @@ -129,7 +129,7 @@ cc_library( ":similarity_to_quality_mapper", ":tflite_quality_mapper", ":visqol_config_cc_proto", - "@armadillo_headers//:armadillo_header", + "@eigen//:eigen", "@com_google_absl//absl/flags:flag", "@com_google_absl//absl/flags:parse", "@com_google_absl//absl/flags:usage", diff --git a/WORKSPACE b/WORKSPACE index f3c5aeb..7e0019c 100644 --- a/WORKSPACE +++ b/WORKSPACE @@ -158,20 +158,21 @@ cc_library( urls = ["https://github.com/cjlin1/libsvm/archive/v324.zip"], ) -# Armadillo Headers +# Eigen dependency http_archive( - name = "armadillo_headers", + name = "eigen", build_file_content = """ cc_library( - name = "armadillo_header", - hdrs = glob(["include/armadillo", "include/armadillo_bits/*.hpp"]), - includes = ["include/"], - visibility = ["//visibility:public"], + name = 'eigen', + srcs = [], + includes = ['Eigen'], + hdrs = glob(['Eigen/**']), + visibility = ['//visibility:public'], ) """, - sha256 = "53d7ad6124d06fdede8d839c091c649c794dae204666f1be0d30d7931737d635", - strip_prefix = "armadillo-9.900.1", - urls = ["http://sourceforge.net/projects/arma/files/armadillo-9.900.1.tar.xz"], + sha256 = "0fa5cafe78f66d2b501b43016858070d52ba47bd9b1016b0165a7b8e04675677", + strip_prefix = "eigen-3.3.9", + urls = ["https://gitlab.com/libeigen/eigen/-/archive/3.3.9/eigen-3.3.9.tar.bz2",], ) # PFFFT diff --git a/src/amatrix.cc b/src/amatrix.cc index 5697821..83a25fb 100644 --- a/src/amatrix.cc +++ b/src/amatrix.cc @@ -14,8 +14,8 @@ #include "amatrix.h" -#include #include +#include #include #include #include @@ -25,95 +25,97 @@ namespace Visqol { template inline AMatrix::AMatrix(const AMatrix& other) { - matrix_ = arma::Mat(other.matrix_); + matrix_ = BackendMatrix(other.matrix_); } template -inline AMatrix::AMatrix(const arma::Mat& mat) { +inline AMatrix::AMatrix(const BackendMatrix& mat) { matrix_ = mat; } template inline AMatrix::AMatrix(const std::vector& col) { - matrix_ = arma::Mat(col); + matrix_ = BackendMatrix::Map(col.data(), col.size(), 1); } template inline AMatrix::AMatrix(const absl::Span& col) { - matrix_ = arma::Mat(col.data(), col.size(), 1); + matrix_ = BackendMatrix::Map(col.data(), col.size(), 1); } template inline AMatrix::AMatrix(const std::valarray& va) { std::vector v; - v.reserve(va.size()); + v.resize(va.size()); v.assign(std::begin(va), std::end(va)); - matrix_ = arma::Mat(v); + matrix_ = BackendMatrix::Map(v.data(), v.size(), 1); } template inline AMatrix::AMatrix(const std::vector>& vecOfCols) { // assumes all vectors all of equal length - matrix_ = arma::Mat(vecOfCols[0].size(), vecOfCols.size()); + matrix_ = BackendMatrix(vecOfCols[0].size(), vecOfCols.size()); for (size_t i = 0; i < vecOfCols.size(); i++) { - matrix_.col(i) = arma::Col(vecOfCols[i]); + matrix_.col(i) = BackendMatrix::Map(vecOfCols[i].data(), vecOfCols[i].size(), 1); } } template inline AMatrix::AMatrix(size_t rows, size_t cols, std::vector&& data) { - matrix_ = arma::Mat(&data[0], (arma::uword)rows, (arma::uword)cols, false); + // TODO: We should avoid copying here + matrix_ = BackendMatrix::Map(data.data(), rows, cols); } template inline AMatrix::AMatrix(size_t rows, size_t cols, const std::vector& data) { - matrix_ = arma::Mat(&data[0], (arma::uword)rows, (arma::uword)cols); + matrix_ = BackendMatrix::Map(data.data(), rows, cols); } template inline AMatrix::AMatrix(size_t rows, size_t cols) { - matrix_ = arma::Mat((arma::uword)rows, (arma::uword)cols); + matrix_ = BackendMatrix(rows, cols); } template -inline AMatrix::AMatrix(arma::Mat&& matrix) { +inline AMatrix::AMatrix(BackendMatrix&& matrix) { matrix_ = std::move(matrix); } template AMatrix AMatrix::GetSpan(size_t rowStart, size_t rowEnd, size_t colStart, size_t colEnd) const { - arma::Mat m = matrix_.submat(rowStart, colStart, rowEnd, colEnd); + size_t num_rows = rowEnd - rowStart + 1; + size_t num_cols = colEnd - colStart + 1; + BackendMatrix m = matrix_.block(rowStart, colStart, num_rows, num_cols); return AMatrix(m); } template inline AMatrix AMatrix::Filled(size_t rows, size_t cols, T initialValue) { - arma::Mat m((arma::uword)rows, (arma::uword)cols); + BackendMatrix m (rows, cols); m.fill(initialValue); - AMatrix am(std::move(m)); - return am; + return AMatrix(std::move(m)); } template inline T& AMatrix::operator()(size_t row, size_t column) { - return matrix_((arma::uword)row, (arma::uword)column); + return matrix_(row, column); } template inline T AMatrix::operator()(size_t row, size_t column) const { - return matrix_((arma::uword)row, (arma::uword)column); + return matrix_(row, column); } template inline T& AMatrix::operator()(size_t elementIndex) { - return matrix_.at((arma::uword)elementIndex); + return *(matrix_.data() + elementIndex); } template inline T AMatrix::operator()(size_t elementIndex) const { - return matrix_((arma::uword)elementIndex); + return *(matrix_.data() + elementIndex); } template @@ -124,281 +126,309 @@ inline AMatrix& AMatrix::operator=(const AMatrix& other) { template inline bool AMatrix::operator==(const AMatrix& other) const { - if (matrix_.n_rows != other.matrix_.n_rows) return false; - if (matrix_.n_cols != other.matrix_.n_cols) return true; - return std::equal(matrix_.cbegin(), matrix_.cend(), other.matrix_.cbegin()); - // Why doesn't == work when documentation says it does? matrix_ == - // other.matrix_; + if (matrix_.rows() != other.matrix_.rows()) return false; + if (matrix_.cols() != other.matrix_.cols()) return false; + return matrix_ == other.matrix_; } template inline AMatrix AMatrix::operator+(const AMatrix& other) const { - return AMatrix(matrix_ + other.matrix_); + return AMatrix((matrix_.array() + other.matrix_.array()).matrix()); } template inline AMatrix AMatrix::operator+(T v) const { - return AMatrix(matrix_ + v); + return AMatrix(matrix_.array() + v); } template inline AMatrix AMatrix::operator*(T v) const { - return AMatrix(matrix_ * v); + return AMatrix(matrix_.array() * v); } template inline AMatrix AMatrix::operator-(T v) const { - return AMatrix(matrix_ - v); + return AMatrix(matrix_.array() - v); } template inline AMatrix AMatrix::operator/(T v) const { - return AMatrix(matrix_ / v); + return AMatrix(matrix_.array() / v); } template inline AMatrix AMatrix::operator-(const AMatrix& m) const { - return AMatrix(matrix_ - m.matrix_); + return AMatrix((matrix_.array() - m.matrix_.array()).matrix()); } template inline AMatrix AMatrix::PointWiseProduct(const AMatrix& m) const { - return AMatrix(std::move(matrix_ % m.matrix_)); + return AMatrix(std::move(matrix_.cwiseProduct(m.matrix_))); } template inline AMatrix AMatrix::PointWiseDivide(const AMatrix& m) const { - return AMatrix(std::move(matrix_ / m.matrix_)); + return AMatrix(std::move(matrix_.array() / m.matrix_.array())); } template inline AMatrix AMatrix::Transpose() const { - return AMatrix{std::move(trans(matrix_))}; + return AMatrix{std::move(matrix_.transpose())}; } template inline std::vector AMatrix::RowSubset(size_t rowIndex, size_t startColumnIndex, size_t endColumnIndex) const { - std::vector vec(arma::conv_to>::from( - matrix_.row(rowIndex).subvec(startColumnIndex, endColumnIndex))); + std::vector vec; + vec.resize(endColumnIndex - startColumnIndex + 1); + for (size_t colIndex = startColumnIndex; colIndex <= endColumnIndex; ++colIndex) { + vec[colIndex - startColumnIndex] = matrix_(rowIndex, colIndex); + } return vec; } template inline size_t AMatrix::LongestDimensionLength() const { - return (matrix_.n_rows > matrix_.n_cols) ? matrix_.n_rows : matrix_.n_cols; + return (matrix_.rows() > matrix_.cols()) ? matrix_.rows() : matrix_.cols(); } template <> inline AMatrix AMatrix::Abs() const { - AMatrix absMatrix(matrix_.n_rows, matrix_.n_cols); - for (size_t i = 0; i < absMatrix.NumElements(); i++) { - absMatrix.matrix_(i) = std::abs(matrix_(i)); - } - return absMatrix; + return AMatrix(matrix_.cwiseAbs()); } template <> inline AMatrix AMatrix>::Abs() const { - arma::Mat absMatrix; - absMatrix.set_size(matrix_.n_rows, matrix_.n_cols); - for (size_t i = 0; i < absMatrix.n_elem; i++) { - absMatrix(i) = std::abs(matrix_(i)); + BackendMatrix absMatrix (matrix_.rows(), matrix_.cols()); + for (size_t rowIndex = 0; rowIndex < matrix_.rows(); rowIndex++) { + for (size_t colIndex = 0; colIndex < matrix_.rows(); colIndex++) { + absMatrix(rowIndex, colIndex) = std::abs(matrix_(rowIndex, colIndex)); + } } return AMatrix(std::move(absMatrix)); } template inline std::vector AMatrix::GetRow(size_t row) const { - std::vector vec(arma::conv_to>::from(matrix_.row(row))); + std::vector vec; + vec.resize(matrix_.cols()); + for (size_t i = 0; i < matrix_.cols(); ++i) + { + vec[i] = matrix_(row, i); + } return vec; } template inline AMatrix AMatrix::GetRows(size_t rowStart, size_t rowEnd) const { - return AMatrix(std::move(matrix_.rows(rowStart, rowEnd))); + return AMatrix(std::move(matrix_.middleRows(rowStart, rowEnd - rowStart + 1))); } template inline AMatrix AMatrix::GetColumn(size_t column) const { - return AMatrix(matrix_.cols(column, column)); + return AMatrix(matrix_.col(column)); } template inline AMatrix AMatrix::GetColumns(size_t colStart, size_t colEnd) const { - return AMatrix(matrix_.cols(colStart, colEnd)); + return AMatrix(matrix_.middleCols(colStart, colEnd - colStart + 1)); } template inline void AMatrix::SetRow(size_t rowIndex, const std::vector& row) { - matrix_.row(rowIndex) = std::move(arma::Row(row)); + matrix_.row(rowIndex) = BackendMatrix::Map(row.data(), 1, row.size()); } template inline void AMatrix::SetColumn(size_t colIndex, AMatrix&& col) { - matrix_.col(colIndex) = std::move(col.GetArmaMat()); + matrix_.col(colIndex) = std::move(col.matrix_); } template inline void AMatrix::SetColumn(size_t colIndex, std::vector&& col) { - matrix_.col(colIndex) = std::move(arma::Col(col)); + matrix_.col(colIndex) = std::move(BackendMatrix::Map(col.data(), col.size(), 1)); } template inline void AMatrix::SetColumn(size_t colIndex, std::valarray&& col) { - std::vector v; - v.reserve(col.size()); - v.assign(std::begin(col), std::end(col)); - matrix_.col(colIndex) = std::move(arma::Col(v)); + for (size_t i = 0; i < col.size(); ++i) { + matrix_(i, colIndex) = col[i]; + } } template inline void AMatrix::SetRow(size_t rowIndex, std::valarray&& row) { - std::vector v; - v.reserve(row.size()); - v.assign(std::begin(row), std::end(row)); - matrix_.row(rowIndex) = std::move(arma::Row(v)); + for (size_t colIndex = 0; colIndex < matrix_.cols(); ++colIndex) { + matrix_(rowIndex, colIndex) = row[colIndex]; + } } template void AMatrix::Print() const { printf("\n"); - matrix_.print(); + std::cout << matrix_ << std::endl; } template <> void AMatrix>::PrintSummary(const char* s) const { printf("%s\n", s); - size_t outMatSize = matrix_.n_rows; + size_t outMatSize = matrix_.rows(); printf("first five\n"); for (size_t i = 0; i < 5; i++) { - printf("complex[%2zu] = %9.20f , %9.20f\n", i, matrix_(i).real(), - matrix_(i).imag()); + printf("complex[%2zu] = %9.20f , %9.20f\n", i, (matrix_.data() + i)->real(), + (matrix_.data() + i)->imag()); } printf("middle \n"); for (size_t i = (outMatSize / 2) - 4; i < (outMatSize / 2) + 6; i++) { - printf("complex[%2zu] = %9.20f , %9.20f\n", i, matrix_(i).real(), - matrix_(i).imag()); + printf("complex[%2zu] = %9.20f , %9.20f\n", i, (matrix_.data() + i)->real(), + (matrix_.data() + i)->imag()); } printf("last five\n"); for (size_t i = outMatSize - 6; i < outMatSize; i++) { - printf("complex[%2zu] = %9.20f , %9.20f\n", i, matrix_(i).real(), - matrix_(i).imag()); + printf("complex[%2zu] = %9.20f , %9.20f\n", i, (matrix_.data() + i)->real(), + (matrix_.data() + i)->imag()); } } template <> void AMatrix::PrintSummary(const char* s) const { printf("%s\n", s); - size_t outMatSize = matrix_.n_rows; + size_t outMatSize = matrix_.rows(); printf("first five\n"); for (size_t i = 0; i < 5; i++) { - printf("double[%2zu] = %9.20f\n", i, matrix_(i)); + printf("double[%2zu] = %9.20f\n", i, *(matrix_.data() + i)); } printf("middle \n"); for (size_t i = (outMatSize / 2) - 4; i < (outMatSize / 2) + 6; i++) { - printf("double[%2zu] = %9.20f\n", i, matrix_(i)); + printf("double[%2zu] = %9.20f\n", i, *(matrix_.data() + i)); } printf("last five\n"); for (size_t i = outMatSize - 6; i < outMatSize; i++) { - printf("double[%2zu] = %9.20f\n", i, matrix_(i)); + printf("double[%2zu] = %9.20f\n", i, *(matrix_.data() + i)); } } template inline const T* AMatrix::MemPtr() const { - return matrix_.memptr(); + return matrix_.data(); } template inline void AMatrix::Resize(const size_t rows, const size_t cols) { - matrix_.resize(rows, cols); + // We need to try to keep the same data + matrix_.conservativeResize(rows, cols); } template inline AMatrix AMatrix::JoinVertically(const AMatrix& other) const { - // Having a separate line for return is needed to cast to Mat - arma::Mat m = arma::join_vert(matrix_, other.GetArmaMat()); + // Asserting other matrix has the same number of columns + if (matrix_.cols() != other.matrix_.cols()) return AMatrix (); + size_t numCols = matrix_.cols(); + size_t numRows = matrix_.rows() + other.matrix_.rows(); + BackendMatrix m (numRows, numCols); + m.topRows(matrix_.rows()) = matrix_; + m.bottomRows(other.matrix_.rows()) = other.matrix_; return m; } template inline AMatrix AMatrix::FlipUpDown() const { - return AMatrix(arma::flipud(matrix_)); + return AMatrix(matrix_.colwise().reverse()); } template inline AMatrix AMatrix::Mean(kDimension dim) const { - return AMatrix(arma::mean(matrix_, static_cast(dim))); + if (dim == kDimension::COLUMN) { + return AMatrix(matrix_.colwise().mean()); + } + else { // if (dim == kDimension::ROW) + return AMatrix(matrix_.rowwise().mean()); + } } template inline AMatrix AMatrix::StdDev(kDimension dim) const { // The second parameter of arma::stddev is norm_type. // norm_type=0 means that stddev should use the unbiased estimate. - return AMatrix(arma::stddev(matrix_, 0, static_cast(dim))); + if (dim == kDimension::COLUMN) { + if (matrix_.rows() == 1) { + return AMatrix(BackendMatrix::Zero(1, matrix_.cols())); + } + else { + return AMatrix((matrix_.rowwise() - matrix_.colwise().mean()).colwise().norm()/std::sqrt((double)matrix_.rows()-1.0)); + } + } + else { // if (dim == kDimension::ROW) + if (matrix_.cols() == 1) { + return AMatrix(BackendMatrix::Zero(matrix_.rows(), 1)); + } + else { + return AMatrix((matrix_.colwise() - matrix_.rowwise().mean()).rowwise().norm()/std::sqrt((double)matrix_.cols()-1.0)); + } + } } template -inline const arma::Mat& AMatrix::GetArmaMat() const { +inline const BackendMatrix& AMatrix::GetBackendMat() const { return matrix_; } template inline std::vector AMatrix::ToVector() const { - return arma::conv_to>::from(matrix_.col(0)); + std::vector v (static_cast(matrix_.rows())); + for (size_t rowIndex = 0; rowIndex < matrix_.rows(); ++rowIndex) { + v[rowIndex] = matrix_(rowIndex, 0); + } + return std::move(v); } template inline std::valarray AMatrix::ToValArray() const { - std::vector v = arma::conv_to>::from(matrix_.col(0)); - return std::valarray(v.data(), v.size()); -} - -template <> -inline std::vector> -AMatrix>::ToVector() const { - return arma::conv_to>>::from(matrix_.col(0)); + std::valarray v (static_cast(matrix_.rows())); + for (size_t rowIndex = 0; rowIndex < matrix_.rows(); ++rowIndex) { + v[rowIndex] = matrix_(rowIndex, 0); + } + return std::move(v); } template inline size_t AMatrix::NumRows() const { - return matrix_.n_rows; + return matrix_.rows(); } template inline size_t AMatrix::NumCols() const { - return matrix_.n_cols; + return matrix_.cols(); } template inline size_t AMatrix::NumElements() const { - return matrix_.n_elem; + return matrix_.size(); } template inline typename AMatrix::iterator AMatrix::begin() { - return matrix_.begin(); + return matrix_.data(); } template inline typename AMatrix::iterator AMatrix::end() { - return matrix_.end(); + return matrix_.data() + matrix_.size(); } template inline typename AMatrix::const_iterator AMatrix::cbegin() const { - return matrix_.cbegin(); + return matrix_.data(); } template inline typename AMatrix::const_iterator AMatrix::cend() const { - return matrix_.cend(); + return matrix_.data() + matrix_.size(); } template inline const T* AMatrix::data() const { - return matrix_.mem; + return matrix_.data(); } template inline T* AMatrix::mutData() const { - return const_cast(matrix_.mem); + return const_cast(matrix_.data()); } template class AMatrix; diff --git a/src/include/amatrix.h b/src/include/amatrix.h index 0838f8b..341c821 100644 --- a/src/include/amatrix.h +++ b/src/include/amatrix.h @@ -17,13 +17,12 @@ #ifndef VISQOL_INCLUDE_AMATRIX_H #define VISQOL_INCLUDE_AMATRIX_H -#include +#include "Eigen/Dense" #include #include #include #include "absl/types/span.h" -#define ARMA_64BIT_WORD #if defined(_MSC_VER) && _MSC_VER >= 1400 #pragma warning(push) @@ -33,6 +32,10 @@ namespace Visqol { enum class kDimension { COLUMN = 0, ROW = 1 }; +// Define an alias for the Eigen matrix +template +using BackendMatrix = Eigen::Matrix; + template class AMatrix; template @@ -43,7 +46,7 @@ template class AMatrix { public: AMatrix() {} - AMatrix(const arma::Mat& mat); + AMatrix(const BackendMatrix& mat); AMatrix(const AMatrix& other); AMatrix(const std::vector& col); AMatrix(const absl::Span& col); @@ -52,7 +55,7 @@ class AMatrix { AMatrix(size_t rows, size_t cols); AMatrix(size_t rows, size_t cols, std::vector&& data); AMatrix(size_t rows, size_t cols, const std::vector& data); - AMatrix(arma::Mat&& matrix); // could this be private? + AMatrix(BackendMatrix&& matrix); // could this be private? T& operator()(size_t row, size_t column); T operator()(size_t row, size_t column) const; T& operator()(size_t elementIndex); @@ -104,11 +107,11 @@ class AMatrix { std::valarray ToValArray() const; // hack to access private member - const arma::Mat& GetArmaMat() const; + const BackendMatrix& GetBackendMat() const; public: - typedef typename arma::Mat::iterator iterator; - typedef typename arma::Mat::const_iterator const_iterator; + using iterator = T*; + using const_iterator = const T*; iterator begin(); iterator end(); const_iterator cbegin() const; @@ -118,7 +121,7 @@ class AMatrix { T* mutData() const; // bit of a hack private: - arma::Mat matrix_; + BackendMatrix matrix_; }; } // namespace Visqol From a475fa6faff47c2342548967c124daca8605bbdf Mon Sep 17 00:00:00 2001 From: Stergios Terpinas Date: Thu, 15 Jun 2023 13:02:51 +0300 Subject: [PATCH 2/2] 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 8683226..44fff6c 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 6314f8e..c539c81 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 fcf01c8..41c3a05 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 4102bff..4ac75d1 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 1ab8ae0..0940f4b 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