diff carfac/carfac.cc @ 610:01986636257a

Second check-in of Alex Brandmeyer's C++ implementation of CARFAC. Addressed style issues and completed implementation of remaining functions. Still needs proper testing of the output stages against the MATLAB version, and runtime functions need improvements in efficiency.

--- a/carfac/carfac.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/carfac.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,61 +22,138 @@
 
 #include "carfac.h"
 void CARFAC::Design(int n_ears, long fs, CARParams car_params,
-                    IHCParams ihc_params, AGCParams agc_params){
+                    IHCParams ihc_params, AGCParams agc_params) {
   n_ears_ = n_ears;
   fs_ = fs;
   ears_ = new Ear[n_ears_];
-  for (int i = 0; i < n_ears_; i++){
-    ears_[i].InitEar(fs_, car_params,ihc_params,agc_params);
+  
+  n_ch_ = 0;
+  FPType pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
+  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_);
   }
-  n_ch_ = ears_[0].car_coeffs_.n_ch_;
+  FloatArray pole_freqs(n_ch_);
+  pole_hz = car_params.first_pole_theta_ * fs / (2 * PI);
+  for(int ch=0;ch < n_ch_; ch++) {
+    pole_freqs(ch) = pole_hz;
+    pole_hz = pole_hz - car_params.erb_per_step_ *
+    ERBHz(pole_hz, car_params.erb_break_freq_, car_params.erb_q_);
+  }
+  max_channels_per_octave_ = log(2) / log(pole_freqs(0) / pole_freqs(1));
+  // Once we have the basic information about the pole frequencies and the
+  // number of channels, we initialize the ear(s).
+  for (int i = 0; i < n_ears_; i++) {
+    ears_[i].InitEar(n_ch_, fs_, pole_freqs, car_params, ihc_params, agc_params);
+  }
 }
 
