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"
alexbrandmeyer@609: void CARFAC::Design(int n_ears, long fs, CARParams car_params,
alexbrandmeyer@610:                     IHCParams ihc_params, AGCParams agc_params) {
alexbrandmeyer@609:   n_ears_ = n_ears;
alexbrandmeyer@609:   fs_ = fs;
alexbrandmeyer@609:   ears_ = new Ear[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@610:     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@610:   FloatArray pole_freqs(n_ch_);
alexbrandmeyer@610:   pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
alexbrandmeyer@610:   for(int ch=0;ch < n_ch_; ch++) {
alexbrandmeyer@610:     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@610:   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@610:   for (int i = 0; i < n_ears_; i++) {
alexbrandmeyer@610:     ears_[i].InitEar(n_ch_, fs_, pole_freqs, car_params, ihc_params, agc_params);
alexbrandmeyer@610:   }
alexbrandmeyer@609: }
alexbrandmeyer@609: 
alexbrandmeyer@610: CARFACOutput CARFAC::Run(FloatArray2d sound_data, bool open_loop) {
alexbrandmeyer@610:   // We initialize one output object to store the final output.
alexbrandmeyer@609:   CARFACOutput *output = new CARFACOutput();
alexbrandmeyer@610:   // A second object is used to store the output for the individual segments.
alexbrandmeyer@609:   CARFACOutput *seg_output = new CARFACOutput();
alexbrandmeyer@609:   int n_audio_channels = int(sound_data.cols());
alexbrandmeyer@610:   long seg_len = 441;  // We use a fixed segment length for now.
alexbrandmeyer@609:   long n_timepoints = sound_data.rows();
alexbrandmeyer@610:   long n_segs = ceil((n_timepoints * 1.0) / seg_len);
alexbrandmeyer@609:   output->InitOutput(n_audio_channels, n_ch_, n_timepoints);
alexbrandmeyer@609:   seg_output->InitOutput(n_audio_channels, n_ch_, seg_len);
alexbrandmeyer@610:   // These values store the start and endpoints for each segment
alexbrandmeyer@610:   long start;
alexbrandmeyer@610:   long length = seg_len;
alexbrandmeyer@610:   // This section loops over the individual audio segments.
alexbrandmeyer@610:   for (long 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@610:     // Once we've determined the start point and segment length, we run the
alexbrandmeyer@610:     // CARFAC model on the current segment.
alexbrandmeyer@610:     RunSegment(sound_data.block(start, 0, length, n_audio_channels),
alexbrandmeyer@610:                seg_output, open_loop);
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:     output->MergeOutput(*seg_output, start, length);
alexbrandmeyer@609:   }
alexbrandmeyer@609:   return *output;
alexbrandmeyer@609: }
alexbrandmeyer@609: 
alexbrandmeyer@610: void CARFAC::RunSegment(FloatArray2d sound_data, CARFACOutput *seg_output,
alexbrandmeyer@610:                         bool open_loop) {
alexbrandmeyer@610:   // The number of timepoints is determined from the length of the audio
alexbrandmeyer@610:   // segment.
alexbrandmeyer@609:   long n_timepoints = sound_data.rows();
alexbrandmeyer@610:   // The number of ears is equal to the number of audio channels. This could
alexbrandmeyer@610:   // potentially be removed since we already know the n_ears_ during the design
alexbrandmeyer@610:   // stage. 
alexbrandmeyer@609:   int n_ears = int(sound_data.cols());
alexbrandmeyer@610:   // A nested loop structure is used to iterate through the individual samples
alexbrandmeyer@610:   // for each ear (audio channel).
alexbrandmeyer@610:   bool updated;  // This variable is used by the AGC stage.
alexbrandmeyer@610:   for (long i = 0; i < n_timepoints; i++) {
alexbrandmeyer@610:     for (int j = 0; j < n_ears; j++) {
alexbrandmeyer@610:       // This stores the audio sample currently being processed.
alexbrandmeyer@610:       FPType input = sound_data(i, j);
alexbrandmeyer@610:       // Now we apply the three stages of the model in sequence to the current
alexbrandmeyer@610:       // audio sample.
alexbrandmeyer@609:       FloatArray car_out = ears_[j].CARStep(input);
alexbrandmeyer@609:       FloatArray ihc_out = ears_[j].IHCStep(car_out);
alexbrandmeyer@610:       updated = ears_[j].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@610:       seg_output->ears_[j].nap_.block(0, i, n_ch_, 1) = ihc_out;
alexbrandmeyer@610:       // TODO alexbrandmeyer: Check with Dick to determine the C++ strategy for
alexbrandmeyer@610:       // storing optional output structures.
alexbrandmeyer@610:       seg_output->ears_[j].bm_.block(0, i, n_ch_, 1) = car_out;
alexbrandmeyer@610:       seg_output->ears_[j].ohc_.block(0, 1, n_ch_, 1) =
alexbrandmeyer@610:         ears_[j].ReturnZAMemory();
alexbrandmeyer@610:       seg_output->ears_[j].agc_.block(0, i, n_ch_, 1) =
alexbrandmeyer@610:         ears_[j].ReturnZBMemory();
alexbrandmeyer@610:     }
alexbrandmeyer@610:     if (updated && n_ears > 1) { CrossCouple(); }
alexbrandmeyer@610:     if (! open_loop) { CloseAGCLoop(); }
alexbrandmeyer@610:   }
alexbrandmeyer@610: }
alexbrandmeyer@610: 
alexbrandmeyer@610: void CARFAC::CrossCouple() {
alexbrandmeyer@610:   for (int stage = 0; stage < ears_[0].ReturnAGCNStages(); stage++) {
alexbrandmeyer@610:     if (ears_[0].ReturnAGCStateDecimPhase(stage) > 0) {
alexbrandmeyer@610:       break;
alexbrandmeyer@610:     } else {
alexbrandmeyer@610:       FPType mix_coeff = ears_[0].ReturnAGCMixCoeff(stage);
alexbrandmeyer@610:       if (mix_coeff > 0) {
alexbrandmeyer@610:         FloatArray stage_state;
alexbrandmeyer@610:         FloatArray this_stage_values = FloatArray::Zero(n_ch_);
alexbrandmeyer@610:         for (int ear = 0; ear < n_ears_; ear++) {
alexbrandmeyer@610:           stage_state = ears_[ear].ReturnAGCStateMemory(stage);
alexbrandmeyer@610:           this_stage_values += stage_state;
alexbrandmeyer@610:         }
alexbrandmeyer@610:         this_stage_values /= n_ears_;
alexbrandmeyer@610:         for (int ear = 0; ear < n_ears_; ear++) {
alexbrandmeyer@610:           stage_state = ears_[ear].ReturnAGCStateMemory(stage);
alexbrandmeyer@610:           ears_[ear].SetAGCStateMemory(stage,
alexbrandmeyer@610:                                          stage_state + mix_coeff *
alexbrandmeyer@610:                                          (this_stage_values - stage_state));
alexbrandmeyer@610:         }
alexbrandmeyer@610:       }
alexbrandmeyer@609:     }
alexbrandmeyer@609:   }
alexbrandmeyer@609: }
alexbrandmeyer@610: 
alexbrandmeyer@610: void CARFAC::CloseAGCLoop() {
alexbrandmeyer@610:   for (int ear = 0; ear < n_ears_; ear++) {
alexbrandmeyer@610:     FloatArray undamping = 1 - ears_[ear].ReturnAGCStateMemory(1);
alexbrandmeyer@610:     // This updates the target stage gain for the new damping.
alexbrandmeyer@610:     ears_[ear].SetCARStateDZBMemory(ears_[ear].ReturnZRCoeffs() * undamping -
alexbrandmeyer@610:                                     ears_[ear].ReturnZBMemory() /
alexbrandmeyer@610:                                     ears_[ear].ReturnAGCDecimation(1));
alexbrandmeyer@610:     ears_[ear].SetCARStateDGMemory((ears_[ear].StageGValue(undamping) -
alexbrandmeyer@610:                                    ears_[ear].ReturnGMemory()) /
alexbrandmeyer@610:                                    ears_[ear].ReturnAGCDecimation(1));
alexbrandmeyer@610:   }
alexbrandmeyer@610: }