Mercurial > hg > aimc
changeset 679:594b410c2aed
Sixth revision of Alex Brandmeyer's C++ implementation. Only small changes in response to Lyon's comments from r285. Note: I tried to use a consistent indentation with two spaces, but also preserving parenthetical structure to make reading the longer equations easier. Please advise if this is OK. Additional documentation and tests with non-standard parameter sets will be added in a subsequent revision.
| author | alexbrandmeyer |
|---|---|
| date | Tue, 28 May 2013 15:54:54 +0000 |
| parents | 7f424c1a8b78 |
| children | be55786eeb04 |
| files | trunk/carfac/SConstruct trunk/carfac/agc_params.cc trunk/carfac/agc_params.h trunk/carfac/car_params.cc trunk/carfac/car_params.h trunk/carfac/carfac.cc trunk/carfac/carfac.h trunk/carfac/carfac_common.h trunk/carfac/carfac_output.cc trunk/carfac/carfac_output.h trunk/carfac/carfac_test.cc trunk/carfac/ear.cc trunk/carfac/ear.h trunk/carfac/ihc_params.cc trunk/carfac/ihc_params.h |
| diffstat | 15 files changed, 204 insertions(+), 315 deletions(-) [+] |
line wrap: on
line diff
--- a/trunk/carfac/SConstruct Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/SConstruct Tue May 28 15:54:54 2013 +0000 @@ -43,14 +43,14 @@ import os carfac_sources = [ - 'carfac_common.cc', - 'agc_params.cc', + 'carfac_common.cc', + 'agc_params.h', 'agc_coeffs.h', 'agc_state.h', 'carfac_output.cc', - 'car_params.cc', + 'car_params.h', 'car_coeffs.h', - 'ihc_params.cc', + 'ihc_params.h', 'car_state.h', 'ihc_coeffs.h', 'ihc_state.h',
--- a/trunk/carfac/agc_params.cc Mon May 27 16:36:54 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,43 +0,0 @@ -// -// agc_params.cc -// CARFAC Open Source C++ Library -// -// Created by Alex Brandmeyer on 5/10/13. -// -// This C++ file is part of an implementation of Lyon's cochlear model: -// "Cascade of Asymmetric Resonators with Fast-Acting Compression" -// to supplement Lyon's upcoming book "Human and Machine Hearing" -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -#include "agc_params.h" - -// The default constructor for AGCParams initializes with the settings from -// Lyon's book 'Human and Machine Hearing' -AGCParams::AGCParams() { - n_stages_ = 4; - agc_stage_gain_ = 2.0; - time_constants_.resize(n_stages_); - agc1_scales_.resize(n_stages_); - agc2_scales_.resize(n_stages_); - agc1_scales_[0] = 1.0; - agc2_scales_[0] = 1.65; - time_constants_[0] = 0.002; - for (int i = 1; i < n_stages_; ++i) { - agc1_scales_[i] = agc1_scales_[i - 1] * sqrt(2.0); - agc2_scales_[i] = agc2_scales_[i - 1] * sqrt(2.0); - time_constants_[i] = time_constants_[i - 1] * 4.0; - } - decimation_ = {8, 2, 2, 2}; - agc_mix_coeff_ = 0.5; -} \ No newline at end of file
--- a/trunk/carfac/agc_params.h Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/agc_params.h Tue May 28 15:54:54 2013 +0000 @@ -26,7 +26,23 @@ #include "carfac_common.h" struct AGCParams { - AGCParams(); + AGCParams() { + n_stages_ = 4; + agc_stage_gain_ = 2.0; + time_constants_.resize(n_stages_); + agc1_scales_.resize(n_stages_); + agc2_scales_.resize(n_stages_); + agc1_scales_[0] = 1.0; + agc2_scales_[0] = 1.65; + time_constants_[0] = 0.002; + for (int i = 1; i < n_stages_; ++i) { + agc1_scales_[i] = agc1_scales_[i - 1] * sqrt(2.0); + agc2_scales_[i] = agc2_scales_[i - 1] * sqrt(2.0); + time_constants_[i] = time_constants_[i - 1] * 4.0; + } + decimation_ = {8, 2, 2, 2}; + agc_mix_coeff_ = 0.5; + } int n_stages_; FPType agc_stage_gain_; FPType agc_mix_coeff_;
--- a/trunk/carfac/car_params.cc Mon May 27 16:36:54 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,38 +0,0 @@ -// -// car_params.cc -// CARFAC Open Source C++ Library -// -// Created by Alex Brandmeyer on 5/10/13. -// -// This C++ file is part of an implementation of Lyon's cochlear model: -// "Cascade of Asymmetric Resonators with Fast-Acting Compression" -// to supplement Lyon's upcoming book "Human and Machine Hearing" -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -#include "car_params.h" - -CARParams::CARParams() { - velocity_scale_ = 0.1; - v_offset_ = 0.04; - min_zeta_ = 0.1; - max_zeta_ = 0.35; - first_pole_theta_ = 0.85 * PI; - zero_ratio_ = sqrt(2.0); - high_f_damping_compression_ = 0.5; - erb_per_step_ = 0.5; - min_pole_hz_ = 30; - erb_break_freq_ = 165.3; // This is the Greenwood map's break frequency. - // This represents Glassberg and Moore's high-cf ratio. - erb_q_ = 1000/(24.7*4.37); -}; \ No newline at end of file
--- a/trunk/carfac/car_params.h Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/car_params.h Tue May 28 15:54:54 2013 +0000 @@ -26,7 +26,21 @@ #include "carfac_common.h" struct CARParams { - CARParams(); // The constructor initializes using default parameter values. + // The constructor initializes using default parameter values. + CARParams() { + velocity_scale_ = 0.1; + v_offset_ = 0.04; + min_zeta_ = 0.1; + max_zeta_ = 0.35; + first_pole_theta_ = 0.85 * kPi; + zero_ratio_ = sqrt(2.0); + high_f_damping_compression_ = 0.5; + erb_per_step_ = 0.5; + min_pole_hz_ = 30; + erb_break_freq_ = 165.3; // This is the Greenwood map's break frequency. + // This represents Glassberg and Moore's high-cf ratio. + erb_q_ = 1000/(24.7*4.37); + }; FPType velocity_scale_; // This is used for the velocity nonlinearity. FPType v_offset_; // The offset gives us quadratic part. FPType min_zeta_; // This is the minimum damping factor in mid-freq channels.
--- a/trunk/carfac/carfac.cc Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/carfac.cc Tue May 28 15:54:54 2013 +0000 @@ -30,14 +30,14 @@ fs_ = fs; ears_.resize(n_ears_); n_ch_ = 0; - FPType pole_hz = car_params.first_pole_theta_ * fs / (2 * PI); + FPType pole_hz = car_params.first_pole_theta_ * fs / (2 * kPi); while (pole_hz > car_params.min_pole_hz_) { ++n_ch_; pole_hz = pole_hz - car_params.erb_per_step_ * ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_); } pole_freqs_.resize(n_ch_); - pole_hz = car_params.first_pole_theta_ * fs / (2 * PI); + pole_hz = car_params.first_pole_theta_ * fs / (2 * kPi); for (int ch = 0; ch < n_ch_; ++ch) { pole_freqs_(ch) = pole_hz; pole_hz = pole_hz - car_params.erb_per_step_ * @@ -53,7 +53,7 @@ DesignAGCCoeffs(agc_params, fs, &agc_coeffs); // Once we have the coefficient structure we can design the ears. for (auto& ear : ears_) { - ear.DesignEar(n_ch_, fs_, car_coeffs, ihc_coeffs, + ear.InitEar(n_ch_, fs_, car_coeffs, ihc_coeffs, agc_coeffs); } } @@ -87,8 +87,6 @@ const bool open_loop, CARFACOutput* seg_output) { // A nested loop structure is used to iterate through the individual samples // for each ear (audio channel). - FloatArray car_out(n_ch_); - FloatArray ihc_out(n_ch_); bool updated; // This variable is used by the AGC stage. for (int32_t i = 0; i < length; ++i) { for (int j = 0; j < n_ears_; ++j) { @@ -98,11 +96,11 @@ FPType input = sound_data[j][start+i]; // Now we apply the three stages of the model in sequence to the current // audio sample. - ear.CARStep(input, &car_out); - ear.IHCStep(car_out, &ihc_out); - updated = ear.AGCStep(ihc_out); + ear.CARStep(input); + ear.IHCStep(ear.car_out()); + updated = ear.AGCStep(ear.ihc_out()); } - seg_output->StoreOutput(&ears_); + seg_output->StoreOutput(ears_); if (updated) { if (n_ears_ > 1) { CrossCouple(); @@ -161,12 +159,12 @@ car_coeffs->h_coeffs_.resize(n_ch_); car_coeffs->g0_coeffs_.resize(n_ch_); FPType f = car_params.zero_ratio_ * car_params.zero_ratio_ - 1.0; - FloatArray theta = pole_freqs * ((2.0 * PI) / fs); + FloatArray theta = pole_freqs * ((2.0 * kPi) / fs); car_coeffs->c0_coeffs_ = theta.sin(); car_coeffs->a0_coeffs_ = theta.cos(); FPType ff = car_params.high_f_damping_compression_; - FloatArray x = theta / PI; - car_coeffs->zr_coeffs_ = PI * (x - (ff * (x*x*x))); + FloatArray x = theta / kPi; + car_coeffs->zr_coeffs_ = kPi * (x - (ff * (x*x*x))); FPType max_zeta = car_params.max_zeta_; FPType min_zeta = car_params.min_zeta_; car_coeffs->r1_coeffs_ = (1.0 - (car_coeffs->zr_coeffs_ * max_zeta)); @@ -183,9 +181,8 @@ FloatArray r = car_coeffs->r1_coeffs_ + (car_coeffs->zr_coeffs_ * relative_undamping); car_coeffs->g0_coeffs_ = (1.0 - (2.0 * r * car_coeffs->a0_coeffs_) + (r*r)) / - (1 - (2 * r * car_coeffs->a0_coeffs_) + - (car_coeffs->h_coeffs_ * r * car_coeffs->c0_coeffs_) - + (r*r)); + (1 - (2 * r * car_coeffs->a0_coeffs_) + + (car_coeffs->h_coeffs_ * r * car_coeffs->c0_coeffs_) + (r*r)); } void CARFAC::DesignIHCCoeffs(const IHCParams& ihc_params, const FPType fs, @@ -242,7 +239,7 @@ ihc_coeffs->rest_cap2_ = ihc_coeffs->cap2_voltage_; } } - ihc_coeffs->ac_coeff_ = 2 * PI * ihc_params.ac_corner_hz_ / fs; + ihc_coeffs->ac_coeff_ = 2 * kPi * ihc_params.ac_corner_hz_ / fs; } void CARFAC::DesignAGCCoeffs(const AGCParams& agc_params, const FPType fs,
--- a/trunk/carfac/carfac.h Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/carfac.h Tue May 28 15:54:54 2013 +0000 @@ -78,4 +78,4 @@ FloatArray pole_freqs_; }; -#endif +#endif \ No newline at end of file
--- a/trunk/carfac/carfac_common.h Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/carfac_common.h Tue May 28 15:54:54 2013 +0000 @@ -51,32 +51,27 @@ // occur. // <math.h> is used during coefficient calculations and runtime operations. #include <math.h> -// <vector> is used in place of 2d Eigen Arrays for the AGC memory +//<vector> is used to store two dimensional data structures. #include <vector> -// The Eigen library is used extensively for 1d and 2d floating point arrays. +// The Eigen library is used extensively for floating point arrays. // For more information, see: http://eigen.tuxfamily.org #include <Eigen/Dense> + using namespace Eigen; -// One constant value is defined here, but see my TODO regarding style issues. -// A fixed value of PI is defined throughout the project. -// TODO alexbrandmeyer: verify that this is OK with Google Style. -#define PI 3.141592653589793238 - -// Three typedefs are used for enabling quick switching of precision and array -// usage. These names may be changed based on the discussions being had within -// the AIMC group. // The 'FPType' typedef is used to enable easy switching in precision level. // It's currently set to double for during the unit testing phase of the // project. typedef double FPType; // A typedef is used to define a one-dimensional Eigen array with the same // precision level as FPType. -typedef Eigen::Array<FPType, Dynamic, 1> FloatArray; +typedef Eigen::Array<FPType, Dynamic, 1> FloatArray; + +// A fixed value of PI is defined throughout the project. +static const FPType kPi = 3.141592653589793238; // Two helper functions are defined here for use by the different model stages // in calculating coeffecients and during model runtime. - // Function: ERBHz // Auditory filter nominal Equivalent Rectangular Bandwidth // Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138 @@ -88,4 +83,4 @@ // values. This is here because it is called both in design and run phases. FloatArray CARFACDetect (const FloatArray& x); -#endif +#endif \ No newline at end of file
--- a/trunk/carfac/carfac_output.cc Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/carfac_output.cc Tue May 28 15:54:54 2013 +0000 @@ -21,6 +21,7 @@ // limitations under the License. #include "carfac_output.h" + using std::vector; void CARFACOutput::Init(const int n_ears, const bool store_nap, @@ -35,29 +36,28 @@ } -void CARFACOutput::StoreOutput(std::vector<Ear>* ears) { +void CARFACOutput::StoreOutput(const vector<Ear>& ears) { if (store_nap_) { nap_.push_back(vector<FloatArray>()); - for (auto& ear : *ears) { + for (auto ear : ears) { nap_.back().push_back(ear.ihc_out()); } } if (store_ohc_) { ohc_.push_back(vector<FloatArray>()); - for (auto& ear : *ears) { + for (auto ear : ears) { ohc_.back().push_back(ear.za_memory()); } } if (store_agc_) { agc_.push_back(vector<FloatArray>()); - for (auto& ear : *ears) { - + for (auto ear : ears) { agc_.back().push_back(ear.zb_memory()); } } if (store_bm_) { bm_.push_back(vector<FloatArray>()); - for (auto& ear : *ears) { + for (auto ear : ears) { bm_.back().push_back(ear.zy_memory()); } }
--- a/trunk/carfac/carfac_output.h Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/carfac_output.h Tue May 28 15:54:54 2013 +0000 @@ -44,17 +44,14 @@ public: void Init(const int n_ears, const bool store_nap, const bool store_nap_decim, const bool store_bm, const bool store_ohc, const bool store_agc); - void StoreOutput(std::vector<Ear>* ears); - FPType nap(const int ear, const int32_t timepoint, const int channel) { - return nap_[timepoint][ear](channel); } - FPType bm(const int ear, const int32_t timepoint, const int channel) { - return bm_[timepoint][ear](channel); } - std::deque<std::vector<FloatArray>> nap_; - std::deque<std::vector<FloatArray>> nap_decim_; - std::deque<std::vector<FloatArray>> bm_; - std::deque<std::vector<FloatArray>> ohc_; - std::deque<std::vector<FloatArray>> agc_; - + void StoreOutput(const std::vector<Ear>& ears); + // Here we define several acessors for the data members. + const std::deque<std::vector<FloatArray>>& nap() { return nap_;} + const std::deque<std::vector<FloatArray>>& bm() { return bm_; } + const std::deque<std::vector<FloatArray>>& nap_decim() { return nap_decim_;} + const std::deque<std::vector<FloatArray>>& ohc() { return ohc_; } + const std::deque<std::vector<FloatArray>>& agc() { return agc_;} + private: int n_ears_; bool store_nap_; @@ -62,6 +59,11 @@ bool store_bm_; bool store_ohc_; bool store_agc_; + std::deque<std::vector<FloatArray>> nap_; + std::deque<std::vector<FloatArray>> nap_decim_; + std::deque<std::vector<FloatArray>> bm_; + std::deque<std::vector<FloatArray>> ohc_; + std::deque<std::vector<FloatArray>> agc_; }; #endif \ No newline at end of file
--- a/trunk/carfac/carfac_test.cc Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/carfac_test.cc Tue May 28 15:54:54 2013 +0000 @@ -25,27 +25,31 @@ // GoogleTest is now included for running unit tests #include <gtest/gtest.h> #include "carfac.h" + using std::vector; using std::string; using std::ifstream; using std::ofstream; -#define TEST_SRC_DIR "./test_data/" +// This is the 'test_data' subdirectory of aimc/carfac that specifies where to +// locate the text files produced by 'CARFAC_GenerateTestData.m' for comparing +// the ouput of the Matlab version of CARFAC with this C++ version. +static const char* kTestSourceDir= "./test_data/"; // Here we specify the level to which the output should match (7 decimals). -#define PRECISION_LEVEL 1.0e-07 +static const float kPrecisionLevel = 1.0e-7; // Three helper functions are defined here for loading the test data generated // by the Matlab version of CARFAC. // This loads one-dimensional FloatArrays from single-column text files. void WriteNAPOutput(CARFACOutput output, const string filename, int ear) { - string fullfile = TEST_SRC_DIR + filename; + string fullfile = kTestSourceDir + filename; ofstream ofile(fullfile.c_str()); - int32_t n_timepoints = output.nap_.size(); - int channels = output.nap_[0][0].size(); + int32_t n_timepoints = output.nap().size(); + int channels = output.nap()[0][0].size(); if (ofile.is_open()) { for (int32_t i = 0; i < n_timepoints; ++i) { for (int j = 0; j < channels; ++j) { - ofile << output.nap(ear, i, j); + ofile << output.nap()[i][ear](j); if ( j < channels - 1) { ofile << " "; } @@ -56,8 +60,8 @@ ofile.close(); } -FloatArray LoadTestData(const std::string filename, const int number_points) { - string fullfile = TEST_SRC_DIR + filename; +FloatArray LoadTestData(const string filename, const int number_points) { + string fullfile = kTestSourceDir + filename; ifstream file(fullfile.c_str()); FPType myarray[number_points]; FloatArray output(number_points); @@ -74,7 +78,7 @@ // This loads a vector of FloatArrays from multi-column text files. vector<FloatArray> Load2dTestData(const string filename, const int rows, const int columns) { - string fullfile = TEST_SRC_DIR + filename; + string fullfile = kTestSourceDir + filename; ifstream file(fullfile.c_str()); FPType myarray[rows][columns]; vector<FloatArray> output; @@ -98,7 +102,7 @@ // Matlab using the wavread() function. vector<vector<float>> Load2dAudioVector(string filename, int timepoints, int channels) { - string fullfile = TEST_SRC_DIR + filename; + string fullfile = kTestSourceDir + filename; ifstream file(fullfile.c_str()); vector<vector<float>> output; output.resize(channels); @@ -120,22 +124,17 @@ int n_timepoints = 882; int n_channels = 71; int n_ears = 2; - std::string filename = "binaural_test_nap1.txt"; - std::vector<FloatArray> nap1 = Load2dTestData(filename, n_timepoints, - n_channels); + string filename = "binaural_test_nap1.txt"; + vector<FloatArray> nap1 = Load2dTestData(filename, n_timepoints, n_channels); filename = "binaural_test_bm1.txt"; - std::vector<FloatArray> bm1 = Load2dTestData(filename, n_timepoints, - n_channels); + vector<FloatArray> bm1 = Load2dTestData(filename, n_timepoints, n_channels); filename = "binaural_test_nap2.txt"; - std::vector<FloatArray> nap2 = Load2dTestData(filename, n_timepoints, - n_channels); + vector<FloatArray> nap2 = Load2dTestData(filename, n_timepoints, n_channels); filename = "binaural_test_bm2.txt"; - std::vector<FloatArray> bm2 = Load2dTestData(filename, n_timepoints, - n_channels); + vector<FloatArray> bm2 = Load2dTestData(filename, n_timepoints, n_channels); filename = "file_signal_binaural_test.txt"; - std::vector<std::vector<float>> sound_data = Load2dAudioVector(filename, - n_timepoints, - n_ears); + vector<vector<float>> sound_data = Load2dAudioVector(filename, n_timepoints, + n_ears); CARParams car_params; IHCParams ihc_params; AGCParams agc_params; @@ -153,23 +152,23 @@ int n_ch = 71; for (int timepoint = 0; timepoint < n_timepoints; ++timepoint) { for (int channel = 0; channel < n_ch; ++channel) { - FPType cplusplus = my_output.nap(ear, timepoint, channel); + FPType cplusplus = my_output.nap()[timepoint][ear](channel); FPType matlab = nap1[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); - cplusplus = my_output.bm(ear, timepoint, channel); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); + cplusplus = my_output.bm()[timepoint][ear](channel); matlab = bm1[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); } } ear = 1; for (int timepoint = 0; timepoint < n_timepoints; ++timepoint) { for (int channel = 0; channel < n_ch; ++channel) { - FPType cplusplus = my_output.nap(ear, timepoint, channel); + FPType cplusplus = my_output.nap()[timepoint][ear](channel); FPType matlab = nap2[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); - cplusplus = my_output.bm(ear, timepoint, channel); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); + cplusplus = my_output.bm()[timepoint][ear](channel); matlab = bm2[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); } } } @@ -179,21 +178,16 @@ int n_channels = 83; int n_ears = 2; string filename = "long_test_nap1.txt"; - vector<FloatArray> nap1 = Load2dTestData(filename, n_timepoints, - n_channels); + vector<FloatArray> nap1 = Load2dTestData(filename, n_timepoints, n_channels); filename = "long_test_bm1.txt"; - vector<FloatArray> bm1 = Load2dTestData(filename, n_timepoints, - n_channels); + vector<FloatArray> bm1 = Load2dTestData(filename, n_timepoints, n_channels); filename = "long_test_nap2.txt"; - vector<FloatArray> nap2 = Load2dTestData(filename, n_timepoints, - n_channels); + vector<FloatArray> nap2 = Load2dTestData(filename, n_timepoints, n_channels); filename = "long_test_bm2.txt"; - vector<FloatArray> bm2 = Load2dTestData(filename, n_timepoints, - n_channels); + vector<FloatArray> bm2 = Load2dTestData(filename, n_timepoints, n_channels); filename = "file_signal_long_test.txt"; - vector<vector<float>> sound_data = Load2dAudioVector(filename, - n_timepoints, - n_ears); + vector<vector<float>> sound_data = Load2dAudioVector(filename, n_timepoints, + n_ears); CARParams car_params; IHCParams ihc_params; AGCParams agc_params; @@ -210,30 +204,23 @@ int ear = 0; for (int timepoint = 0; timepoint < n_timepoints; ++timepoint) { for (int channel = 0; channel < n_channels; ++channel) { - FPType cplusplus = my_output.nap(ear, timepoint, channel); + FPType cplusplus = my_output.nap()[timepoint][ear](channel); FPType matlab = nap1[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); - cplusplus = my_output.bm(ear, timepoint, channel); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); + cplusplus = my_output.bm()[timepoint][ear](channel); matlab = bm1[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); } } ear = 1; for (int timepoint = 0; timepoint < n_timepoints; ++timepoint) { for (int channel = 0; channel < n_channels; ++channel) { - FPType cplusplus = my_output.nap(ear, timepoint, channel); + FPType cplusplus = my_output.nap()[timepoint][ear](channel); FPType matlab = nap2[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); - cplusplus = my_output.bm(ear, timepoint, channel); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); + cplusplus = my_output.bm()[timepoint][ear](channel); matlab = bm2[timepoint](channel); - ASSERT_NEAR(cplusplus, matlab, PRECISION_LEVEL); + ASSERT_NEAR(cplusplus, matlab, kPrecisionLevel); } } -} - -int main(int argc, char **argv) { - // This initializes the GoogleTest unit testing framework. - ::testing::InitGoogleTest(&argc, argv); - // This runs all of the tests that we've defined above. - return RUN_ALL_TESTS(); } \ No newline at end of file
--- a/trunk/carfac/ear.cc Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/ear.cc Tue May 28 15:54:54 2013 +0000 @@ -21,16 +21,14 @@ // limitations under the License. #include <assert.h> - #include "ear.h" // The 'InitEar' function takes a set of model parameters and initializes the // design coefficients and model state variables needed for running the model // on a single audio channel. -void Ear::DesignEar(const int n_ch, const FPType fs, - const CARCoeffs& car_coeffs, - const IHCCoeffs& ihc_coeffs, - const std::vector<AGCCoeffs>& agc_coeffs) { +void Ear::InitEar(const int n_ch, const FPType fs, const CARCoeffs& car_coeffs, + const IHCCoeffs& ihc_coeffs, + const std::vector<AGCCoeffs>& agc_coeffs) { // The first section of code determines the number of channels that will be // used in the model on the basis of the sample rate and the CAR parameters // that have been passed to this function. @@ -81,7 +79,7 @@ } } -void Ear::CARStep(const FPType input, FloatArray* car_out) { +void Ear::CARStep(const FPType input) { // This interpolates g. car_state_.g_memory_ = car_state_.g_memory_ + car_state_.dg_memory_; // This calculates the AGC interpolation state. @@ -99,8 +97,8 @@ FloatArray z1 = r * ((car_coeffs_.a0_coeffs_ * car_state_.z1_memory_) - (car_coeffs_.c0_coeffs_ * car_state_.z2_memory_)); car_state_.z2_memory_ = r * - ((car_coeffs_.c0_coeffs_ * car_state_.z1_memory_) + - (car_coeffs_.a0_coeffs_ * car_state_.z2_memory_)); + ((car_coeffs_.c0_coeffs_ * car_state_.z1_memory_) + + (car_coeffs_.a0_coeffs_ * car_state_.z2_memory_)); car_state_.zy_memory_ = car_coeffs_.h_coeffs_ * car_state_.z2_memory_; // This section ripples the input-output path, to avoid delay... // It's the only part that doesn't get computed "in parallel": @@ -112,7 +110,6 @@ car_state_.zy_memory_(ch) = in_out; } car_state_.z1_memory_ = z1; - *car_out = car_state_.zy_memory_; } // We start with a quadratic nonlinear function, and limit it via a @@ -127,46 +124,44 @@ // This step is a one sample-time update of the inner-hair-cell (IHC) model, // including the detection nonlinearity and either one or two capacitor state // variables. -void Ear::IHCStep(const FloatArray& car_out, FloatArray* ihc_out) { +void Ear::IHCStep(const FloatArray& car_out) { FloatArray ac_diff = car_out - ihc_state_.ac_coupler_; ihc_state_.ac_coupler_ = ihc_state_.ac_coupler_ + - (ihc_coeffs_.ac_coeff_ * ac_diff); + (ihc_coeffs_.ac_coeff_ * ac_diff); if (ihc_coeffs_.just_half_wave_rectify_) { FloatArray output(n_ch_); for (int ch = 0; ch < n_ch_; ++ch) { FPType a = (ac_diff(ch) > 0.0) ? ac_diff(ch) : 0.0; output(ch) = (a < 2) ? a : 2; } - *ihc_out = output; + ihc_state_.ihc_out_ = output; } else { FloatArray conductance = CARFACDetect(ac_diff); if (ihc_coeffs_.one_capacitor_) { - *ihc_out = conductance * ihc_state_.cap1_voltage_; + ihc_state_.ihc_out_ = conductance * ihc_state_.cap1_voltage_; ihc_state_.cap1_voltage_ = ihc_state_.cap1_voltage_ - - (*ihc_out * ihc_coeffs_.out1_rate_) + - ((1 - ihc_state_.cap1_voltage_) * ihc_coeffs_.in1_rate_); + (ihc_state_.ihc_out_ * ihc_coeffs_.out1_rate_) + + ((1 - ihc_state_.cap1_voltage_) * ihc_coeffs_.in1_rate_); } else { - *ihc_out = conductance * ihc_state_.cap2_voltage_; + ihc_state_.ihc_out_ = conductance * ihc_state_.cap2_voltage_; ihc_state_.cap1_voltage_ = ihc_state_.cap1_voltage_ - - ((ihc_state_.cap1_voltage_ - ihc_state_.cap2_voltage_) - * ihc_coeffs_.out1_rate_) + - ((1 - ihc_state_.cap1_voltage_) * - ihc_coeffs_.in1_rate_); + ((ihc_state_.cap1_voltage_ - ihc_state_.cap2_voltage_) + * ihc_coeffs_.out1_rate_) + ((1 - ihc_state_.cap1_voltage_) * + ihc_coeffs_.in1_rate_); ihc_state_.cap2_voltage_ = ihc_state_.cap2_voltage_ - - (*ihc_out * ihc_coeffs_.out2_rate_) + - ((ihc_state_.cap1_voltage_ - ihc_state_.cap2_voltage_) - * ihc_coeffs_.in2_rate_); + (ihc_state_.ihc_out_ * ihc_coeffs_.out2_rate_) + + ((ihc_state_.cap1_voltage_ - ihc_state_.cap2_voltage_) + * ihc_coeffs_.in2_rate_); } // Here we smooth the output twice using a LPF. - *ihc_out *= ihc_coeffs_.output_gain_; + ihc_state_.ihc_out_ *= ihc_coeffs_.output_gain_; ihc_state_.lpf1_state_ += ihc_coeffs_.lpf_coeff_ * - (*ihc_out - ihc_state_.lpf1_state_); + (ihc_state_.ihc_out_ - ihc_state_.lpf1_state_); ihc_state_.lpf2_state_ += ihc_coeffs_.lpf_coeff_ * - (ihc_state_.lpf1_state_ - ihc_state_.lpf2_state_); - *ihc_out = ihc_state_.lpf2_state_ - ihc_coeffs_.rest_output_; + (ihc_state_.lpf1_state_ - ihc_state_.lpf2_state_); + ihc_state_.ihc_out_ = ihc_state_.lpf2_state_ - ihc_coeffs_.rest_output_; } - ihc_state_.ihc_out_ = *ihc_out; - ihc_state_.ihc_accum_ += *ihc_out; + ihc_state_.ihc_accum_ += ihc_state_.ihc_out_; } bool Ear::AGCStep(const FloatArray& ihc_out) { @@ -203,14 +198,12 @@ // Afterwards we add its output to this stage input, whether it updated or // not. agc_in += agc_coeffs.agc_stage_gain_ * - agc_state_[stage + 1].agc_memory_; + agc_state_[stage + 1].agc_memory_; } - FloatArray agc_stage_state = agc_state.agc_memory_; // This performs a first-order recursive smoothing filter update, in time. - agc_stage_state += agc_coeffs.agc_epsilon_ * - (agc_in - agc_stage_state); - agc_stage_state = AGCSpatialSmooth(stage, agc_stage_state); - agc_state.agc_memory_ = agc_stage_state; + agc_state.agc_memory_ += agc_coeffs.agc_epsilon_ * + (agc_in - agc_state.agc_memory_); + AGCSpatialSmooth(agc_coeffs_[stage], &agc_state.agc_memory_); updated = true; } else { updated = false; @@ -218,90 +211,81 @@ return updated; } -// TODO (alexbrandmeyer): figure out how to operate directly on stage_state. -// Using a pointer breaks the () indexing of the Eigen arrays, but there must -// be a way around this. -FloatArray Ear::AGCSpatialSmooth(const int stage, FloatArray stage_state) { - int n_iterations = agc_coeffs_[stage].agc_spatial_iterations_; +void Ear::AGCSpatialSmooth(const AGCCoeffs& agc_coeffs, + FloatArray* stage_state) { + int n_iterations = agc_coeffs.agc_spatial_iterations_; bool use_fir; use_fir = (n_iterations < 4) ? true : false; if (use_fir) { - FPType fir_coeffs_left = agc_coeffs_[stage].agc_spatial_fir_left_; - FPType fir_coeffs_mid = agc_coeffs_[stage].agc_spatial_fir_mid_; - FPType fir_coeffs_right = agc_coeffs_[stage].agc_spatial_fir_right_; + FPType fir_coeffs_left = agc_coeffs.agc_spatial_fir_left_; + FPType fir_coeffs_mid = agc_coeffs.agc_spatial_fir_mid_; + FPType fir_coeffs_right = agc_coeffs.agc_spatial_fir_right_; FloatArray ss_tap1(n_ch_); FloatArray ss_tap2(n_ch_); FloatArray ss_tap3(n_ch_); FloatArray ss_tap4(n_ch_); - int n_taps = agc_coeffs_[stage].agc_spatial_n_taps_; + int n_taps = agc_coeffs.agc_spatial_n_taps_; // This initializes the first two taps of stage state, which are used for // both possible cases. - ss_tap1(0) = stage_state(0); - ss_tap1.block(1, 0, n_ch_ - 1, 1) = stage_state.block(0, 0, n_ch_ - 1, 1); - ss_tap2(n_ch_ - 1) = stage_state(n_ch_ - 1); - ss_tap2.block(0, 0, n_ch_ - 1, 1) = stage_state.block(1, 0, n_ch_ - 1, 1); + ss_tap1(0) = (*stage_state)(0); + ss_tap1.block(1, 0, n_ch_ - 1, 1) = stage_state->block(0, 0, n_ch_ - 1, 1); + ss_tap2(n_ch_ - 1) = (*stage_state)(n_ch_ - 1); + ss_tap2.block(0, 0, n_ch_ - 1, 1) = stage_state->block(1, 0, n_ch_ - 1, 1); switch (n_taps) { case 3: - stage_state = (fir_coeffs_left * ss_tap1) + - (fir_coeffs_mid * stage_state) + - (fir_coeffs_right * ss_tap2); + *stage_state = (fir_coeffs_left * ss_tap1) + + (fir_coeffs_mid * *stage_state) + (fir_coeffs_right * ss_tap2); break; case 5: // Now we initialize last two taps of stage state, which are only used // for the 5-tap case. - ss_tap3(0) = stage_state(0); - ss_tap3(1) = stage_state(1); + ss_tap3(0) = (*stage_state)(0); + ss_tap3(1) = (*stage_state)(1); ss_tap3.block(2, 0, n_ch_ - 2, 1) = - stage_state.block(0, 0, n_ch_ - 2, 1); - ss_tap4(n_ch_ - 2) = stage_state(n_ch_ - 1); - ss_tap4(n_ch_ - 1) = stage_state(n_ch_ - 2); + stage_state->block(0, 0, n_ch_ - 2, 1); + ss_tap4(n_ch_ - 2) = (*stage_state)(n_ch_ - 1); + ss_tap4(n_ch_ - 1) = (*stage_state)(n_ch_ - 2); ss_tap4.block(0, 0, n_ch_ - 2, 1) = - stage_state.block(2, 0, n_ch_ - 2, 1); - stage_state = (fir_coeffs_left * (ss_tap3 + ss_tap1)) + - (fir_coeffs_mid * stage_state) + - (fir_coeffs_right * (ss_tap2 + ss_tap4)); + stage_state->block(2, 0, n_ch_ - 2, 1); + *stage_state = (fir_coeffs_left * (ss_tap3 + ss_tap1)) + + (fir_coeffs_mid * *stage_state) + + (fir_coeffs_right * (ss_tap2 + ss_tap4)); break; default: assert(true && "Bad n_taps in AGCSpatialSmooth; should be 3 or 5."); break; } } else { - stage_state = AGCSmoothDoubleExponential(stage_state, - agc_coeffs_[stage].agc_pole_z1_, - agc_coeffs_[stage].agc_pole_z2_); + AGCSmoothDoubleExponential(agc_coeffs.agc_pole_z1_, + agc_coeffs.agc_pole_z2_, stage_state); } - return stage_state; } -// TODO (alexbrandmeyer): figure out how to operate directly on stage_state. -// Same point as above for AGCSpatialSmooth. -FloatArray Ear::AGCSmoothDoubleExponential(FloatArray stage_state, - const FPType pole_z1, - const FPType pole_z2) { - int32_t n_pts = stage_state.size(); +void Ear::AGCSmoothDoubleExponential(const FPType pole_z1, const FPType pole_z2, + FloatArray* stage_state) { + int32_t n_pts = stage_state->size(); FPType input; FPType state = 0.0; // TODO (alexbrandmeyer): I'm assuming one dimensional input for now, but this // should be verified with Dick for the final version - for (int i = n_pts - 11; i < n_pts; ++i){ - input = stage_state(i); + for (int i = n_pts - 11; i < n_pts; ++i) { + input = (*stage_state)(i); state = state + (1 - pole_z1) * (input - state); } - for (int i = n_pts - 1; i > -1; --i){ - input = stage_state(i); + for (int i = n_pts - 1; i > -1; --i) { + input = (*stage_state)(i); state = state + (1 - pole_z2) * (input - state); } - for (int i = 0; i < n_pts; ++i){ - input = stage_state(i); + for (int i = 0; i < n_pts; ++i) { + input = (*stage_state)(i); state = state + (1 - pole_z1) * (input - state); - stage_state(i) = state; + (*stage_state)(i) = state; } - return stage_state; } FloatArray Ear::StageGValue(const FloatArray& undamping) { FloatArray r = car_coeffs_.r1_coeffs_ + car_coeffs_.zr_coeffs_ * undamping; return (1 - 2 * r * car_coeffs_.a0_coeffs_ + (r * r)) / - (1 - 2 * r * car_coeffs_.a0_coeffs_ + - car_coeffs_.h_coeffs_ * r * car_coeffs_.c0_coeffs_ + (r * r)); -} + (1 - 2 * r * car_coeffs_.a0_coeffs_ + car_coeffs_.h_coeffs_ * r * + car_coeffs_.c0_coeffs_ + (r * r)); +} \ No newline at end of file
--- a/trunk/carfac/ear.h Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/ear.h Tue May 28 15:54:54 2013 +0000 @@ -31,19 +31,22 @@ public: // This is the primary initialization function that is called for each // Ear object in the CARFAC 'Design' method. - void DesignEar(const int n_ch, const FPType fs, + void InitEar(const int n_ch, const FPType fs, const CARCoeffs& car_params, const IHCCoeffs& ihc_params, const std::vector<AGCCoeffs>& agc_params); // These three methods apply the different stages of the model in sequence // to individual audio samples. - void CARStep(const FPType input, FloatArray* car_out); - void IHCStep(const FloatArray& car_out, FloatArray* ihc_out); + void CARStep(const FPType input); + void IHCStep(const FloatArray& car_out); bool AGCStep(const FloatArray& ihc_out); // These accessor functions return portions of the CAR state for storage in // the CAROutput structures. const FloatArray& za_memory() { return car_state_.za_memory_; } const FloatArray& zb_memory() { return car_state_.zb_memory_; } + // The zy_memory_ of the CARState is equivalent to the CAR output. A second + // accessor function is included for documentation purposes. const FloatArray& zy_memory() { return car_state_.zy_memory_; } + const FloatArray& car_out() { return car_state_.zy_memory_; } const FloatArray& g_memory() { return car_state_.g_memory_; } // This returns the IHC output for storage. const FloatArray& ihc_out() { return ihc_state_.ihc_out_; } @@ -82,10 +85,9 @@ void OHCNonlinearFunction(const FloatArray& velocities, FloatArray* nonlinear_fun); bool AGCRecurse(const int stage, FloatArray agc_in); - FloatArray AGCSpatialSmooth(const int stage, FloatArray stage_state); - FloatArray AGCSmoothDoubleExponential(FloatArray stage_state, - const FPType pole_z1, - const FPType pole_z2); + void AGCSpatialSmooth(const AGCCoeffs& agc_coeffs , FloatArray* stage_state); + void AGCSmoothDoubleExponential(const FPType pole_z1, const FPType pole_z2, + FloatArray* stage_state); // These are the private data members that store the state and coefficient // information. CARCoeffs car_coeffs_;
--- a/trunk/carfac/ihc_params.cc Mon May 27 16:36:54 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,36 +0,0 @@ -// -// ihc_params.cc -// CARFAC Open Source C++ Library -// -// Created by Alex Brandmeyer on 5/10/13. -// -// This C++ file is part of an implementation of Lyon's cochlear model: -// "Cascade of Asymmetric Resonators with Fast-Acting Compression" -// to supplement Lyon's upcoming book "Human and Machine Hearing" -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -#include "ihc_params.h" - -// The default constructor for IHCParams initializes with the settings from -// Lyon's book 'Human and Machine Hearing' -IHCParams::IHCParams() { - just_half_wave_rectify_ = false; - one_capacitor_ = true; - tau_lpf_ = 0.000080; - tau1_out_ = 0.0005; - tau1_in_ = 0.010; - tau2_out_ = 0.0025; - tau2_in_ = 0.005; - ac_corner_hz_ = 20.0; -} \ No newline at end of file
--- a/trunk/carfac/ihc_params.h Mon May 27 16:36:54 2013 +0000 +++ b/trunk/carfac/ihc_params.h Tue May 28 15:54:54 2013 +0000 @@ -26,7 +26,16 @@ #include "carfac_common.h" struct IHCParams { - IHCParams(); + IHCParams() { + just_half_wave_rectify_ = false; + one_capacitor_ = true; + tau_lpf_ = 0.000080; + tau1_out_ = 0.0005; + tau1_in_ = 0.010; + tau2_out_ = 0.0025; + tau2_in_ = 0.005; + ac_corner_hz_ = 20.0; + }; bool just_half_wave_rectify_; bool one_capacitor_; FPType tau_lpf_;