alexbrandmeyer@609: //
alexbrandmeyer@609: //  carfac.cc
alexbrandmeyer@609: //  CARFAC Open Source C++ Library
alexbrandmeyer@609: //
alexbrandmeyer@609: //  Created by Alex Brandmeyer on 5/10/13.
alexbrandmeyer@609: //
alexbrandmeyer@609: // This C++ file is part of an implementation of Lyon's cochlear model:
alexbrandmeyer@609: // "Cascade of Asymmetric Resonators with Fast-Acting Compression"
alexbrandmeyer@609: // to supplement Lyon's upcoming book "Human and Machine Hearing"
alexbrandmeyer@609: //
alexbrandmeyer@609: // Licensed under the Apache License, Version 2.0 (the "License");
alexbrandmeyer@609: // you may not use this file except in compliance with the License.
alexbrandmeyer@609: // You may obtain a copy of the License at
alexbrandmeyer@609: //
alexbrandmeyer@609: //     http://www.apache.org/licenses/LICENSE-2.0
alexbrandmeyer@609: //
alexbrandmeyer@609: // Unless required by applicable law or agreed to in writing, software
alexbrandmeyer@609: // distributed under the License is distributed on an "AS IS" BASIS,
alexbrandmeyer@609: // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
alexbrandmeyer@609: // See the License for the specific language governing permissions and
alexbrandmeyer@609: // limitations under the License.
alexbrandmeyer@609: 
alexbrandmeyer@609: #include "carfac.h"
ronw@625: 
alexbrandmeyer@626: void CARFAC::Design(const int n_ears, const FPType fs,
alexbrandmeyer@626:                     const CARParams& car_params, const IHCParams& ihc_params,
alexbrandmeyer@626:                     const AGCParams& agc_params) {
alexbrandmeyer@609:   n_ears_ = n_ears;
alexbrandmeyer@609:   fs_ = fs;
alexbrandmeyer@611:   ears_.resize(n_ears_);
alexbrandmeyer@610:   
alexbrandmeyer@610:   n_ch_ = 0;
alexbrandmeyer@610:   FPType pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
alexbrandmeyer@610:   while (pole_hz > car_params.min_pole_hz_) {
alexbrandmeyer@626:     ++n_ch_;
alexbrandmeyer@610:     pole_hz = pole_hz - car_params.erb_per_step_ *
alexbrandmeyer@610:     ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_);
alexbrandmeyer@609:   }
alexbrandmeyer@626:   pole_freqs_.resize(n_ch_);
alexbrandmeyer@610:   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@610:     pole_hz = pole_hz - car_params.erb_per_step_ *
alexbrandmeyer@610:     ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_);
alexbrandmeyer@610:   }
alexbrandmeyer@626:   max_channels_per_octave_ = log(2) / log(pole_freqs_(0) / pole_freqs_(1));
alexbrandmeyer@610:   // Once we have the basic information about the pole frequencies and the
alexbrandmeyer@610:   // number of channels, we initialize the ear(s).
alexbrandmeyer@626:   for (auto& ear : ears_) {
alexbrandmeyer@626:     ear.InitEar(n_ch_, fs_, pole_freqs_, car_params, ihc_params,
alexbrandmeyer@626:                 agc_params);
alexbrandmeyer@610:   }
alexbrandmeyer@609: }
alexbrandmeyer@609: 
alexbrandmeyer@626: CARFACOutput CARFAC::Run(const std::vector<std::vector<float>>& sound_data) {
alexbrandmeyer@610:   // We initialize one output object to store the final output.
alexbrandmeyer@626:   CARFACOutput seg_output;
alexbrandmeyer@626:   int n_audio_channels = sound_data.size();
alexbrandmeyer@611:   int32_t seg_len = 441;  // We use a fixed segment length for now.
alexbrandmeyer@626:   int32_t n_timepoints = sound_data[0].size();
alexbrandmeyer@611:   int32_t n_segs = ceil((n_timepoints * 1.0) / seg_len);
alexbrandmeyer@626:   seg_output.InitOutput(n_audio_channels, n_ch_, n_timepoints);
alexbrandmeyer@610:   // These values store the start and endpoints for each segment
alexbrandmeyer@611:   int32_t start;
alexbrandmeyer@611:   int32_t length = seg_len;
alexbrandmeyer@610:   // This section loops over the individual audio segments.
alexbrandmeyer@626:   for (int32_t i = 0; i < n_segs; ++i) {
alexbrandmeyer@610:     // For each segment we calculate the start point and the segment length.
alexbrandmeyer@610:     start = i * seg_len;
alexbrandmeyer@610:     if (i == n_segs - 1) {
alexbrandmeyer@610:       // The last segment can be shorter than the rest.
alexbrandmeyer@610:       length = n_timepoints - start;
alexbrandmeyer@609:     }
alexbrandmeyer@626:     std::vector<std::vector<float>> segment_data;
alexbrandmeyer@626:     segment_data.resize(n_audio_channels);
alexbrandmeyer@626:     for (int channel = 0; channel < n_audio_channels; ++channel) {
alexbrandmeyer@626:       segment_data[channel].resize(length);
alexbrandmeyer@626:       for (int32_t timepoint = 0; timepoint < length; ++timepoint) {
alexbrandmeyer@626:         segment_data[channel][timepoint] =
alexbrandmeyer@626:         sound_data[channel][start + timepoint];
alexbrandmeyer@626:       }
alexbrandmeyer@626:     }
alexbrandmeyer@610:     // Once we've determined the start point and segment length, we run the
alexbrandmeyer@610:     // CARFAC model on the current segment.
alexbrandmeyer@626:     RunSegment(segment_data, start,
alexbrandmeyer@626:                length, &seg_output, true);
alexbrandmeyer@610:     // Afterwards we merge the output for the current segment into the larger
alexbrandmeyer@610:     // output structure for the entire audio file.
alexbrandmeyer@609:   }
alexbrandmeyer@626:   return seg_output;
alexbrandmeyer@609: }
alexbrandmeyer@609: 
alexbrandmeyer@626: void CARFAC::RunSegment(const std::vector<std::vector<float>>& sound_data,
alexbrandmeyer@626:                         const int32_t start, const int32_t length,
alexbrandmeyer@626:                         CARFACOutput* seg_output, const bool open_loop) {
alexbrandmeyer@626:   // The number of ears is equal to the number of audio channels. This could
alexbrandmeyer@626:   // potentially be removed since we already know the n_ears_ during the design
alexbrandmeyer@626:   // stage.
alexbrandmeyer@626:   int n_ears = sound_data.size();
alexbrandmeyer@610:   // The number of timepoints is determined from the length of the audio
alexbrandmeyer@610:   // segment.
alexbrandmeyer@626:   int32_t n_timepoints = sound_data[0].size();
alexbrandmeyer@610:   // A nested loop structure is used to iterate through the individual samples
alexbrandmeyer@610:   // for each ear (audio channel).
alexbrandmeyer@611:   FloatArray car_out(n_ch_);
alexbrandmeyer@611:   FloatArray 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:     for (int j = 0; j < n_ears; ++j) {
alexbrandmeyer@626:       // First we create a reference to the current Ear object.
alexbrandmeyer@626:       Ear& ear = ears_[j];
alexbrandmeyer@610:       // This stores the audio sample currently being processed.
alexbrandmeyer@626:       FPType input = sound_data[j][i];
alexbrandmeyer@610:       // Now we apply the three stages of the model in sequence to the current
alexbrandmeyer@610:       // audio sample.
alexbrandmeyer@626:       ear.CARStep(input, &car_out);
alexbrandmeyer@626:       ear.IHCStep(car_out, &ihc_out);
alexbrandmeyer@626:       updated = ear.AGCStep(ihc_out);
alexbrandmeyer@610:       // These lines assign the output of the model for the current sample
alexbrandmeyer@610:       // to the appropriate data members of the current ear in the output
alexbrandmeyer@610:       // object.
alexbrandmeyer@626:       seg_output->StoreNAPOutput(start + i, j, ihc_out);
alexbrandmeyer@626:       seg_output->StoreBMOutput(start + i, j, car_out);
alexbrandmeyer@626:       seg_output->StoreOHCOutput(start + i, j, ear.za_memory());
alexbrandmeyer@626:       seg_output->StoreAGCOutput(start + i, j, ear.zb_memory());
alexbrandmeyer@610:     }
alexbrandmeyer@626:     if (updated && n_ears > 1) {
alexbrandmeyer@626:       CrossCouple();
alexbrandmeyer@626:     }
alexbrandmeyer@626:     if (! open_loop) {
alexbrandmeyer@626:       CloseAGCLoop();
alexbrandmeyer@626:     }
alexbrandmeyer@610:   }
alexbrandmeyer@610: }
alexbrandmeyer@610: 
alexbrandmeyer@610: void CARFAC::CrossCouple() {
alexbrandmeyer@626:   for (int stage = 0; stage < ears_[0].agc_nstages(); ++stage) {
alexbrandmeyer@626:     if (ears_[0].agc_decim_phase(stage) > 0) {
alexbrandmeyer@610:       break;
alexbrandmeyer@610:     } else {
alexbrandmeyer@626:       FPType mix_coeff = ears_[0].agc_mix_coeff(stage);
alexbrandmeyer@610:       if (mix_coeff > 0) {
alexbrandmeyer@610:         FloatArray stage_state;
alexbrandmeyer@610:         FloatArray this_stage_values = FloatArray::Zero(n_ch_);
alexbrandmeyer@626:         for (auto& ear : ears_) {
alexbrandmeyer@626:           stage_state = ear.agc_memory(stage);
alexbrandmeyer@610:           this_stage_values += stage_state;
alexbrandmeyer@610:         }
alexbrandmeyer@610:         this_stage_values /= n_ears_;
alexbrandmeyer@626:         for (auto& ear : ears_) {
alexbrandmeyer@626:           stage_state = ear.agc_memory(stage);
alexbrandmeyer@626:           ear.set_agc_memory(stage, stage_state + mix_coeff *
alexbrandmeyer@626:                              (this_stage_values - stage_state));
alexbrandmeyer@610:         }
alexbrandmeyer@610:       }
alexbrandmeyer@609:     }
alexbrandmeyer@609:   }
alexbrandmeyer@609: }
alexbrandmeyer@610: 
alexbrandmeyer@610: void CARFAC::CloseAGCLoop() {
alexbrandmeyer@626:   for (auto& ear: ears_) {
alexbrandmeyer@626:     FloatArray undamping = 1 - ear.agc_memory(0);
alexbrandmeyer@610:     // 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@610:   }
alexbrandmeyer@626: }