-CARFACOutput CARFAC::Run(FloatArray2d sound_data){
-  //to store the final output
+CARFACOutput CARFAC::Run(FloatArray2d sound_data, bool open_loop) {
+  // We initialize one output object to store the final output.
   CARFACOutput *output = new CARFACOutput();
-  //to store the output of the individual segments
+  // A second object is used to store the output for the individual segments.
   CARFACOutput *seg_output = new CARFACOutput();
-  
   int n_audio_channels = int(sound_data.cols());
-  long seg_len = 441; //Fixed segment length for now
+  long seg_len = 441;  // We use a fixed segment length for now.
   long n_timepoints = sound_data.rows();
-  double n_segs = ceil(double(n_timepoints)/double(seg_len));
+  long n_segs = ceil((n_timepoints * 1.0) / seg_len);
   output->InitOutput(n_audio_channels, n_ch_, n_timepoints);
   seg_output->InitOutput(n_audio_channels, n_ch_, seg_len);
-  //loop over individual audio segments
-  long start, length; //to store the start and endpoints for each segment
-  for (long i = 0; i < long(n_segs); i++){
-    //determine start and end points
-    start = (i * seg_len);
-    if (i < n_segs - 1){
-      length = seg_len;
-    } else {
-      length = n_timepoints - start; //the last segment can be shorter than the rest
+  // These values store the start and endpoints for each segment
+  long start;
+  long length = seg_len;
+  // This section loops over the individual audio segments.
+  for (long i = 0; i < n_segs; i++) {
+    // For each segment we calculate the start point and the segment length.
+    start = i * seg_len;
+    if (i == n_segs - 1) {
+      // The last segment can be shorter than the rest.
+      length = n_timepoints - start;
     }
-    RunSegment(sound_data.block(start,0,length,n_audio_channels),seg_output);
+    // Once we've determined the start point and segment length, we run the
+    // CARFAC model on the current segment.
+    RunSegment(sound_data.block(start, 0, length, n_audio_channels),
+               seg_output, open_loop);
+    // Afterwards we merge the output for the current segment into the larger
+    // output structure for the entire audio file.
     output->MergeOutput(*seg_output, start, length);
   }
-  
   return *output;
 }
 
-void CARFAC::RunSegment(FloatArray2d sound_data, CARFACOutput *seg_output){
+void CARFAC::RunSegment(FloatArray2d sound_data, CARFACOutput *seg_output,
+                        bool open_loop) {
+  // The number of timepoints is determined from the length of the audio
+  // segment.
   long n_timepoints = sound_data.rows();
+  // The number of ears is equal to the number of audio channels. This could
+  // potentially be removed since we already know the n_ears_ during the design
+  // stage. 
   int n_ears = int(sound_data.cols());
-  for (long i = 0; i < n_timepoints; i++){
-    for (int j = 0; j < n_ears; j++){
-      FPType input = sound_data(i,j);
+  // A nested loop structure is used to iterate through the individual samples
+  // for each ear (audio channel).
+  bool updated;  // This variable is used by the AGC stage.
+  for (long i = 0; i < n_timepoints; i++) {
+    for (int j = 0; j < n_ears; j++) {
+      // This stores the audio sample currently being processed.
+      FPType input = sound_data(i, j);
+      // Now we apply the three stages of the model in sequence to the current
+      // audio sample.
       FloatArray car_out = ears_[j].CARStep(input);
       FloatArray ihc_out = ears_[j].IHCStep(car_out);
-      bool updated = ears_[j].AGCStep(ihc_out);
-      seg_output->ears_[j].nap_.block(0,i,n_ch_,1) = ihc_out;
-      seg_output->ears_[j].bm_.block(0,i,n_ch_,1) = car_out;
-      seg_output->ears_[j].ohc_.block(0,1,n_ch_,1) =
-        ears_[j].car_state_.za_memory_;
-      seg_output->ears_[j].agc_.block(0,i,n_ch_,1) =
-        ears_[j].car_state_.zb_memory_;
+      updated = ears_[j].AGCStep(ihc_out);
+      // These lines assign the output of the model for the current sample
+      // to the appropriate data members of the current ear in the output
+      // object.
+      seg_output->ears_[j].nap_.block(0, i, n_ch_, 1) = ihc_out;
+      // TODO alexbrandmeyer: Check with Dick to determine the C++ strategy for
+      // storing optional output structures.
+      seg_output->ears_[j].bm_.block(0, i, n_ch_, 1) = car_out;
+      seg_output->ears_[j].ohc_.block(0, 1, n_ch_, 1) =
+        ears_[j].ReturnZAMemory();
+      seg_output->ears_[j].agc_.block(0, i, n_ch_, 1) =
+        ears_[j].ReturnZBMemory();
+    }
+    if (updated && n_ears > 1) { CrossCouple(); }
+    if (! open_loop) { CloseAGCLoop(); }
+  }
+}
+
+void CARFAC::CrossCouple() {
+  for (int stage = 0; stage < ears_[0].ReturnAGCNStages(); stage++) {
+    if (ears_[0].ReturnAGCStateDecimPhase(stage) > 0) {
+      break;
+    } else {
+      FPType mix_coeff = ears_[0].ReturnAGCMixCoeff(stage);
+      if (mix_coeff > 0) {
+        FloatArray stage_state;
+        FloatArray this_stage_values = FloatArray::Zero(n_ch_);
+        for (int ear = 0; ear < n_ears_; ear++) {
+          stage_state = ears_[ear].ReturnAGCStateMemory(stage);
+          this_stage_values += stage_state;
+        }
+        this_stage_values /= n_ears_;
+        for (int ear = 0; ear < n_ears_; ear++) {
+          stage_state = ears_[ear].ReturnAGCStateMemory(stage);
+          ears_[ear].SetAGCStateMemory(stage,
+                                         stage_state + mix_coeff *
+                                         (this_stage_values - stage_state));
+        }
+      }
     }
   }
-  
 }
+
+void CARFAC::CloseAGCLoop() {
+  for (int ear = 0; ear < n_ears_; ear++) {
+    FloatArray undamping = 1 - ears_[ear].ReturnAGCStateMemory(1);
+    // This updates the target stage gain for the new damping.
+    ears_[ear].SetCARStateDZBMemory(ears_[ear].ReturnZRCoeffs() * undamping -
+                                    ears_[ear].ReturnZBMemory() /
+                                    ears_[ear].ReturnAGCDecimation(1));
+    ears_[ear].SetCARStateDGMemory((ears_[ear].StageGValue(undamping) -
+                                   ears_[ear].ReturnGMemory()) /
+                                   ears_[ear].ReturnAGCDecimation(1));
+  }
+}