alexbrandmeyer@626: //
alexbrandmeyer@626: //  carfac_test.cc
alexbrandmeyer@626: //  CARFAC Open Source C++ Library
alexbrandmeyer@626: //
alexbrandmeyer@626: //  Created by Alex Brandmeyer on 5/22/13.
alexbrandmeyer@626: //
alexbrandmeyer@626: // This C++ file is part of an implementation of Lyon's cochlear model:
alexbrandmeyer@626: // "Cascade of Asymmetric Resonators with Fast-Acting Compression"
alexbrandmeyer@626: // to supplement Lyon's upcoming book "Human and Machine Hearing"
alexbrandmeyer@626: //
alexbrandmeyer@626: // Licensed under the Apache License, Version 2.0 (the "License");
alexbrandmeyer@626: // you may not use this file except in compliance with the License.
alexbrandmeyer@626: // You may obtain a copy of the License at
alexbrandmeyer@626: //
alexbrandmeyer@626: //     http://www.apache.org/licenses/LICENSE-2.0
alexbrandmeyer@626: //
alexbrandmeyer@626: // Unless required by applicable law or agreed to in writing, software
alexbrandmeyer@626: // distributed under the License is distributed on an "AS IS" BASIS,
alexbrandmeyer@626: // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
alexbrandmeyer@626: // See the License for the specific language governing permissions and
alexbrandmeyer@626: // limitations under the License.
alexbrandmeyer@626: 
alexbrandmeyer@626: #include "carfac_test.h"
alexbrandmeyer@626: // Three helper functions are defined here for loading the test data generated
alexbrandmeyer@626: // by the Matlab version of CARFAC.
alexbrandmeyer@626: // This loads one-dimensional FloatArrays from single-column text files.
alexbrandmeyer@626: FloatArray LoadTestData(const std::string filename, const int number_points) {
alexbrandmeyer@626:   std::string fullfile = TEST_SRC_DIR + filename;
alexbrandmeyer@626:   std::ifstream file(fullfile.c_str());
alexbrandmeyer@626:   FPType myarray[number_points];
alexbrandmeyer@626:   FloatArray output(number_points);
alexbrandmeyer@626:   if (file.is_open()) {
alexbrandmeyer@626:     for (int i = 0; i < number_points; ++i) {
alexbrandmeyer@626:       file >> myarray[i];
alexbrandmeyer@626:       output(i) = myarray[i];
alexbrandmeyer@626:     }
alexbrandmeyer@626:   }
alexbrandmeyer@626:   return output;
alexbrandmeyer@626: }
alexbrandmeyer@626: 
alexbrandmeyer@626: // This loads two-dimensional FloatArrays from multi-column text files.
alexbrandmeyer@626: std::vector<FloatArray> Load2dTestData(const std::string filename, const int rows,
alexbrandmeyer@626:                             const int columns) {
alexbrandmeyer@626:   std::string fullfile = TEST_SRC_DIR + filename;
alexbrandmeyer@626:   std::ifstream file(fullfile.c_str());
alexbrandmeyer@626:   FPType myarray[rows][columns];
alexbrandmeyer@626:   std::vector<FloatArray> output;
alexbrandmeyer@626:   output.resize(rows);
alexbrandmeyer@626:   for (auto& timepoint : output) {
alexbrandmeyer@626:     timepoint.resize(columns);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   if (file.is_open()) {
alexbrandmeyer@626:     for (int i = 0; i < rows; ++i) {
alexbrandmeyer@626:       for (int j = 0; j < columns; ++j) {
alexbrandmeyer@626:         file >> myarray[i][j];
alexbrandmeyer@626:         output[i](j) = myarray[i][j];
alexbrandmeyer@626:       }
alexbrandmeyer@626:     }
alexbrandmeyer@626:   }
alexbrandmeyer@626:   return output;
alexbrandmeyer@626: }
alexbrandmeyer@626: 
alexbrandmeyer@626: // This loads two dimensional vectors of audio data using data generated in
alexbrandmeyer@626: // Matlab using the wavread() function.
alexbrandmeyer@626: std::vector<std::vector<float>> Load2dAudioVector(std::string filename,
alexbrandmeyer@626:                                                   int timepoints,
alexbrandmeyer@626:                                                   int channels) {
alexbrandmeyer@626:   std::string fullfile = TEST_SRC_DIR + filename;
alexbrandmeyer@626:   std::ifstream file(fullfile.c_str());
alexbrandmeyer@626:   std::vector<std::vector<float>> output;
alexbrandmeyer@626:   output.resize(channels);
alexbrandmeyer@626:   for (auto& channel : output) {
alexbrandmeyer@626:     channel.resize(timepoints);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   if (file.is_open()) {
alexbrandmeyer@626:     for (int i = 0; i < timepoints; ++i) {
alexbrandmeyer@626:       for (int j = 0; j < channels; ++j) {
alexbrandmeyer@626:         file >> output[j][i];
alexbrandmeyer@626:       }
alexbrandmeyer@626:     }
alexbrandmeyer@626:   }
alexbrandmeyer@626:   return output;
alexbrandmeyer@626: }
alexbrandmeyer@626: 
alexbrandmeyer@626: // The first test verifies that the resulting CAR coefficients are the same as
alexbrandmeyer@626: // in Matlab when using the default CAR parameter set.
alexbrandmeyer@626: TEST(CARFACTest, CARCoeffs_Test){
alexbrandmeyer@626:   // These initialze the CAR Params and Coeffs objects needed for this test.
alexbrandmeyer@626:   CARParams car_params;
alexbrandmeyer@626:   CARCoeffs car_coeffs;
alexbrandmeyer@626:   FPType fs = 22050.0;
alexbrandmeyer@626:   // We calculate the pole frequencies and number of channels in the same way
alexbrandmeyer@626:   // as in the CARFAC 'Design' method.
alexbrandmeyer@626:   int n_ch = 0;
alexbrandmeyer@626:   FPType pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
alexbrandmeyer@626:   while (pole_hz > car_params.min_pole_hz_) {
alexbrandmeyer@626:     n_ch++;
alexbrandmeyer@626:     pole_hz = pole_hz - car_params.erb_per_step_ *
alexbrandmeyer@626:     ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   FloatArray pole_freqs(n_ch);
alexbrandmeyer@626:   pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
alexbrandmeyer@626:   for (int ch = 0; ch < n_ch; ++ch) {
alexbrandmeyer@626:     pole_freqs(ch) = pole_hz;
alexbrandmeyer@626:     pole_hz = pole_hz - car_params.erb_per_step_ *
alexbrandmeyer@626:     ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   // This initializes the CAR coeffecients object and runs the design method.
alexbrandmeyer@626:   car_coeffs.Design(car_params, 22050, pole_freqs);
alexbrandmeyer@626:   // Now we go through each set of coefficients to verify that the values are
alexbrandmeyer@626:   // the same as in MATLAB.
alexbrandmeyer@626:   std::string filename;
alexbrandmeyer@626:   FloatArray output;
alexbrandmeyer@626:   
alexbrandmeyer@626:   ASSERT_EQ(car_coeffs.v_offset_, 0.04);
alexbrandmeyer@626:   ASSERT_EQ(car_coeffs.velocity_scale_, 0.1);
alexbrandmeyer@626:   
alexbrandmeyer@626:   filename = "r1_coeffs.txt";
alexbrandmeyer@626:   output = LoadTestData(filename, n_ch);
alexbrandmeyer@626:   for (int i = 0; i < n_ch; ++i) {
alexbrandmeyer@626:     ASSERT_NEAR(output(i), car_coeffs.r1_coeffs_(i), PRECISION_LEVEL);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   
alexbrandmeyer@626:   filename = "a0_coeffs.txt";
alexbrandmeyer@626:   output = LoadTestData(filename, n_ch);
alexbrandmeyer@626:   for (int i = 0; i < n_ch; ++i) {
alexbrandmeyer@626:     ASSERT_NEAR(output(i), car_coeffs.a0_coeffs_(i), PRECISION_LEVEL);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   
alexbrandmeyer@626:   filename = "c0_coeffs.txt";
alexbrandmeyer@626:   output = LoadTestData(filename, n_ch);
alexbrandmeyer@626:   for (int i = 0; i < n_ch; ++i) {
alexbrandmeyer@626:     ASSERT_NEAR(output(i), car_coeffs.c0_coeffs_(i), PRECISION_LEVEL);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   
alexbrandmeyer@626:   filename = "zr_coeffs.txt";
alexbrandmeyer@626:   output = LoadTestData(filename, n_ch);
alexbrandmeyer@626:   for (int i = 0; i < n_ch; ++i) {
alexbrandmeyer@626:     ASSERT_NEAR(output(i), car_coeffs.zr_coeffs_(i), PRECISION_LEVEL);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   
alexbrandmeyer@626:   filename = "h_coeffs.txt";
alexbrandmeyer@626:   output = LoadTestData(filename, n_ch);
alexbrandmeyer@626:   for (int i = 0; i < n_ch; ++i) {
alexbrandmeyer@626:     ASSERT_NEAR(output(i), car_coeffs.h_coeffs_(i), PRECISION_LEVEL);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   
alexbrandmeyer@626:   filename = "g0_coeffs.txt";
alexbrandmeyer@626:   output = LoadTestData(filename, n_ch);
alexbrandmeyer@626:   for (int i = 0; i < n_ch; ++i) {
alexbrandmeyer@626:     ASSERT_NEAR(output(i), car_coeffs.g0_coeffs_(i), PRECISION_LEVEL);
alexbrandmeyer@626:   }
alexbrandmeyer@626: }
alexbrandmeyer@626: 
alexbrandmeyer@626: // The second test verifies that the IHC coefficient calculations result in the
alexbrandmeyer@626: // same set of values as in the Matlab version of the CARFAC.
alexbrandmeyer@626: TEST(CARFACTest, IHCCoeffs_Test){
alexbrandmeyer@626:   IHCParams ihc_params;
alexbrandmeyer@626:   IHCCoeffs ihc_coeffs;
alexbrandmeyer@626:   FPType fs = 22050.0;
alexbrandmeyer@626:   ihc_coeffs.Design(ihc_params, fs);
alexbrandmeyer@626:   
alexbrandmeyer@626:   std::string filename = "ihc_coeffs.txt";
alexbrandmeyer@626:   FloatArray output = LoadTestData(filename, 9);
alexbrandmeyer@626:   
alexbrandmeyer@626:   // The sequence of the individual coefficients is determined using the
alexbrandmeyer@626:   // CARFAC_GenerateTestData() function in the Matlab version, with all of the
alexbrandmeyer@626:   // parameters placed in a single output file for convenience.
alexbrandmeyer@626:   bool just_hwr = output(0);
alexbrandmeyer@626:   FPType lpf_coeff = output(1);
alexbrandmeyer@626:   FPType out_rate = output(2);
alexbrandmeyer@626:   FPType in_rate = output(3);
alexbrandmeyer@626:   bool one_cap = output(4);
alexbrandmeyer@626:   FPType output_gain = output(5);
alexbrandmeyer@626:   FPType rest_output = output(6);
alexbrandmeyer@626:   FPType rest_cap = output(7);
alexbrandmeyer@626:   FPType ac_coeff = output(8);
alexbrandmeyer@626:   
alexbrandmeyer@626:   // Once we have the Matlab values initialized, we can compare them to the
alexbrandmeyer@626:   // output of the IHCCoeffs 'Design' method.
alexbrandmeyer@626:   ASSERT_EQ(just_hwr, ihc_coeffs.just_hwr_);
alexbrandmeyer@626:   ASSERT_NEAR(lpf_coeff, ihc_coeffs.lpf_coeff_, PRECISION_LEVEL);
alexbrandmeyer@626:   ASSERT_NEAR(out_rate, ihc_coeffs.out1_rate_, PRECISION_LEVEL);
alexbrandmeyer@626:   ASSERT_NEAR(in_rate, ihc_coeffs.in1_rate_, PRECISION_LEVEL);
alexbrandmeyer@626:   ASSERT_EQ(one_cap, ihc_coeffs.one_cap_);
alexbrandmeyer@626:   ASSERT_NEAR(output_gain, ihc_coeffs.output_gain_, PRECISION_LEVEL);
alexbrandmeyer@626:   ASSERT_NEAR(rest_output, ihc_coeffs.rest_output_, PRECISION_LEVEL);
alexbrandmeyer@626:   ASSERT_NEAR(rest_cap, ihc_coeffs.rest_cap1_, PRECISION_LEVEL);
alexbrandmeyer@626:   ASSERT_NEAR(ac_coeff, ihc_coeffs.ac_coeff_, PRECISION_LEVEL);
alexbrandmeyer@626: }
alexbrandmeyer@626: 
alexbrandmeyer@626: 
alexbrandmeyer@626: TEST(CARFACTest, AGCCoeffs_Test) {
alexbrandmeyer@626:   AGCParams agc_params;
alexbrandmeyer@626:   std::vector<AGCCoeffs> agc_coeffs;
alexbrandmeyer@626:   std::vector<FloatArray> output;
alexbrandmeyer@626:   output.resize(agc_params.n_stages_);
alexbrandmeyer@626:   std::string filename = "agc_coeffs_1.txt";
alexbrandmeyer@626:   output[0] = LoadTestData(filename, 14);
alexbrandmeyer@626:   filename = "agc_coeffs_2.txt";
alexbrandmeyer@626:   output[1] = LoadTestData(filename, 14);
alexbrandmeyer@626:   filename = "agc_coeffs_3.txt";
alexbrandmeyer@626:   output[2] = LoadTestData(filename, 14);
alexbrandmeyer@626:   filename = "agc_coeffs_4.txt";
alexbrandmeyer@626:   output[3] = LoadTestData(filename, 14);
alexbrandmeyer@626:   agc_coeffs.resize(agc_params.n_stages_);
alexbrandmeyer@626:   // We initialize the AGC stages in the same was as in Ear::Init.
alexbrandmeyer@626:   FPType fs = 22050.0;
alexbrandmeyer@626:   FPType previous_stage_gain = 0.0;
alexbrandmeyer@626:   FPType decim = 1.0;
alexbrandmeyer@626:   for (int stage = 0; stage < agc_params.n_stages_; ++stage) {
alexbrandmeyer@626:     agc_coeffs[stage].Design(agc_params, stage, fs, previous_stage_gain, decim);
alexbrandmeyer@626:     previous_stage_gain = agc_coeffs[stage].agc_gain_;
alexbrandmeyer@626:     decim = agc_coeffs[stage].decim_;
alexbrandmeyer@626:   }
alexbrandmeyer@626:   // Now we run through the individual coefficients and verify that they're the
alexbrandmeyer@626:   // same as in Matlab.
alexbrandmeyer@626:   for (int stage = 0; stage < agc_params.n_stages_; ++stage) {
alexbrandmeyer@626:     int n_agc_stages = output[stage](1);
alexbrandmeyer@626:     FPType agc_stage_gain = output[stage](2);
alexbrandmeyer@626:     int decimation = output[stage](3);
alexbrandmeyer@626:     FPType agc_epsilon = output[stage](4);
alexbrandmeyer@626:     FPType agc_polez1 = output[stage](5);
alexbrandmeyer@626:     FPType agc_polez2 = output[stage](6);
alexbrandmeyer@626:     int agc_spatial_iterations = output[stage](7);
alexbrandmeyer@626:     FPType agc_spatial_fir_1 = output[stage](8);
alexbrandmeyer@626:     FPType agc_spatial_fir_2 = output[stage](9);
alexbrandmeyer@626:     FPType agc_spatial_fir_3 = output[stage](10);
alexbrandmeyer@626:     int agc_spatial_n_taps = output[stage](11);
alexbrandmeyer@626:     FPType agc_mix_coeffs = output[stage](12);
alexbrandmeyer@626:     FPType detect_scale = output[stage](13);
alexbrandmeyer@626:     
alexbrandmeyer@626:     ASSERT_EQ(n_agc_stages, agc_coeffs[stage].n_agc_stages_);
alexbrandmeyer@626:     ASSERT_NEAR(agc_stage_gain, agc_coeffs[stage].agc_stage_gain_,
alexbrandmeyer@626:                 PRECISION_LEVEL);
alexbrandmeyer@626:     ASSERT_EQ(decimation, agc_coeffs[stage].decimation_);
alexbrandmeyer@626:     ASSERT_NEAR(agc_epsilon, agc_coeffs[stage].agc_epsilon_, PRECISION_LEVEL);
alexbrandmeyer@626:     ASSERT_NEAR(agc_polez1, agc_coeffs[stage].agc_pole_z1_, PRECISION_LEVEL);
alexbrandmeyer@626:     ASSERT_NEAR(agc_polez2, agc_coeffs[stage].agc_pole_z2_, PRECISION_LEVEL);
alexbrandmeyer@626:     ASSERT_EQ(agc_spatial_iterations,
alexbrandmeyer@626:               agc_coeffs[stage].agc_spatial_iterations_);
alexbrandmeyer@626:     ASSERT_NEAR(agc_spatial_fir_1, agc_coeffs[stage].agc_spatial_fir_[0],
alexbrandmeyer@626:                 PRECISION_LEVEL);
alexbrandmeyer@626:     ASSERT_NEAR(agc_spatial_fir_2, agc_coeffs[stage].agc_spatial_fir_[1],
alexbrandmeyer@626:                 PRECISION_LEVEL);
alexbrandmeyer@626:     ASSERT_EQ(agc_spatial_n_taps,
alexbrandmeyer@626:               agc_coeffs[stage].agc_spatial_n_taps_);
alexbrandmeyer@626:     ASSERT_NEAR(agc_spatial_fir_3, agc_coeffs[stage].agc_spatial_fir_[2],
alexbrandmeyer@626:                 PRECISION_LEVEL);
alexbrandmeyer@626:     ASSERT_NEAR(agc_mix_coeffs, agc_coeffs[stage].agc_mix_coeffs_,
alexbrandmeyer@626:                 PRECISION_LEVEL);
alexbrandmeyer@626:     
alexbrandmeyer@626:     // The last stage will have the correct detect_scale_ value on the basis of
alexbrandmeyer@626:     // the total gain accumlated over the stages.
alexbrandmeyer@626:     if (stage == agc_params.n_stages_ - 1) {
alexbrandmeyer@626:       ASSERT_NEAR(detect_scale, agc_coeffs[stage].detect_scale_,
alexbrandmeyer@626:                   PRECISION_LEVEL);
alexbrandmeyer@626:     }
alexbrandmeyer@626:   }
alexbrandmeyer@626: }
alexbrandmeyer@626: 
alexbrandmeyer@626: // This test verifies the output of the C++ code relative to that of the Matlab
alexbrandmeyer@626: // version using a single segment (441 samples) of audio from the "plan.wav"
alexbrandmeyer@626: // file. The single-channel audio data and different output matrices from Matlab
alexbrandmeyer@626: // are stored in text files and then read into 2d Eigen arrays (for now, this
alexbrandmeyer@626: // should be changed to a vector of FloatArrays... TODO (alexbrandmeyer)). For
alexbrandmeyer@626: // reference, see the CARFAC_GenerateTestData() function in the Matlab branch
alexbrandmeyer@626: // of the repository.
alexbrandmeyer@626: //
alexbrandmeyer@626: // A single Ear object is used along with the code from CARFAC.RunSegment() to
alexbrandmeyer@626: // evaluate the output of the CAR and IHC steps on a sample by sample basis
alexbrandmeyer@626: // relative to the output read in from Matlab. The test passes with 11 degrees
alexbrandmeyer@626: // of precision, with the Matlab data stored using 12 decimals.
alexbrandmeyer@626: //
alexbrandmeyer@626: // TODO (alexbrandmeyer): A subseqent version of this test will operate directly
alexbrandmeyer@626: // on the CARFACOutput structure and will evaluate binaural data.
alexbrandmeyer@626: TEST(CARFACTest, Monaural_Output_Test) {
alexbrandmeyer@626:   std::string filename = "monaural_test_nap.txt";
alexbrandmeyer@626:   std::vector<FloatArray> nap = Load2dTestData(filename, 441, 71);
alexbrandmeyer@626:   filename = "monaural_test_bm.txt";
alexbrandmeyer@626:   std::vector<FloatArray> bm = Load2dTestData(filename, 441, 71);
alexbrandmeyer@626:   filename = "monaural_test_ohc.txt";
alexbrandmeyer@626:   std::vector<FloatArray> ohc = Load2dTestData(filename, 441, 71);
alexbrandmeyer@626:   filename = "monaural_test_agc.txt";
alexbrandmeyer@626:   std::vector<FloatArray> agc = Load2dTestData(filename, 441, 71);
alexbrandmeyer@626:   filename = "file_signal_monaural_test.txt";
alexbrandmeyer@626:   std::vector<std::vector<float>> sound_data = Load2dAudioVector(filename, 441,
alexbrandmeyer@626:                                                                  1);
alexbrandmeyer@626:   // The number of timepoints is determined from the length of the audio
alexbrandmeyer@626:   // segment.
alexbrandmeyer@626:   int32_t n_timepoints = sound_data[0].size();
alexbrandmeyer@626:   
alexbrandmeyer@626:   CARParams car_params;
alexbrandmeyer@626:   IHCParams ihc_params;
alexbrandmeyer@626:   AGCParams agc_params;
alexbrandmeyer@626:   FPType fs = 22050.0;
alexbrandmeyer@626:   int n_ch = 0;
alexbrandmeyer@626:   FPType pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
alexbrandmeyer@626:   while (pole_hz > car_params.min_pole_hz_) {
alexbrandmeyer@626:     n_ch++;
alexbrandmeyer@626:     pole_hz = pole_hz - car_params.erb_per_step_ *
alexbrandmeyer@626:     ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   FloatArray pole_freqs(n_ch);
alexbrandmeyer@626:   pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
alexbrandmeyer@626:   for (int ch = 0; ch < n_ch; ++ch) {
alexbrandmeyer@626:     pole_freqs(ch) = pole_hz;
alexbrandmeyer@626:     pole_hz = pole_hz - car_params.erb_per_step_ *
alexbrandmeyer@626:     ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_);
alexbrandmeyer@626:   }
alexbrandmeyer@626:   
alexbrandmeyer@626:   // This initializes the CARFAC object and runs the design method.
alexbrandmeyer@626:   Ear ear;
alexbrandmeyer@626:   ear.InitEar(n_ch, fs, pole_freqs, car_params, ihc_params, agc_params);
alexbrandmeyer@626:   
alexbrandmeyer@626:   CARFACOutput seg_output;
alexbrandmeyer@626:   seg_output.InitOutput(1, n_ch, n_timepoints);
alexbrandmeyer@626:   
alexbrandmeyer@626:   // A nested loop structure is used to iterate through the individual samples
alexbrandmeyer@626:   // for each ear (audio channel).
alexbrandmeyer@626:   FloatArray car_out(n_ch);
alexbrandmeyer@626:   FloatArray ihc_out(n_ch);
alexbrandmeyer@626:   FloatArray matlab_car_out(n_ch);
alexbrandmeyer@626:   FloatArray matlab_ihc_out(n_ch);
alexbrandmeyer@626:   bool updated;  // This variable is used by the AGC stage.
alexbrandmeyer@626:   for (int32_t i = 0; i < n_timepoints; ++i) {
alexbrandmeyer@626:     int j = 0;
alexbrandmeyer@626:     // First we create a reference to the current Ear object.
alexbrandmeyer@626:     // This stores the audio sample currently being processed.
alexbrandmeyer@626:     FPType input = sound_data[j][i];
alexbrandmeyer@626:     // Now we apply the three stages of the model in sequence to the current
alexbrandmeyer@626:     // audio sample.
alexbrandmeyer@626:     ear.CARStep(input, &car_out);
alexbrandmeyer@626:     matlab_car_out = bm[i];
alexbrandmeyer@626:     // This step verifies that the ouput of the CAR step is the same at each
alexbrandmeyer@626:     // timepoint and channel as that of the Matlab version.
alexbrandmeyer@626:     for (int channel = 0; channel < n_ch; ++channel) {
alexbrandmeyer@626:       FPType a = matlab_car_out(channel);
alexbrandmeyer@626:       FPType b = car_out(channel);
alexbrandmeyer@626:       ASSERT_NEAR(a, b, PRECISION_LEVEL);
alexbrandmeyer@626:     }
alexbrandmeyer@626:     ear.IHCStep(car_out, &ihc_out);
alexbrandmeyer@626:     matlab_ihc_out = nap[i];
alexbrandmeyer@626:     // This step verifies that the ouput of the IHC step is the same at each
alexbrandmeyer@626:     // timepoint and channel as that of the Matlab version.
alexbrandmeyer@626:     for (int channel = 0; channel < n_ch; ++channel) {
alexbrandmeyer@626:       FPType a = matlab_ihc_out(channel);
alexbrandmeyer@626:       FPType b = ihc_out(channel);
alexbrandmeyer@626:       ASSERT_NEAR(a, b, PRECISION_LEVEL);
alexbrandmeyer@626:     }
alexbrandmeyer@626:     
alexbrandmeyer@626:     updated = ear.AGCStep(ihc_out);
alexbrandmeyer@626:     // These lines assign the output of the model for the current sample
alexbrandmeyer@626:     // to the appropriate data members of the current ear in the output
alexbrandmeyer@626:     // object.
alexbrandmeyer@626:     seg_output.StoreNAPOutput(i, j, ihc_out);
alexbrandmeyer@626:     seg_output.StoreBMOutput(i, j, car_out);
alexbrandmeyer@626:     seg_output.StoreOHCOutput(i, j, ear.za_memory());
alexbrandmeyer@626:     seg_output.StoreAGCOutput(i, j, ear.zb_memory());
alexbrandmeyer@626:     if (updated) {
alexbrandmeyer@626:       FloatArray undamping = 1 - ear.agc_memory(0);
alexbrandmeyer@626:       // This updates the target stage gain for the new damping.
alexbrandmeyer@626:       ear.set_dzb_memory((ear.zr_coeffs() * undamping - ear.zb_memory()) /
alexbrandmeyer@626:                          ear.agc_decimation(0));
alexbrandmeyer@626:       ear.set_dg_memory((ear.StageGValue(undamping) - ear.g_memory()) /
alexbrandmeyer@626:                         ear.agc_decimation(0));
alexbrandmeyer@626:     }
alexbrandmeyer@626:   }
alexbrandmeyer@626: }