# HG changeset patch
# User alexbrandmeyer
# Date 1368725603 0
# Node ID 01986636257a8a7010ddf75d9bc6b392c6b491dc
# Parent  aefe2ca0674f4397914aec7880f0f1240fb63f7b
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.

diff -r aefe2ca0674f -r 01986636257a carfac/agc_coeffs.cc
--- a/carfac/agc_coeffs.cc	Mon May 13 22:51:15 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,167 +0,0 @@
-//
-//  agc_coeffs.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_coeffs.h"
-
-//This method has been created for debugging purposes and depends on <iostream>.
-//Could possibly be removed in the final version to reduce dependencies.
-void AGCCoeffs::OutputCoeffs(){
-  std::cout << "AGCCoeffs Values" << std::endl;
-  std::cout << "****************" << std::endl;
-  std::cout << "n_ch_ = " << n_ch_ << std::endl;
-  std::cout << "n_agc_stages_ = " << n_agc_stages_ << std::endl;
-  std::cout << "agc_stage_gain_ = " << agc_stage_gain_ << std::endl;
-  std::cout << "agc_epsilon_ = " << agc_epsilon_ << std::endl;
-  std::cout << "decimation_ = " << decimation_ << std::endl;
-  std::cout << "agc_pole_z1_ = " << agc_pole_z1_ << std::endl;
-  std::cout << "agc_pole_z2_ = " << agc_pole_z2_ << std::endl;
-  std::cout << "agc_spatial_iterations_ = " << agc_spatial_iterations_
-            << std::endl;
-  std::cout << "agc_spatial_fir_ = " << agc_spatial_fir_ << std::endl;
-  std::cout << "agc_spatial_n_taps_ = " << agc_spatial_n_taps_ << std::endl;
-  std::cout << "agc_mix_coeffs_ = " << agc_mix_coeffs_ << std::endl;
-  std::cout << "agc1_scales_ = " << agc1_scales_ << std::endl;
-  std::cout << "agc2_scales_ = " << agc2_scales_ << std::endl;
-  std::cout << "time_constants_ = " << time_constants_
-            << std::endl << std::endl;
-}
-
-void AGCCoeffs::DesignAGC(AGCParams agc_params, long fs, int n_ch){
-  n_ch_ = n_ch;
-  n_agc_stages_ = agc_params.n_stages_;
-  agc_stage_gain_ = agc_params.agc_stage_gain_;
-  agc1_scales_ = agc_params.agc1_scales_;
-  agc2_scales_ = agc_params.agc2_scales_;
-  time_constants_ = agc_params.time_constants_;
-  agc_epsilon_.resize(n_agc_stages_);
-  agc_pole_z1_.resize(n_agc_stages_);
-  agc_pole_z2_.resize(n_agc_stages_);
-  agc_spatial_iterations_.resize(n_agc_stages_);
-  agc_spatial_n_taps_.resize(n_agc_stages_);
-  agc_spatial_fir_.resize(3,n_agc_stages_);
-  agc_mix_coeffs_.resize(n_agc_stages_);
-  mix_coeff_ = agc_params.agc_mix_coeff_;
-  fir_.resize(3);
-  decim_ = 1;
-  decimation_ = agc_params.decimation_;
-  total_dc_gain_ = 0;
-  
-  for (int stage=0; stage < n_agc_stages_; stage++){
-    tau_ = time_constants_(stage);
-    decim_ = decim_ * decimation_(stage);
-    agc_epsilon_(stage) = 1 - exp((-1 * decim_) / (tau_ * fs));
-    n_times_ = tau_ * (fs / decim_);
-    delay_ = (agc2_scales_(stage) - agc1_scales_(stage)) / n_times_;
-    spread_sq_ = (pow(agc1_scales_(stage),2) + pow(agc2_scales_(stage),2)) /
-    n_times_;
-    u_ = 1 + (1 / spread_sq_);
-    p_ = u_ - sqrt(pow(u_,2) - 1);
-    dp_ = delay_ * (1 - (2 * p_) + pow(p_,2)) / 2;
-    pole_z1_ = p_ - dp_;
-    pole_z2_ = p_ + dp_;
-    agc_pole_z1_(stage) = pole_z1_;
-    agc_pole_z2_(stage) = pole_z2_;
-    n_taps_ = 0;
-    fir_ok_ = 0;
-    n_iterations_ = 1;
-    //initialize FIR coeffs settings
-    try {
-      while (! fir_ok_){
-        switch (n_taps_){
-          case 0:
-            n_taps_ = 3;
-            break;
-          case 3:
-            n_taps_ = 5;
-            break;
-          case 5:
-            n_iterations_ ++;
-            if (n_iterations_ > 16){
-              throw 10; //too many iterations
-            }
-            break;
-          default:
-            throw 20; //bad n_taps
-        }
-        //Design FIR Coeffs
-        FPType var = spread_sq_ / n_iterations_;
-        FPType mn = delay_ / n_iterations_;
-        switch (n_taps_){
-          case 3:
-            a = (var + pow(mn,2) - mn) / 2;
-            b = (var + pow(mn,2) + mn) / 2;
-            fir_ << a, 1 - a - b, b;
-            if (fir_(2) >= 0.2) {
-              fir_ok_ = true;
-            } else {
-              fir_ok_ = false;
-            }
-            break;
-          case 5:
-            a = (((var + pow(mn,2)) * 2/5) - (mn * 2/3)) / 2;
-            b = (((var + pow(mn,2)) * 2/5) + (mn * 2/3)) / 2;
-            fir_ << a/2, 1 - a - b, b/2;
-            if (fir_(2) >= 0.1) {
-              fir_ok_ = true;
-            } else {
-              fir_ok_ = false;
-            }
-            break;
-          default:
-            throw 30; //bad n_taps in FIR design
-        }
-      }
-    }
-    catch (int e) {
-      switch (e) {
-        case 10:
-          std::cout << "ERROR: Too many n_iterations in agc_coeffs.DesignAGC"
-          << std::endl;
-          break;
-        case 20:
-          std::cout << "ERROR: Bad n_taps in agc_coeffs.DesignAGC" << std::endl;
-          break;
-        case 30:
-          std::cout << "ERROR: Bad n_taps in agc_coeffs.DesignAGC/FIR"
-          << std::endl;
-          break;
-        default:
-          std::cout << "ERROR: unknown error in agc_coeffs.DesignAGC"
-          << std::endl;
-      }
-    }
-    //assign output of filter design
-    agc_spatial_iterations_(stage) = n_iterations_;
-    agc_spatial_n_taps_(stage) = n_taps_;
-    agc_spatial_fir_(0,stage) = fir_(0);
-    agc_spatial_fir_(1,stage) = fir_(1);
-    agc_spatial_fir_(2,stage) = fir_(2);
-    total_dc_gain_ = total_dc_gain_ + pow(agc_stage_gain_,(stage));
-    if (stage == 0) {
-      agc_mix_coeffs_(stage) = 0;
-    } else {
-      agc_mix_coeffs_(stage) = mix_coeff_ / (tau_ * (fs /decim_));
-    }
-  }
-  agc_gain_ = total_dc_gain_;
-  detect_scale_ = 1 / total_dc_gain_;
-}
diff -r aefe2ca0674f -r 01986636257a carfac/agc_coeffs.h
--- a/carfac/agc_coeffs.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/agc_coeffs.h	Thu May 16 17:33:23 2013 +0000
@@ -25,47 +25,19 @@
 
 #include "agc_params.h"
 
-class AGCCoeffs {
-public:
-  int n_ch_;
+struct AGCCoeffs {
   int n_agc_stages_;
   int agc_stage_gain_;
-  FloatArray agc_epsilon_; //FloatArray
-  FloatArray decimation_; //FloatArray
-  FloatArray agc_pole_z1_; //FloatArray
-  FloatArray agc_pole_z2_; //FloatArray
-  FloatArray agc_spatial_iterations_; //FloatArray
-  FloatArray2d agc_spatial_fir_; //2-d FloatArray
-  FloatArray agc_spatial_n_taps_; //FloatArray
-  FloatArray agc_mix_coeffs_; //FloatArray
+  FloatArray agc_epsilon_;
+  std::vector<int> decimation_;
+  FloatArray agc_pole_z1_;
+  FloatArray agc_pole_z2_;
+  FloatArray agc_spatial_iterations_;
+  FloatArray2d agc_spatial_fir_;
+  FloatArray agc_spatial_n_taps_;
+  FloatArray agc_mix_coeffs_;
   FPType agc_gain_;
   FPType detect_scale_;
-  
-  FloatArray agc1_scales_;
-  FloatArray agc2_scales_;
-  FloatArray time_constants_;
-  FPType tau_;
-  FPType decim_;
-  FPType n_times_;
-  FPType delay_;
-  FPType spread_sq_;
-  FPType u_;
-  FPType p_;
-  FPType dp_;
-  FPType pole_z1_;
-  FPType pole_z2_;
-  int n_taps_;
-  bool fir_ok_;
-  int n_iterations_;
-  FPType total_dc_gain_;
-  FPType a, b;
-  FPType mix_coeff_;
-  FloatArray fir_;
-
-  
-  void OutputCoeffs(); //Method to view coeffs, could go in final version
-  void DesignAGC(AGCParams agc_params, long fs, int n_ch);
 };
 
-
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/agc_params.cc
--- a/carfac/agc_params.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/agc_params.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,15 +22,16 @@
 
 #include "agc_params.h"
 
-AGCParams::AGCParams(){
+// 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;
   FPType agc1f = 1.0;
   FPType agc2f = 1.65;
   time_constants_.resize(n_stages_);
-  time_constants_ << 1*0.002, 4*0.002, 16*0.002, 64*0.002;
-  decimation_.resize(n_stages_);
-  decimation_ << 8, 2, 2, 2;
+  time_constants_ << 1 * 0.002, 4 * 0.002, 16 * 0.002, 64 * 0.002;
+  decimation_ = {8, 2, 2, 2};
   agc1_scales_.resize(n_stages_);
   agc1_scales_ << 1.0 * agc1f, 1.4 * agc1f, 2.0 * agc1f, 2.8 * agc1f;
   agc2_scales_.resize(n_stages_);
@@ -38,8 +39,10 @@
   agc_mix_coeff_ = 0.5;
 }
 
-void AGCParams::SetParams(int ns, FPType agcsg, FPType agcmc, FloatArray tc,
-                          FloatArray dec, FloatArray agc1sc, FloatArray agc2sc){
+// The overloaded constructor allows for use of different AGC parameters.
+AGCParams::AGCParams(int ns, FPType agcsg, FPType agcmc, FloatArray tc,
+                     std::vector<int> dec, FloatArray agc1sc,
+                          FloatArray agc2sc) {
   n_stages_ = ns;
   agc_stage_gain_ = agcsg;
   agc_mix_coeff_ = agcmc;
@@ -47,16 +50,4 @@
   decimation_ = dec;
   agc1_scales_ = agc1sc;
   agc2_scales_ = agc2sc;
-}
-
-void AGCParams::OutputParams(){
-  std::cout << "AGCParams Values" << std::endl;
-  std::cout << "****************" << std::endl;
-  std::cout << "n_stages_ = " << n_stages_ << std::endl;
-  std::cout << "agc_stage_gain_ = " << agc_stage_gain_ << std::endl;
-  std::cout << "agc_mix_coeff_ = " << agc_mix_coeff_ << std::endl;
-  std::cout << "time_constants_ = " << time_constants_ << std::endl;
-  std::cout << "decimation_ = " << decimation_ << std::endl;
-  std::cout << "agc1_scales_ = " << agc1_scales_ << std::endl;
-  std::cout << "agc2_scales_ = " << agc2_scales_ << std::endl << std::endl;
 }
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/agc_params.h
--- a/carfac/agc_params.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/agc_params.h	Thu May 16 17:33:23 2013 +0000
@@ -25,21 +25,17 @@
 
 #include "carfac_common.h"
 
-class AGCParams {
-public:
+struct AGCParams {
+  AGCParams();
+  AGCParams(int ns, FPType agcsg, FPType agcmc, FloatArray tc,
+            std::vector<int> dec, FloatArray agc1sc, FloatArray agc2sc);
   int n_stages_;
   FPType agc_stage_gain_;
   FPType agc_mix_coeff_;
   FloatArray time_constants_;
-  FloatArray decimation_;
+  std::vector<int> decimation_;
   FloatArray agc1_scales_;
   FloatArray agc2_scales_;
-  
-  
-  void OutputParams();
-  void SetParams(int ns, FPType agcsg,FPType agcmc, FloatArray tc,
-                 FloatArray dec, FloatArray agc1sc, FloatArray agc2sc);
-  AGCParams();
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/agc_state.cc
--- a/carfac/agc_state.cc	Mon May 13 22:51:15 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-//
-//  agc_state.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_state.h"
-
-void AGCState::InitAGCState(AGCCoeffs agc_coeffs){
-  n_ch_ = agc_coeffs.n_ch_;
-  n_agc_stages_ = agc_coeffs.n_agc_stages_;
-  agc_memory_ = FloatArray2d::Zero(n_ch_, n_agc_stages_);
-  input_accum_ = FloatArray2d::Zero(n_ch_, n_agc_stages_);
-  decim_phase_ = FloatArray::Zero(n_agc_stages_);
-}
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/agc_state.h
--- a/carfac/agc_state.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/agc_state.h	Thu May 16 17:33:23 2013 +0000
@@ -25,14 +25,12 @@
 
 #include "agc_coeffs.h"
 
-class AGCState {
-public:
+struct AGCState {
   int n_ch_;
   int n_agc_stages_;
-  FloatArray2d agc_memory_;
-  FloatArray2d input_accum_;
-  FloatArray decim_phase_;
-  void InitAGCState(AGCCoeffs agc_coeffs);
+  std::vector<FloatArray> agc_memory_;
+  std::vector<FloatArray> input_accum_;
+  std::vector<int> decim_phase_;
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/car_coeffs.cc
--- a/carfac/car_coeffs.cc	Mon May 13 22:51:15 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,77 +0,0 @@
-//
-//  car_coeffs.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_coeffs.h"
-
-//This method has been created for debugging purposes and depends on <iostream>.
-//Could possibly be removed in the final version to reduce dependencies. 
-void CARCoeffs::OutputCoeffs(){
-  std::cout << "CARCoeffs Values" << std::endl;
-  std::cout << "****************" << std::endl;
-  std::cout << "n_ch_ = " << n_ch_ << std::endl;
-  std::cout << "velocity_scale_ = " << velocity_scale_ << std::endl;
-  std::cout << "v_offset_ = " << v_offset_ << std::endl;
-  std::cout << "r1_coeffs_ = " << r1_coeffs_ << std::endl;
-  std::cout << "a0_coeffs_ = " << a0_coeffs_ << std::endl;
-  std::cout << "c0_coeffs_ = " << c0_coeffs_ << std::endl;
-  std::cout << "h_coeffs_ = " << h_coeffs_ << std::endl;
-  std::cout << "g0_coeffs_ = " << g0_coeffs_ << std::endl;
-  std::cout << "zr_coeffs_ = " << zr_coeffs_ << std::endl << std::endl;
-}
-
-void CARCoeffs::DesignFilters(CARParams car_params, long fs,
-                              FloatArray pole_freqs[]){
-  
-  n_ch_ = int(pole_freqs->size());
-  velocity_scale_ = car_params.velocity_scale_;
-  v_offset_ = car_params.v_offset_;
-  FloatArray p_freqs = *pole_freqs;
-  r1_coeffs_.resize(n_ch_);
-  a0_coeffs_.resize(n_ch_);
-  c0_coeffs_.resize(n_ch_);
-  h_coeffs_.resize(n_ch_);
-  g0_coeffs_.resize(n_ch_);
-  FPType f = car_params.zero_ratio_ * car_params.zero_ratio_ - 1;
-  FloatArray theta = p_freqs * ((2 * PI) / fs);
-  c0_coeffs_ = sin(theta);
-  a0_coeffs_ = cos(theta);
-  FPType ff = car_params.high_f_damping_compression_;
-  FloatArray x = theta / PI;
-  zr_coeffs_ = PI * (x - (ff * (x * x * x)));//change to exponet
-  FPType max_zeta = car_params.max_zeta_;
-  FPType min_zeta = car_params.min_zeta_;
-  r1_coeffs_ = (1 - (zr_coeffs_ * max_zeta));
-  FPType curfreq;
-  FloatArray erb_freqs(n_ch_);
-  for (int ch=0; ch < n_ch_; ch++){
-    curfreq = p_freqs(ch);
-    erb_freqs(ch) = ERBHz(curfreq, car_params.erb_break_freq_,
-                           car_params.erb_q_);
-  }
-  FloatArray min_zetas = min_zeta + (0.25 * ((erb_freqs / p_freqs) - min_zeta));
-  zr_coeffs_ = zr_coeffs_ * (max_zeta - min_zetas);
-  h_coeffs_ = c0_coeffs_ * f;
-  FloatArray relative_undamping = FloatArray::Ones(n_ch_);
-  FloatArray r = r1_coeffs_ + (zr_coeffs_ * relative_undamping);
-  g0_coeffs_ = (1 - (2 * r * a0_coeffs_) + (r * r)) /
-            (1 - (2 * r * a0_coeffs_) + (h_coeffs_ * r * c0_coeffs_) + (r * r));
-}
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/car_coeffs.h
--- a/carfac/car_coeffs.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/car_coeffs.h	Thu May 16 17:33:23 2013 +0000
@@ -25,9 +25,7 @@
 
 #include "car_params.h"
 
-class CARCoeffs {
-public:
-  int n_ch_;
+struct CARCoeffs {
   FPType velocity_scale_;
   FPType v_offset_;
   FloatArray r1_coeffs_;
@@ -36,11 +34,6 @@
   FloatArray h_coeffs_;
   FloatArray g0_coeffs_;
   FloatArray zr_coeffs_;
-  
-  void OutputCoeffs(); //Method to view coeffs, could go in final version
-  void DesignFilters(CARParams car_params, long fs,
-                     FloatArray pole_freqs[]);
 };
 
-
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/car_params.cc
--- a/carfac/car_params.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/car_params.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,40 +22,21 @@
 
 #include "car_params.h"
 
-CARParams::CARParams(){
-  velocity_scale_ = 0.1; //for the velocity nonlinearity
-  v_offset_ = 0.04; //offset gives quadratic part
-  min_zeta_ = 0.1; //minimum damping factor in mid-freq channels
-  max_zeta_ = 0.35; //maximum damping factor in mid-freq channels
+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_ = 1.4142; //how far zero is above pole
-  high_f_damping_compression_ = 0.5; //0 to 1 to compress theta
-  erb_per_step_ = 0.5; //assume G&M's ERB formula
+  zero_ratio_ = 1.4142;
+  high_f_damping_compression_ = 0.5;
+  erb_per_step_ = 0.5;
   min_pole_hz_ = 30;
-  erb_break_freq_ = 165.3; //Greenwood map's break frequency
+  erb_break_freq_ = 165.3;
   erb_q_ = 1000/(24.7*4.37);
 };
 
-//This method has been created for debugging purposes and depends on <iostream>.
-//Could possibly be removed in the final version to reduce dependencies.
-void CARParams::OutputParams(){
-  std::cout << "CARParams Values" << std::endl;
-  std::cout << "****************" << std::endl;
-  std::cout << "velocity_scale_ = " << velocity_scale_ << std::endl;
-  std::cout << "v_offset_ = " << v_offset_ << std::endl;
-  std::cout << "min_zeta_ = " << min_zeta_ << std::endl;
-  std::cout << "max_zeta_ = " << max_zeta_ << std::endl;
-  std::cout << "first_pole_theta_ = " << first_pole_theta_ << std::endl;
-  std::cout << "zero_ratio_ = " << zero_ratio_ << std::endl;
-  std::cout << "high_f_damping_compression_ = " << high_f_damping_compression_
-            << std::endl;
-  std::cout << "erb_per_step_ = " << erb_per_step_ << std::endl;
-  std::cout << "min_pole_hz_ = " << min_pole_hz_ << std::endl;
-  std::cout << "erb_break_freq_ = " << erb_break_freq_ << std::endl;
-  std::cout << "erb_q_ = " << erb_q_ << std::endl << std::endl;
-}
-
-void CARParams::SetParams(FPType vs, FPType voff, FPType min_z, FPType max_z,
+CARParams::CARParams(FPType vs, FPType voff, FPType min_z, FPType max_z,
                           FPType fpt, FPType zr, FPType hfdc, FPType erbps,
                           FPType mph, FPType erbbf, FPType erbq) {  
   velocity_scale_ = vs;
diff -r aefe2ca0674f -r 01986636257a carfac/car_params.h
--- a/carfac/car_params.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/car_params.h	Thu May 16 17:33:23 2013 +0000
@@ -25,26 +25,22 @@
 
 #include "carfac_common.h"
 
-class CARParams {
-public:
-  FPType velocity_scale_; //for the velocity nonlinearity
-  FPType v_offset_; //offset gives quadratic part
-  FPType min_zeta_; //minimum damping factor in mid-freq channels
-  FPType max_zeta_; //maximum damping factor in mid-freq channels
+struct CARParams {
+  CARParams();  // The constructor initializes using default parameter values.
+  CARParams(FPType vs, FPType voff, FPType min_z, FPType max_z, FPType fpt,
+                 FPType zr, FPType hfdc, FPType erbps, FPType mph, FPType erbbf,
+                 FPType erbq);  // This is a method to set non-default params.
+  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.
+  FPType max_zeta_;  // This is the maximum damping factor in mid-freq channels.
   FPType first_pole_theta_;
-  FPType zero_ratio_; //how far zero is above pole
-  FPType high_f_damping_compression_; //0 to 1 to compress theta
-  FPType erb_per_step_; //assume G&M's ERB formula
+  FPType zero_ratio_;  // This is how far zero is above the pole.
+  FPType high_f_damping_compression_;  // A range from 0 to 1 to compress theta.
+  FPType erb_per_step_;  // Here we assume G&M's ERB formula.
   FPType min_pole_hz_;
-  FPType erb_break_freq_; //Greenwood map's break frequency
-  FPType erb_q_; //G&M's high-cf ratio
-  
-  CARParams(); //Constructor initializes default parameter values
-  void OutputParams(); //Output current parameter values using std::cout
-  void SetParams(FPType vs, FPType voff, FPType min_z, FPType max_z, FPType fpt,
-                 FPType zr, FPType hfdc, FPType erbps, FPType mph, FPType erbbf,
-                 FPType erbq); //Method to set non-default params
+  FPType erb_break_freq_;  // This is the Greenwood map's break frequency.
+  FPType erb_q_;  // This represents Glassberg and Moore's high-cf ratio.
 };
 
-
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/car_state.cc
--- a/carfac/car_state.cc	Mon May 13 22:51:15 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,37 +0,0 @@
-//
-//  car_state.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_state.h"
-
-void CARState::InitCARState(CARCoeffs car_coeffs){
-  n_ch_ = car_coeffs.n_ch_;
-  z1_memory_ = FloatArray::Zero(n_ch_);
-  z2_memory_ = FloatArray::Zero(n_ch_);
-  za_memory_ = FloatArray::Zero(n_ch_);
-  zb_memory_ = car_coeffs.zr_coeffs_;
-  dzb_memory_ = FloatArray::Zero(n_ch_);
-  zy_memory_ = FloatArray::Zero(n_ch_);
-  g_memory_ = car_coeffs.g0_coeffs_;
-  dg_memory_ = FloatArray::Zero(n_ch_);
-  
-  
-}
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/car_state.h
--- a/carfac/car_state.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/car_state.h	Thu May 16 17:33:23 2013 +0000
@@ -25,9 +25,7 @@
 
 #include "car_coeffs.h"
 
-class CARState {
-public:
-  int n_ch_;
+struct CARState {
   FloatArray z1_memory_;
   FloatArray z2_memory_;
   FloatArray za_memory_;
@@ -36,10 +34,6 @@
   FloatArray zy_memory_;
   FloatArray g_memory_;
   FloatArray dg_memory_;
-  
-  void InitCARState(CARCoeffs car_coeffs);
-  
 };
 
-
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/carfac.cc
--- 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));
+  }
+}
diff -r aefe2ca0674f -r 01986636257a carfac/carfac.h
--- a/carfac/carfac.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/carfac.h	Thu May 16 17:33:23 2013 +0000
@@ -43,12 +43,7 @@
 #include "carfac_output.h"
 
 class CARFAC {
-public:
-  int n_ears_; //number of ears
-  long fs_; //sample rate
-  int n_ch_; //number of channels in the CARFAC model
-  Ear *ears_; //array of Ear objects for mono/stereo/multichannel processing
-  
+ public:
   // The 'Design' method takes a set of CAR, IHC and AGC parameters along with
   // arguments specifying the number of 'ears' (audio file channels) and sample
   // rate. This initializes a vector of 'Ear' objects -- one for mono, two for
@@ -57,15 +52,23 @@
   // model.
   void Design(int n_ears, long fs, CARParams car_params, IHCParams ihc_params,
               AGCParams agc_params);
+  // The 'Run' method processes an entire file with the current model, using
+  // subsequent calls to the 'RunSegment' method
+  CARFACOutput Run(FloatArray2d sound_data, bool open_loop);
+  // The 'RunSegment' method processes individual sound segments
+  void RunSegment(FloatArray2d sound_data, CARFACOutput *seg_output,
+                  bool open_loop);
   
-  //The 'Run' method processes an entire file with the current model, using
-  //subsequent calls to the 'RunSegment' method
-  CARFACOutput Run(FloatArray2d sound_data);
+ private:
+  void CrossCouple();
+  void CloseAGCLoop();
   
-  //The 'RunSegment' method processes individual sound segments
-  void RunSegment(FloatArray2d sound_data, CARFACOutput *seg_output);
+  int n_ears_;  // This is the number of ears.
+  long fs_;  // This is our current sample rate.
+  int n_ch_;  // This is the number of channels in the CARFAC model.
+  FPType max_channels_per_octave_;
+  // We store an array of Ear objects for mono/stereo/multichannel processing:
+  Ear *ears_; 
 };
 
-
-
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/carfac_common.cc
--- a/carfac/carfac_common.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/carfac_common.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,39 +22,22 @@
 
 #include "carfac_common.h"
 
-//Auditory filter nominal Equivalent Rectangular Bandwidth
-//Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138
+// Auditory filter nominal Equivalent Rectangular Bandwidth
+// Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138
 FPType ERBHz (FPType cf_hz, FPType erb_break_freq, FPType erb_q) {
-  
   FPType erb;
   erb = (erb_break_freq + cf_hz) / erb_q;
   return erb;
 }
 
-//An IHC-like sigmoidal detection nonlinearity for the CARFAC.
-//Resulting conductance is in about [0...1.3405]
-FPType CARFACDetect (FPType x) {
-  
-  FPType conductance, z;
+FloatArray CARFACDetect (FloatArray x) {
+  FloatArray conductance, z, set;
   FPType a = 0.175;
-  //offset of low-end tail into neg x territory
-  //this parameter is adjusted for the book, to make the 20% DC response
-  //threshold at 0.1
+  // This offsets the low-end tail into negative x territory.
+  // The parameter is adjusted for the book, to make the 20% DC response
+  // threshold at 0.1.
   z  = x + a;
-  conductance = pow(z,3) / (pow(z,3) + pow(z,2) + 0.1);
-  //zero is the final answer for many points:
-  return conductance;
-}
-
-FloatArray CARFACDetect (FloatArray x) {
-  
-  FloatArray conductance, z;
-  FPType a = 0.175;
-  //offset of low-end tail into neg x territory
-  //this parameter is adjusted for the book, to make the 20% DC response
-  //threshold at 0.1
-  z  = x + a;
-  conductance = (z * z * z) / ((z * z * z) + (z * z) + 0.1);
-  //zero is the final answer for many points:
+  // Zero is the final answer for many points.
+  conductance = (z < 0).select(0.0, (z * z * z) / (z * z * z + z * z + 0.1));
   return conductance;
 }
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/carfac_common.h
--- a/carfac/carfac_common.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/carfac_common.h	Thu May 16 17:33:23 2013 +0000
@@ -47,17 +47,34 @@
 #ifndef CARFAC_Open_Source_C__Library_CARFACCommon_h
 #define CARFAC_Open_Source_C__Library_CARFACCommon_h
 
-#include <iostream> //Used for debugging output, could go in final version
-#include <math.h> //Used during coefficient calculations
-#include <Eigen/Dense> //Used for 1d and 2d floating point arrays
+// This section is where the base include operations for the CARFAC project
+// occur.
+// <iostream> is used for debugging output, but it could go in final version.
+#include <iostream>
+// <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
+#include <vector>
+// The Eigen library is used extensively for 1d and 2d floating point arrays.
+// For more information, see: http://eigen.tuxfamily.org
+#include <Eigen/Dense>
 using namespace Eigen;
 
-#define PI 3.141592
-typedef float FPType; //Used to enable easy switching in precision level
+// 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
 
-//typedefs for Eigen floating point arrays
-typedef Eigen::Array<FPType,Dynamic,1> FloatArray; //1d floating point array
-typedef Eigen::Array<FPType,Dynamic,Dynamic> FloatArray2d; //2d fpoint array
+// Three typedefs are used for enabling quick switching of precision and array
+// usage. 
+// The 'FPType' typedef is used to enable easy switching in precision level.
+typedef double FPType;
+// These are the two typedefs for Eigen floating point arrays.
+typedef Eigen::Array<FPType, Dynamic, 1> FloatArray;  // This is a 1d array.
+typedef Eigen::Array<FPType, Dynamic, Dynamic> FloatArray2d;  // This is 2d.
+
+// Two psychoacoustics helper functions are defined here for use by the
+// different processing stages in calculating coeffecients.
 
 // Function: ERBHz
 // Auditory filter nominal Equivalent Rectangular Bandwidth
@@ -66,6 +83,6 @@
 
 // Function CARFACDetect
 // TODO explain a bit more
-FPType CARFACDetect (FPType x);
 FloatArray CARFACDetect (FloatArray x);
-#endif
+
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/carfac_output.cc
--- a/carfac/carfac_output.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/carfac_output.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,15 +22,15 @@
 
 #include "carfac_output.h"
 
-void CARFACOutput::InitOutput(int n_ears, int n_ch, long n_tp){
+void CARFACOutput::InitOutput(int n_ears, int n_ch, long n_tp) {
   n_ears_ = n_ears;
   ears_ = new EarOutput[n_ears_];
-  for (int i = 0; i < n_ears_; i++){
-    ears_[i].InitOutput(n_ch,n_tp);
+  for (int i = 0; i < n_ears_; i++) {
+    ears_[i].InitOutput(n_ch, n_tp);
   }
 }
 
-void CARFACOutput::MergeOutput(CARFACOutput output, long start, long length){
+void CARFACOutput::MergeOutput(CARFACOutput output, long start, long length) {
   for (int i = 0; i < n_ears_; i++){
     ears_[i].MergeOutput(output.ears_[i], start, length);
   }
diff -r aefe2ca0674f -r 01986636257a carfac/carfac_output.h
--- a/carfac/carfac_output.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/carfac_output.h	Thu May 16 17:33:23 2013 +0000
@@ -43,12 +43,11 @@
 
 #include "ear_output.h"
 
-class CARFACOutput {
-public:
+struct CARFACOutput {
+  void InitOutput(int n_ears, int n_ch, long n_tp);
+  void MergeOutput(CARFACOutput output, long start, long length);
   int n_ears_;
   EarOutput *ears_;
-  void InitOutput(int n_ears, int n_ch, long n_tp);
-  void MergeOutput(CARFACOutput output, long start, long length);
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/ear.cc
--- a/carfac/ear.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ear.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,35 +22,263 @@
 
 #include "ear.h"
 
-void Ear::InitEar(long fs, CARParams car_p, IHCParams ihc_p, AGCParams agc_p){
-  car_params_ = car_p;
-  ihc_params_ = ihc_p;
-  agc_params_ = agc_p;
-  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_);
-  }
-  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));
-  car_coeffs_.DesignFilters(car_params_, fs, &pole_freqs);
-  agc_coeffs_.DesignAGC(agc_params_, fs, n_ch_);
-  ihc_coeffs_.DesignIHC(ihc_params_, fs, n_ch_);
-  car_state_.InitCARState(car_coeffs_);
-  agc_state_.InitAGCState(agc_coeffs_);
-  ihc_state_.InitIHCState(ihc_coeffs_);
-
+// 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::InitEar(int n_ch, long fs, FloatArray pole_freqs,
+                  CARParams car_p, IHCParams ihc_p, AGCParams agc_p) {
+  // 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.
+  n_ch_ = n_ch;
+  // These functions use the parameters that have been passed to design the
+  // coefficients for the three stages of the model.
+  DesignFilters(car_p, fs, pole_freqs);
+  DesignIHC(ihc_p, fs);
+  DesignAGC(agc_p, fs);
+  // Once the coefficients have been determined, we can initialize the state
+  // variables that will be used during runtime.
+  InitCARState();
+  InitIHCState();
+  InitAGCState();
 }
 
-FloatArray Ear::CARStep(FPType input){
+void Ear::DesignFilters(CARParams car_params, long fs,
+                              FloatArray pole_freqs) {
+  car_coeffs_.velocity_scale_ = car_params.velocity_scale_;
+  car_coeffs_.v_offset_ = car_params.v_offset_;
+  car_coeffs_.r1_coeffs_.resize(n_ch_);
+  car_coeffs_.a0_coeffs_.resize(n_ch_);
+  car_coeffs_.c0_coeffs_.resize(n_ch_);
+  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;
+  FloatArray theta = pole_freqs * ((2 * PI) / fs);
+  car_coeffs_.c0_coeffs_ = sin(theta);
+  car_coeffs_.a0_coeffs_ = cos(theta);
+  FPType ff = car_params.high_f_damping_compression_;
+  FloatArray x = theta / PI;
+  car_coeffs_.zr_coeffs_ = PI * (x - (ff * (x * x * x)));
+  FPType max_zeta = car_params.max_zeta_;
+  FPType min_zeta = car_params.min_zeta_;
+  car_coeffs_.r1_coeffs_ = (1 - (car_coeffs_.zr_coeffs_ * max_zeta));
+  FPType curfreq;
+  FloatArray erb_freqs(n_ch_);
+  for (int ch=0; ch < n_ch_; ch++) {
+    curfreq = pole_freqs(ch);
+    erb_freqs(ch) = ERBHz(curfreq, car_params.erb_break_freq_,
+                          car_params.erb_q_);
+  }
+  FloatArray min_zetas = min_zeta + (0.25 * ((erb_freqs / pole_freqs) - min_zeta));
+  car_coeffs_.zr_coeffs_ *= max_zeta - min_zetas;
+  car_coeffs_.h_coeffs_ = car_coeffs_.c0_coeffs_ * f;
+  FloatArray relative_undamping = FloatArray::Ones(n_ch_);
+  FloatArray r = car_coeffs_.r1_coeffs_ +
+                 (car_coeffs_.zr_coeffs_ * relative_undamping);
+  car_coeffs_.g0_coeffs_ = (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));
+}
+
+
+void Ear::DesignAGC(AGCParams agc_params, long fs) {
+  // These data members could probably be initialized locally within the design
+  // function.
+  agc_coeffs_.n_agc_stages_ = agc_params.n_stages_;
+  agc_coeffs_.agc_stage_gain_ = agc_params.agc_stage_gain_;
+  FloatArray agc1_scales = agc_params.agc1_scales_;
+  FloatArray agc2_scales = agc_params.agc2_scales_;
+  FloatArray time_constants = agc_params.time_constants_;
+  agc_coeffs_.agc_epsilon_.resize(agc_coeffs_.n_agc_stages_);
+  agc_coeffs_.agc_pole_z1_.resize(agc_coeffs_.n_agc_stages_);
+  agc_coeffs_.agc_pole_z2_.resize(agc_coeffs_.n_agc_stages_);
+  agc_coeffs_.agc_spatial_iterations_.resize(agc_coeffs_.n_agc_stages_);
+  agc_coeffs_.agc_spatial_n_taps_.resize(agc_coeffs_.n_agc_stages_);
+  agc_coeffs_.agc_spatial_fir_.resize(3,agc_coeffs_.n_agc_stages_);
+  agc_coeffs_.agc_mix_coeffs_.resize(agc_coeffs_.n_agc_stages_);
+  FPType mix_coeff = agc_params.agc_mix_coeff_;
+  int decim = 1;
+  agc_coeffs_.decimation_ = agc_params.decimation_;
+  FPType total_dc_gain = 0;
+  // Here we loop through each of the stages of the AGC.
+  for (int stage=0; stage < agc_coeffs_.n_agc_stages_; stage++) {
+    FPType tau = time_constants(stage);
+    decim *= agc_coeffs_.decimation_.at(stage);
+    agc_coeffs_.agc_epsilon_(stage) = 1 - exp((-1 * decim) / (tau * fs));
+    FPType n_times = tau * (fs / decim);
+    FPType delay = (agc2_scales(stage) - agc1_scales(stage)) / n_times;
+    FPType spread_sq = (pow(agc1_scales(stage), 2) +
+                        pow(agc2_scales(stage), 2)) / n_times;
+    FPType u = 1 + (1 / spread_sq);
+    FPType p = u - sqrt(pow(u, 2) - 1);
+    FPType dp = delay * (1 - (2 * p) + (p * p)) / 2;
+    FPType pole_z1 = p - dp;
+    FPType pole_z2 = p + dp;
+    agc_coeffs_.agc_pole_z1_(stage) = pole_z1;
+    agc_coeffs_.agc_pole_z2_(stage) = pole_z2;
+    int n_taps = 0;
+    bool fir_ok = false;
+    int n_iterations = 1;
+    // This section initializes the FIR coeffs settings at each stage.
+    FloatArray fir(3);
+    while (! fir_ok) {
+      switch (n_taps) {
+        case 0:
+          n_taps = 3;
+          break;
+        case 3:
+          n_taps = 5;
+          break;
+        case 5:
+          n_iterations ++;
+          if (n_iterations > 16){
+            // This implies too many iterations, so we shoud indicate and error.
+            // TODO alexbrandmeyer, proper Google error handling.
+          }
+          break;
+        default:
+          // This means a bad n_taps has been provided, so there should again be
+          // an error.
+          // TODO alexbrandmeyer, proper Google error handling.
+          break;
+      }
+      // Now we can design the FIR coefficients.
+      FPType var = spread_sq / n_iterations;
+      FPType mn = delay / n_iterations;
+      FPType a, b;
+      switch (n_taps) {
+        case 3:
+          a = (var + pow(mn, 2) - mn) / 2;
+          b = (var + pow(mn, 2) + mn) / 2;
+          fir << a, 1 - a - b, b;
+          fir_ok = (fir(2) >= 0.2) ? true : false;
+          break;
+        case 5:
+          a = (((var + pow(mn, 2)) * 2/5) - (mn * 2/3)) / 2;
+          b = (((var + pow(mn, 2)) * 2/5) + (mn * 2/3)) / 2;
+          fir << a/2, 1 - a - b, b/2;
+          fir_ok = (fir(2) >= 0.1) ? true : false;
+          break;
+        default:
+          break;  // Again, we've arrived at a bad n_taps in the design.
+          // TODO alexbrandmeyer. Add proper Google error handling.
+      }
+    }
+    // Once we have the FIR design for this stage we can assign it to the
+    // appropriate data members.
+    agc_coeffs_.agc_spatial_iterations_(stage) = n_iterations;
+    agc_coeffs_.agc_spatial_n_taps_(stage) = n_taps;
+    // TODO alexbrandmeyer replace using Eigen block method
+    agc_coeffs_.agc_spatial_fir_(0, stage) = fir(0);
+    agc_coeffs_.agc_spatial_fir_(1, stage) = fir(1);
+    agc_coeffs_.agc_spatial_fir_(2, stage) = fir(2);
+    total_dc_gain += pow(agc_coeffs_.agc_stage_gain_, (stage));
+    agc_coeffs_.agc_mix_coeffs_(stage) =
+      (stage == 0) ? 0 : mix_coeff / (tau * (fs /decim));
+  }
+  agc_coeffs_.agc_gain_ = total_dc_gain;
+  agc_coeffs_.detect_scale_ = 1 / total_dc_gain;
+}
+
+void Ear::DesignIHC(IHCParams ihc_params, long fs) {
+  if (ihc_params.just_hwr_) {
+    ihc_coeffs_.just_hwr_ = ihc_params.just_hwr_;
+  } else {
+    // This section calculates conductance values using two pre-defined scalars.
+    FloatArray x(1);
+    FPType conduct_at_10, conduct_at_0;
+    x << 10;
+    x = CARFACDetect(x);
+    conduct_at_10 = x(0);
+    x << 0;
+    x = CARFACDetect(x);
+    conduct_at_0 = x(0);
+    if (ihc_params.one_cap_) {
+      FPType ro = 1 / conduct_at_10 ;
+      FPType c = ihc_params.tau1_out_ / ro;
+      FPType ri = ihc_params.tau1_in_ / c;
+      FPType saturation_output = 1 / ((2 * ro) + ri);
+      FPType r0 = 1 / conduct_at_0;
+      FPType current = 1 / (ri + r0);
+      ihc_coeffs_.cap1_voltage_ = 1 - (current * ri);
+      ihc_coeffs_.just_hwr_ = false;
+      ihc_coeffs_.lpf_coeff_ = 1 - exp( -1 / (ihc_params.tau_lpf_ * fs));
+      ihc_coeffs_.out1_rate_ = ro / (ihc_params.tau1_out_ * fs);
+      ihc_coeffs_.in1_rate_ = 1 / (ihc_params.tau1_in_ * fs);
+      ihc_coeffs_.one_cap_ = ihc_params.one_cap_;
+      ihc_coeffs_.output_gain_ = 1 / (saturation_output - current);
+      ihc_coeffs_.rest_output_ = current / (saturation_output - current);
+      ihc_coeffs_.rest_cap1_ = ihc_coeffs_.cap1_voltage_;
+    } else {
+      FPType ro = 1 / conduct_at_10;
+      FPType c2 = ihc_params.tau2_out_ / ro;
+      FPType r2 = ihc_params.tau2_in_ / c2;
+      FPType c1 = ihc_params.tau1_out_ / r2;
+      FPType r1 = ihc_params.tau1_in_ / c1;
+      FPType saturation_output = 1 / (2 * ro + r2 + r1);
+      FPType r0 = 1 / conduct_at_0;
+      FPType current = 1 / (r1 + r2 + r0);
+      ihc_coeffs_.cap1_voltage_ = 1 - (current * r1);
+      ihc_coeffs_.cap2_voltage_ = ihc_coeffs_.cap1_voltage_ - (current * r2);
+      ihc_coeffs_.just_hwr_ = false;
+      ihc_coeffs_.lpf_coeff_ = 1 - exp(-1 / (ihc_params.tau_lpf_ * fs));
+      ihc_coeffs_.out1_rate_ = 1 / (ihc_params.tau1_out_ * fs);
+      ihc_coeffs_.in1_rate_ = 1 / (ihc_params.tau1_in_ * fs);
+      ihc_coeffs_.out2_rate_ = ro / (ihc_params.tau2_out_ * fs);
+      ihc_coeffs_.in2_rate_ = 1 / (ihc_params.tau2_in_ * fs);
+      ihc_coeffs_.one_cap_ = ihc_params.one_cap_;
+      ihc_coeffs_.output_gain_ = 1 / (saturation_output - current);
+      ihc_coeffs_.rest_output_ = current / (saturation_output - current);
+      ihc_coeffs_.rest_cap1_ = ihc_coeffs_.cap1_voltage_;
+      ihc_coeffs_.rest_cap2_ = ihc_coeffs_.cap2_voltage_;
+    }
+  }
+  ihc_coeffs_.ac_coeff_ = 2 * PI * ihc_params.ac_corner_hz_ / fs;
+}
+
+void Ear::InitCARState() {
+  car_state_.z1_memory_ = FloatArray::Zero(n_ch_);
+  car_state_.z2_memory_ = FloatArray::Zero(n_ch_);
+  car_state_.za_memory_ = FloatArray::Zero(n_ch_);
+  car_state_.zb_memory_ = car_coeffs_.zr_coeffs_;
+  car_state_.dzb_memory_ = FloatArray::Zero(n_ch_);
+  car_state_.zy_memory_ = FloatArray::Zero(n_ch_);
+  car_state_.g_memory_ = car_coeffs_.g0_coeffs_;
+  car_state_.dg_memory_ = FloatArray::Zero(n_ch_);
+}
+
+void Ear::InitIHCState() {
+  ihc_state_.ihc_accum_ = FloatArray::Zero(n_ch_);
+  if (! ihc_coeffs_.just_hwr_) {
+    ihc_state_.ac_coupler_ = FloatArray::Zero(n_ch_);
+    ihc_state_.lpf1_state_ = ihc_coeffs_.rest_output_ * FloatArray::Ones(n_ch_);
+    ihc_state_.lpf2_state_ = ihc_coeffs_.rest_output_ * FloatArray::Ones(n_ch_);
+    if (ihc_coeffs_.one_cap_) {
+      ihc_state_.cap1_voltage_ = ihc_coeffs_.rest_cap1_ *
+                                 FloatArray::Ones(n_ch_);
+    } else {
+      ihc_state_.cap1_voltage_ = ihc_coeffs_.rest_cap1_ *
+                                 FloatArray::Ones(n_ch_);
+      ihc_state_.cap2_voltage_ = ihc_coeffs_.rest_cap2_ *
+                                 FloatArray::Ones(n_ch_);
+    }
+  }
+}
+
+void Ear::InitAGCState() {
+  agc_state_.n_agc_stages_ = agc_coeffs_.n_agc_stages_;
+  agc_state_.agc_memory_.resize(agc_state_.n_agc_stages_);
+  agc_state_.input_accum_.resize(agc_state_.n_agc_stages_);
+  agc_state_.decim_phase_.resize(agc_state_.n_agc_stages_);
+  for (int i = 0; i < agc_state_.n_agc_stages_; i++) {
+    agc_state_.decim_phase_.at(i) = 0;
+    agc_state_.agc_memory_.at(i) = FloatArray::Zero(n_ch_);
+    agc_state_.input_accum_.at(i) = FloatArray::Zero(n_ch_);
+  }
+}
+
+FloatArray Ear::CARStep(FPType input) {
   FloatArray g(n_ch_);
   FloatArray zb(n_ch_);
   FloatArray za(n_ch_);
@@ -61,28 +289,30 @@
   FloatArray z2(n_ch_);
   FloatArray zy(n_ch_);
   FPType in_out;
-  
-  // do the DOHC stuff:
-  g = car_state_.g_memory_ + car_state_.dg_memory_; //interp g
-  zb = car_state_.zb_memory_ + car_state_.dzb_memory_; //AGC interpolation state
-  // update the nonlinear function of "velocity", and zA (delay of z2):
+  // This performs the DOHC stuff.
+  g = car_state_.g_memory_ + car_state_.dg_memory_;  // This interpolates g.
+  // This calculates the AGC interpolation state.
+  zb = car_state_.zb_memory_ + car_state_.dzb_memory_;
+  // This updates the nonlinear function of 'velocity' along with zA, which is
+  // a delay of z2.
   za = car_state_.za_memory_;
   v = car_state_.z2_memory_ - za;
   nlf = OHC_NLF(v);
-  r = car_coeffs_.r1_coeffs_ + (zb * nlf); // zB * nfl is "undamping" delta r
+  // Here, zb * nfl is "undamping" delta r.
+  r = car_coeffs_.r1_coeffs_ + (zb * nlf);
   za = car_state_.z2_memory_;
-  // now reduce state by r and rotate with the fixed cos/sin coeffs:
+  // Here we reduce the CAR state by r and rotate with the fixed cos/sin coeffs.
   z1 = r * ((car_coeffs_.a0_coeffs_ * car_state_.z1_memory_) -
             (car_coeffs_.c0_coeffs_ * car_state_.z2_memory_));
   z2 = r * ((car_coeffs_.c0_coeffs_ * car_state_.z1_memory_) +
             (car_coeffs_.a0_coeffs_ * car_state_.z2_memory_));
   zy = car_coeffs_.h_coeffs_ * z2;
-  // Ripple input-output path, instead of parallel, to avoid delay...
-  // this is the only part that doesn't get computed "in parallel":
+  // This section ripples the input-output path, to avoid delay...
+  // It's the only part that doesn't get computed "in parallel":
   in_out = input;
-  for (int ch = 0; ch < n_ch_; ch++){
+  for (int ch = 0; ch < n_ch_; ch++) {
     z1(ch) = z1(ch) + in_out;
-    // ripple, saving final channel outputs in zY
+    // This performs the ripple, and saves the final channel outputs in zy.
     in_out = g(ch) * (in_out + zy(ch));
     zy(ch) = in_out;
   }
@@ -92,23 +322,24 @@
   car_state_.zb_memory_ = zb;
   car_state_.zy_memory_ = zy;
   car_state_.g_memory_ = g;
-  // car_out is equal to zy state;
+  // The output of the CAR is equal to the zy state.
   return zy;
 }
 
-// start with a quadratic nonlinear function, and limit it via a
-// rational function; make the result go to zero at high
+// We start with a quadratic nonlinear function, and limit it via a
+// rational function. This makes the result go to zero at high
 // absolute velocities, so it will do nothing there.
-FloatArray Ear::OHC_NLF(FloatArray velocities){
-  FloatArray nlf(n_ch_);
-  nlf = 1 / ((velocities * car_coeffs_.velocity_scale_) +
-             (car_coeffs_.v_offset_ * car_coeffs_.v_offset_));
-  return nlf;
+FloatArray Ear::OHC_NLF(FloatArray velocities) {
+  return 1 / (1 +
+        (((velocities * car_coeffs_.velocity_scale_) + car_coeffs_.v_offset_) *
+         ((velocities * car_coeffs_.velocity_scale_) + car_coeffs_.v_offset_)));
+
 }
 
-// One sample-time update of inner-hair-cell (IHC) model, including the
-// detection nonlinearity and one or two capacitor state variables.
-FloatArray Ear::IHCStep(FloatArray car_out){
+// 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.
+FloatArray Ear::IHCStep(FloatArray car_out) {
   FloatArray ihc_out(n_ch_);
   FloatArray ac_diff(n_ch_);
   FloatArray conductance(n_ch_);
@@ -116,21 +347,11 @@
   ihc_state_.ac_coupler_ = ihc_state_.ac_coupler_ +
                            (ihc_coeffs_.ac_coeff_ * ac_diff);
   if (ihc_coeffs_.just_hwr_) {
-    //TODO Figure out best implementation with Eigen max/min methods
-    for (int ch = 0; ch < n_ch_; ch++){
+    for (int ch = 0; ch < n_ch_; ch++) {
       FPType a;
-      if (ac_diff(ch) > 0){
-        a = ac_diff(ch);
-      } else {
-        a = 0;
-      }
-      if (a < 2){
-        ihc_out(ch) = a;
-      } else {
-        ihc_out(ch) = 2;
-      }
+      a = (ac_diff(ch) > 0) ? ac_diff(ch) : 0;
+      ihc_out(ch) = (a < 2) ? a : 2;
     }
-    
   } else {
     conductance = CARFACDetect(ac_diff);
     if (ihc_coeffs_.one_cap_) {
@@ -150,8 +371,8 @@
                           ((ihc_state_.cap1_voltage_ - ihc_state_.cap2_voltage_)
                           * ihc_coeffs_.in2_rate_);
     }
-    // smooth it twice with LPF:
-    ihc_out = ihc_out * ihc_coeffs_.output_gain_;
+    // Here we smooth the output twice using a LPF.
+    ihc_out *= ihc_coeffs_.output_gain_;
     ihc_state_.lpf1_state_ = ihc_state_.lpf1_state_ +
                   (ihc_coeffs_.lpf_coeff_ * (ihc_out - ihc_state_.lpf1_state_));
     ihc_state_.lpf2_state_ = ihc_state_.lpf2_state_ +
@@ -163,7 +384,7 @@
   return ihc_out;
 }
 
-bool Ear::AGCStep(FloatArray ihc_out){
+bool Ear::AGCStep(FloatArray ihc_out) {
   int stage = 0;
   FloatArray agc_in(n_ch_);
   agc_in = agc_coeffs_.detect_scale_ * ihc_out;
@@ -171,68 +392,61 @@
   return updated;
 }
 
-bool Ear::AGCRecurse(int stage, FloatArray agc_in){
+bool Ear::AGCRecurse(int stage, FloatArray agc_in) {
   bool updated = true;
-  // decim factor for this stage, relative to input or prev. stage:
-  int decim = agc_coeffs_.decimation_(stage);
-  // decim phase of this stage (do work on phase 0 only):
-  //TODO FIX MODULO
-  
-  int decim_phase = agc_state_.decim_phase_(stage);
+  // This is the decim factor for this stage, relative to input or prev. stage:
+  int decim = agc_coeffs_.decimation_.at(stage);
+  // This is the decim phase of this stage (do work on phase 0 only):
+  int decim_phase = agc_state_.decim_phase_.at(stage);
   decim_phase = decim_phase % decim;
-  agc_state_.decim_phase_(stage) = decim_phase;
-  // accumulate input for this stage from detect or previous stage:
-  agc_state_.input_accum_.block(0,stage,n_ch_,1) =
-                        agc_state_.input_accum_.block(0,stage,n_ch_,1) + agc_in;
-  
-  // nothing else to do if it's not the right decim_phase
-  if (decim_phase == 0){
-    // do lots of work, at decimated rate.
-    // decimated inputs for this stage, and to be decimated more for next:
-    agc_in = agc_state_.input_accum_.block(0,stage,n_ch_,1) / decim;
-    // reset accumulator:
-    agc_state_.input_accum_.block(0,stage,n_ch_,1) = FloatArray::Zero(n_ch_);
-    
-    if (stage < (agc_coeffs_.decimation_.size() - 1)){
-      // recurse to evaluate next stage(s)
-      updated = AGCRecurse(stage+1, agc_in);
-      // and add its output to this stage input, whether it updated or not:
-      agc_in = agc_in + (agc_coeffs_.agc_stage_gain_ *
-                         agc_state_.agc_memory_.block(0,stage+1,n_ch_,1));
+  agc_state_.decim_phase_.at(stage) = decim_phase;
+  // Here we accumulate input for this stage from the previous stage:
+  agc_state_.input_accum_.at(stage) += agc_in;
+  // We don't do anything if it's not the right decim_phase.
+  if (decim_phase == 0) {
+    // Now we do lots of work, at the decimated rate.
+    // These are the decimated inputs for this stage, which will be further
+    // decimated at the next stage.
+    agc_in = agc_state_.input_accum_.at(stage) / decim;
+    // This resets the accumulator.
+    agc_state_.input_accum_.at(stage) = FloatArray::Zero(n_ch_);
+    if (stage < (agc_coeffs_.decimation_.size() - 1)) {
+      // Now we recurse to evaluate the next stage(s).
+      updated = AGCRecurse(stage + 1, agc_in);
+      // Afterwards we add its output to this stage input, whether it updated or
+      // not.
+      agc_in += (agc_coeffs_.agc_stage_gain_ *
+                         agc_state_.agc_memory_.at(stage + 1));
     }
-    FloatArray agc_stage_state = agc_state_.agc_memory_.block(0,stage,n_ch_,1);
-    // first-order recursive smoothing filter update, in time:
+    FloatArray agc_stage_state = agc_state_.agc_memory_.at(stage);
+    // This performs a first-order recursive smoothing filter update, in time.
     agc_stage_state = agc_stage_state + (agc_coeffs_.agc_epsilon_(stage) *
                                          (agc_in - agc_stage_state));
     agc_stage_state = AGCSpatialSmooth(stage, agc_stage_state);
-    agc_state_.agc_memory_.block(0,stage,n_ch_,1) = agc_stage_state;
+    agc_state_.agc_memory_.at(stage) = agc_stage_state;
   } else {
     updated = false;
   }
   return updated;
 }
 
-FloatArray Ear::AGCSpatialSmooth(int stage, FloatArray stage_state){
+FloatArray Ear::AGCSpatialSmooth(int stage, FloatArray stage_state) {
   int n_iterations = agc_coeffs_.agc_spatial_iterations_(stage);
   bool use_fir;
-  if (n_iterations < 4){
-    use_fir = true;
-  } else {
-    use_fir = false;
-  }
-  
+  use_fir = (n_iterations < 4) ? true : false;
   if (use_fir) {
-    FloatArray fir_coeffs = agc_coeffs_.agc_spatial_fir_.block(0,stage,3,1);
+    FloatArray fir_coeffs = agc_coeffs_.agc_spatial_fir_.block(0, stage, 3, 1);
     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_.agc_spatial_n_taps_(stage);
-    //Initialize first two taps of stage state, used for both cases
+    // 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.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(0) * ss_tap1) +
@@ -240,29 +454,109 @@
                       (fir_coeffs(2) * ss_tap2);
         break;
       case 5:
-        //Initialize last two taps of stage state, used for 5-tap case
+        // 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.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);
-        ss_tap4.block(0,0,n_ch_-2,1) = stage_state.block(2,0,n_ch_-2,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);
+        ss_tap4.block(0, 0, n_ch_ - 2, 1) =
+          stage_state.block(2, 0, n_ch_ - 2, 1);
         stage_state = (fir_coeffs(0) * (ss_tap3 + ss_tap1)) +
                       (fir_coeffs(1) * stage_state) +
                       (fir_coeffs(2) * (ss_tap2 + ss_tap4));
         break;
       default:
-        //TODO Throw Error
+        // TODO alexbrandmeyer: determine proper error handling implementation.
         std::cout << "Error: bad n-taps in AGCSpatialSmooth" << std::endl;
     }
-    
   } else {
-    stage_state = AGCSmoothDoubleExponential(stage_state);
+    stage_state = AGCSmoothDoubleExponential(stage_state,
+                                             agc_coeffs_.agc_pole_z1_(stage),
+                                             agc_coeffs_.agc_pole_z2_(stage));
   }
   return stage_state;
 }
 
-FloatArray Ear::AGCSmoothDoubleExponential(FloatArray stage_state){
+FloatArray Ear::AGCSmoothDoubleExponential(FloatArray stage_state,
+                                           FPType pole_z1, FPType pole_z2) {
+  int n_pts = stage_state.size();
+  FPType input;
+  FPType state = 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);
+    state = state + (1 - pole_z1) * (input - state);
+  }
+  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);
+    state = state + (1 - pole_z1) * (input - state);
+    stage_state(i) = state;
+  }
   return stage_state;
 }
+
+FloatArray Ear::ReturnZAMemory() {
+  return car_state_.za_memory_;
+}
+
+FloatArray Ear::ReturnZBMemory() {
+  return car_state_.zb_memory_;
+}
+
+FloatArray Ear::ReturnGMemory() {
+  return car_state_.g_memory_;
+};
+
+FloatArray Ear::ReturnZRCoeffs() {
+  return car_coeffs_.zr_coeffs_;
+};
+
+int Ear::ReturnAGCNStages() {
+  return agc_coeffs_.n_agc_stages_;
+}
+
+int Ear::ReturnAGCStateDecimPhase(int stage) {
+  return agc_state_.decim_phase_.at(stage);
+}
+
+FPType Ear::ReturnAGCMixCoeff(int stage) {
+  return agc_coeffs_.agc_mix_coeffs_(stage);
+}
+
+FloatArray Ear::ReturnAGCStateMemory(int stage) {
+  return agc_state_.agc_memory_.at(stage);
+}
+
+int Ear::ReturnAGCDecimation(int stage) {
+  return agc_coeffs_.decimation_.at(stage);
+}
+
+void Ear::SetAGCStateMemory(int stage, FloatArray new_values) {
+  agc_state_.agc_memory_.at(stage) = new_values;
+}
+
+void Ear::SetCARStateDZBMemory(FloatArray new_values) {
+  car_state_.dzb_memory_ = new_values;
+}
+
+void Ear::SetCARStateDGMemory(FloatArray new_values) {
+  car_state_.dg_memory_ = new_values;
+}
+
+FloatArray Ear::StageGValue(FloatArray undamping) {
+  FloatArray r1 = car_coeffs_.r1_coeffs_;
+  FloatArray a0 = car_coeffs_.a0_coeffs_;
+  FloatArray c0 = car_coeffs_.c0_coeffs_;
+  FloatArray h = car_coeffs_.h_coeffs_;
+  FloatArray zr = car_coeffs_.zr_coeffs_;
+  FloatArray r = r1 + zr * undamping;
+  return (1 - 2 * r * a0 + (r * r)) / (1 - 2 * r * a0 + h * r * c0 + (r * r));
+}
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/ear.h
--- a/carfac/ear.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ear.h	Thu May 16 17:33:23 2013 +0000
@@ -27,29 +27,64 @@
 #include "ihc_state.h"
 #include "agc_state.h"
 
-
 class Ear {
-public:
-  int n_ch_;
-  FPType max_channels_per_octave_;
-  CARParams car_params_;
-  IHCParams ihc_params_;
-  AGCParams agc_params_;
+ public:
+  // This is the primary initialization function that is called for each
+  // Ear object in the CARFAC 'Design' method.
+  void InitEar(int n_ch, long fs, FloatArray pole_freqs, CARParams car_p,
+               IHCParams ihc_p, AGCParams agc_p);
+  // These three methods apply the different stages of the model in sequence
+  // to individual audio samples.
+  FloatArray CARStep(FPType input);
+  FloatArray IHCStep(FloatArray car_out);
+  bool AGCStep(FloatArray ihc_out);
+  // These helper functions return portions of the CAR state for storage in the
+  // CAROutput structures.
+  FloatArray ReturnZAMemory();
+  FloatArray ReturnZBMemory();
+  FloatArray ReturnGMemory();
+  FloatArray ReturnZRCoeffs();
+  // These helper functions return portions of the AGC state during the cross
+  // coupling of the ears.
+  int ReturnAGCNStages();
+  int ReturnAGCStateDecimPhase(int stage);
+  FPType ReturnAGCMixCoeff(int stage);
+  int ReturnAGCDecimation(int stage);
+  FloatArray ReturnAGCStateMemory(int stage);
+  // This returns the stage G value during the closing of the AGC loop.
+  FloatArray StageGValue(FloatArray undamping);
+  // This function sets the AGC memory during the cross coupling stage.
+  void SetAGCStateMemory(int stage, FloatArray new_values);
+  // These are two functions to set the CARState dzB and dG memories when
+  // closing the AGC loop
+  void SetCARStateDZBMemory(FloatArray new_values);
+  void SetCARStateDGMemory(FloatArray new_values);
+ private:
+  // These methods carry out the design of the coefficient sets for each of the
+  // three model stages.
+  void DesignFilters(CARParams car_params, long fs, FloatArray pole_freqs);
+  void DesignIHC(IHCParams ihc_params, long fs);
+  void DesignAGC(AGCParams agc_params, long fs);
+  // These are the corresponding methods that initialize the model state
+  // variables before runtime using the model coefficients.
+  void InitIHCState();
+  void InitAGCState();
+  void InitCARState();
+  // These are the various helper functions called during the model runtime.
+  FloatArray OHC_NLF(FloatArray velocities);
+  bool AGCRecurse(int stage, FloatArray agc_in);
+  FloatArray AGCSpatialSmooth(int stage, FloatArray stage_state);
+  FloatArray AGCSmoothDoubleExponential(FloatArray stage_state, FPType pole_z1,
+                                        FPType pole_z2);
+  // These are the private data members that store the state and coefficient
+  // information.
   CARCoeffs car_coeffs_;
   IHCCoeffs ihc_coeffs_;
   AGCCoeffs agc_coeffs_;
   CARState car_state_;
   IHCState ihc_state_;
   AGCState agc_state_;
-  
-  void InitEar(long fs, CARParams car_p, IHCParams ihc_p, AGCParams agc_p);
-  FloatArray CARStep(FPType input);
-  FloatArray OHC_NLF(FloatArray velocities);
-  FloatArray IHCStep(FloatArray car_out);
-  bool AGCStep(FloatArray ihc_out);
-  bool AGCRecurse(int stage, FloatArray agc_in);
-  FloatArray AGCSpatialSmooth(int stage, FloatArray stage_state);
-  FloatArray AGCSmoothDoubleExponential(FloatArray stage_state);
+  int n_ch_;
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/ear_output.cc
--- a/carfac/ear_output.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ear_output.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,16 +22,16 @@
 
 #include "ear_output.h"
 
-void EarOutput::InitOutput(int n_ch, long n_tp){
+void EarOutput::InitOutput(int n_ch, long n_tp) {
   n_ch_ = n_ch;
   n_timepoints_ = n_tp;
-  nap_.resize(n_ch_,n_timepoints_);
-  bm_.resize(n_ch_,n_timepoints_);
-  ohc_.resize(n_ch_,n_timepoints_);
-  agc_.resize(n_ch_,n_timepoints_);
+  nap_.resize(n_ch_, n_timepoints_);
+  bm_.resize(n_ch_, n_timepoints_);
+  ohc_.resize(n_ch_, n_timepoints_);
+  agc_.resize(n_ch_, n_timepoints_);
 }
 
-void EarOutput::MergeOutput(EarOutput ear_output, long start, long length){
+void EarOutput::MergeOutput(EarOutput ear_output, long start, long length) {
   nap_.block(0, start, n_ch_, length) = ear_output.nap_.block(0, 0, n_ch_,
                                                               length);
   bm_.block(0, start, n_ch_, length) = ear_output.bm_.block(0, 0, n_ch_,
diff -r aefe2ca0674f -r 01986636257a carfac/ear_output.h
--- a/carfac/ear_output.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ear_output.h	Thu May 16 17:33:23 2013 +0000
@@ -25,8 +25,9 @@
 
 #include "carfac_common.h"
 
-class EarOutput {
-public:
+struct EarOutput {
+  void InitOutput(int n_ch, long n_tp);
+  void MergeOutput(EarOutput output, long start, long length);
   int n_ch_;
   long n_timepoints_;
   FloatArray2d nap_;
@@ -34,8 +35,6 @@
   FloatArray2d ohc_;
   FloatArray2d agc_;
   FloatArray2d bm_;
-  void InitOutput(int n_ch, long n_tp);
-  void MergeOutput(EarOutput output, long start, long length);
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/ihc_coeffs.cc
--- a/carfac/ihc_coeffs.cc	Mon May 13 22:51:15 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,97 +0,0 @@
-//
-//  ihc_coeffs.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_coeffs.h"
-
-//This method has been created for debugging purposes and depends on <iostream>.
-//Could possibly be removed in the final version to reduce dependencies.
-void IHCCoeffs::OutputCoeffs(){
-  std::cout << "IHCCoeffs Values" << std::endl;
-  std::cout << "****************" << std::endl;
-  std::cout << "n_ch_ = " << n_ch_ << std::endl;
-  std::cout << "just_hwr_ = " << just_hwr_ << std::endl;
-  std::cout << "one_cap_ = " << one_cap_ << std::endl;
-  std::cout << "lpf_coeff_ = " << lpf_coeff_ << std::endl;
-  std::cout << "out1_rate_ = " << out1_rate_ << std::endl;
-  std::cout << "in1_rate_ = " << in1_rate_ << std::endl;
-  std::cout << "out2_rate_ = " << out2_rate_ << std::endl;
-  std::cout << "in2_rate_ = " << in2_rate_ << std::endl;
-  std::cout << "output_gain_ = " << output_gain_ << std::endl;
-  std::cout << "rest_output_ = " << rest_output_ << std::endl;
-  std::cout << "rest_cap1_ = " << rest_cap1_ << std::endl;
-  std::cout << "rest_cap2_ = " << rest_cap2_ << std::endl;
-  std::cout << "ac_coeff_ = " << ac_coeff_ << std::endl << std::endl;
-}
-
-void IHCCoeffs::DesignIHC(IHCParams ihc_params, long fs, int n_ch){
-  if (ihc_params.just_hwr_){
-    n_ch_ = n_ch;
-    just_hwr_ = ihc_params.just_hwr_;
-  } else {
-    if (ihc_params.one_cap_){
-      ro_ = 1 / CARFACDetect(10);
-      c_ = ihc_params.tau1_out_ / ro_;
-      ri_ = ihc_params.tau1_in_ / c_;
-      saturation_output_ = 1 / ((2 * ro_) + ri_);
-      r0_ = 1 / CARFACDetect(0);
-      current_ = 1 / (ri_ + r0_);
-      cap1_voltage_ = 1 - (current_ * ri_);
-      
-      n_ch_ = n_ch;
-      just_hwr_ = false;
-      lpf_coeff_ = 1 - exp( -1 / (ihc_params.tau_lpf_ * fs));
-      out1_rate_ = ro_ / (ihc_params.tau1_out_ * fs);
-      in1_rate_ = 1 / (ihc_params.tau1_in_ * fs);
-      one_cap_ = ihc_params.one_cap_;
-      output_gain_ = 1 / (saturation_output_ - current_);
-      rest_output_ = current_ / (saturation_output_ - current_);
-      rest_cap1_ = cap1_voltage_;
-
-      
-    } else {
-      ro_ = 1 / CARFACDetect(10);
-      c2_ = ihc_params.tau2_out_ / ro_;
-      r2_ = ihc_params.tau2_in_ / c2_;
-      c1_ = ihc_params.tau1_out_ / r2_;
-      r1_ = ihc_params.tau1_in_ / c1_;
-      saturation_output_ = 1 / (2 * ro_ + r2_ + r1_);
-      r0_ = 1 / CARFACDetect(0);
-      current_ = 1 / (r1_ + r2_ + r0_);
-      cap1_voltage_ = 1 - (current_ * r1_);
-      cap2_voltage_ = cap1_voltage_ - (current_ * r2_);
-      
-      n_ch_ = n_ch;
-      just_hwr_ = false;
-      lpf_coeff_ = 1 - exp(-1 / (ihc_params.tau_lpf_ * fs));
-      out1_rate_ = 1 / (ihc_params.tau1_out_ * fs);
-      in1_rate_ = 1 / (ihc_params.tau1_in_ * fs);
-      out2_rate_ = ro_ / (ihc_params.tau2_out_ * fs);
-      in2_rate_ = 1 / (ihc_params.tau2_in_ * fs);
-      one_cap_ = ihc_params.one_cap_;
-      output_gain_ = 1 / (saturation_output_ - current_);
-      rest_output_ = current_ / (saturation_output_ - current_);
-      rest_cap1_ = cap1_voltage_;
-      rest_cap2_ = cap2_voltage_;
-    }
-  }
-  ac_coeff_ = 2 * PI * ihc_params.ac_corner_hz_ / fs;
-}
diff -r aefe2ca0674f -r 01986636257a carfac/ihc_coeffs.h
--- a/carfac/ihc_coeffs.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ihc_coeffs.h	Thu May 16 17:33:23 2013 +0000
@@ -25,9 +25,7 @@
 
 #include "ihc_params.h"
 
-class IHCCoeffs {
-public:
-  int n_ch_;
+struct IHCCoeffs {
   bool just_hwr_;
   bool one_cap_;
   FPType lpf_coeff_;
@@ -40,13 +38,8 @@
   FPType rest_cap1_;
   FPType rest_cap2_;
   FPType ac_coeff_;
-  
-  FPType ro_, c_, ri_, r0_, saturation_output_, current_, cap1_voltage_,
-         cap2_voltage_;
-  FPType c2_, r2_, c1_, r1_;
-  
-  void OutputCoeffs(); //Method to view coeffs, could go in final version
-  void DesignIHC(IHCParams ihc_params, long fs, int n_ch);
+  FPType cap1_voltage_;
+  FPType cap2_voltage_;
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/ihc_params.cc
--- a/carfac/ihc_params.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ihc_params.cc	Thu May 16 17:33:23 2013 +0000
@@ -22,38 +22,24 @@
 
 #include "ihc_params.h"
 
-
-//Default constructor for IHCParams initializes with the settings from Lyon's
-//book 'Human and Machine Hearing'
-IHCParams::IHCParams(){
-  just_hwr_ = false; //not just a simple HWR
-  one_cap_ = true; //uses the Allen model, as in Lyon's book
-  tau_lpf_ = 0.000080; //80 microseconds smoothing twice
-  tau1_out_ = 0.0005; //depletion tau is pretty fast
-  tau1_in_ = 0.010; //recovery tau is slower
+// The default constructor for IHCParams initializes with the settings from
+// Lyon's book 'Human and Machine Hearing'
+IHCParams::IHCParams() {
+  just_hwr_ = false;
+  one_cap_ = 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;
 }
 
-//OutputParams method uses <iostream> for debugging purposes, could go in the
-//final version
-void IHCParams::OutputParams(){
-  std::cout << "IHCParams Values" << std::endl;
-  std::cout << "****************" << std::endl;
-  std::cout << "just_hwr_ = " << just_hwr_ << std::endl;
-  std::cout << "one_cap_ = " << one_cap_ << std::endl;
-  std::cout << "tau_lpf_ = " << tau_lpf_ << std::endl;
-  std::cout << "tau1_out_ = " << tau1_out_ << std::endl;
-  std::cout << "tau1_in_ = " << tau1_in_ << std::endl;
-  std::cout << "tau2_out_ = " << tau2_out_ << std::endl;
-  std::cout << "tau2_in_ = " << tau2_in_ << std::endl;
-  std::cout << "ac_corner_hz_ = " << ac_corner_hz_ << std::endl << std::endl;
-}
-
-//SetParams method allows for use of different inner hair cell parameters
-void IHCParams::SetParams(bool jh, bool oc, FPType tlpf, FPType t1out,
-                          FPType t1in, FPType t2out, FPType t2in, FPType acchz){
+// The overloaded constructor allows for use of different inner hair cell
+// parameters.
+IHCParams::IHCParams(bool jh, bool oc, FPType tlpf, FPType t1out,
+                          FPType t1in, FPType t2out, FPType t2in,
+                          FPType acchz) {
   just_hwr_ = jh;
   one_cap_ = oc;
   tau_lpf_ = tlpf;
diff -r aefe2ca0674f -r 01986636257a carfac/ihc_params.h
--- a/carfac/ihc_params.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ihc_params.h	Thu May 16 17:33:23 2013 +0000
@@ -25,8 +25,10 @@
 
 #include "carfac_common.h"
 
-class IHCParams {
-public:
+struct IHCParams {
+  IHCParams();
+  IHCParams(bool jh, bool oc, FPType tlpf, FPType t1out, FPType t1in,
+                 FPType t2out, FPType t2in, FPType acchz);
   bool just_hwr_;
   bool one_cap_;
   FPType tau_lpf_;
@@ -35,11 +37,6 @@
   FPType tau2_out_;
   FPType tau2_in_;
   FPType ac_corner_hz_;
-  
-  IHCParams();
-  void OutputParams(); //Output current parameter values using std::cout
-  void SetParams(bool jh, bool oc, FPType tlpf, FPType t1out, FPType t1in,
-                 FPType t2out, FPType t2in, FPType acchz);//Method to set non-default params
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/ihc_state.cc
--- a/carfac/ihc_state.cc	Mon May 13 22:51:15 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,39 +0,0 @@
-//
-//  ihc_state.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_state.h"
-
-void IHCState::InitIHCState(IHCCoeffs ihc_coeffs){
-  n_ch_ = ihc_coeffs.n_ch_;
-  ihc_accum_ = FloatArray::Zero(n_ch_);
-  if (! ihc_coeffs.just_hwr_){
-    ac_coupler_ = FloatArray::Zero(n_ch_);
-    lpf1_state_ = ihc_coeffs.rest_output_ * FloatArray::Ones(n_ch_);
-    lpf2_state_ = ihc_coeffs.rest_output_ * FloatArray::Ones(n_ch_);
-    if (ihc_coeffs.one_cap_){
-      cap1_voltage_ = ihc_coeffs.rest_cap1_ * FloatArray::Ones(n_ch_);
-    } else {
-      cap1_voltage_ = ihc_coeffs.rest_cap1_ * FloatArray::Ones(n_ch_);
-      cap2_voltage_ = ihc_coeffs.rest_cap2_ * FloatArray::Ones(n_ch_);
-    }
-  }
-}
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/ihc_state.h
--- a/carfac/ihc_state.h	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/ihc_state.h	Thu May 16 17:33:23 2013 +0000
@@ -25,18 +25,13 @@
 
 #include "ihc_coeffs.h"
 
-class IHCState {
-public:
-  int n_ch_;
+struct IHCState {
   FloatArray ihc_accum_;
   FloatArray cap1_voltage_;
   FloatArray cap2_voltage_;
   FloatArray lpf1_state_;
   FloatArray lpf2_state_;
   FloatArray ac_coupler_;
-  
-  void InitIHCState(IHCCoeffs ihc_coeffs);
-  
 };
 
-#endif
+#endif
\ No newline at end of file
diff -r aefe2ca0674f -r 01986636257a carfac/main.cc
--- a/carfac/main.cc	Mon May 13 22:51:15 2013 +0000
+++ b/carfac/main.cc	Thu May 16 17:33:23 2013 +0000
@@ -23,7 +23,7 @@
 // *****************************************************************************
 // main.cc
 // *****************************************************************************
-// This 'main' file is not currently intended as part of the CARFAC distribution,
+// This 'main' file is not currently intended as part of the CARFAC distribution
 // but serves as a testbed for debugging and implementing various aspects of the
 // library. It currently includes the "libsndfile" API for loading soundfiles:
 //
@@ -35,36 +35,54 @@
 // design stage) and sound data (for running the model).
 #include <sndfile.h>
 #include "carfac.h"
+#include <fstream>
+//GoogleTest is now included for running unit tests
+#include <gtest/gtest.h>
 
 // ReadSound takes a character array (filename string) as an argument and
 // returns a two dimensional (samples x channels) FloatArray (Eigen ArrayXX)
 // containing the sound data
-FloatArray2d ReadSound(const char * filename){
-  FloatArray2d mysnd; //output data
+FloatArray2d ReadSound(const char * filename) {
+  // This output variable stores the sound data in an Eigen Array:
+  FloatArray2d mysnd;
+  // These two structures from the libsndfile library are used to store sound
+  // file information and load the data.
   SNDFILE *sf; 
   SF_INFO info;
+  // Several scalars are used to store relevant information needed for
+  // transfering data from the WAV file to an Eigen Array. 
   long num, num_items;
   double *buf;
-  long f,sr,c;
+  long f, sr, c;
+  // This opens the sound file and prints an error message if the file can't
+  // be found.
   sf = sf_open(filename,SFM_READ,&info);
   if (sf == NULL)
   {
     std::cout << "Failed to open the file" << std::endl;
     return mysnd;
   }
+  // Here we store relevant header information in our scalars and use them to
+  // read into our data buffer.
   f = info.frames;
   sr = info.samplerate;
   c = info.channels;
-  num_items = f*c;
+  num_items = f * c;
   buf = new double[num_items];
-  num = sf_read_double(sf,buf,num_items);
-  mysnd.resize(f,c);
+  num = sf_read_double(sf, buf, num_items);
+  // We resize our Eigen array to hold the correct number of samples and
+  // audio channels.
+  mysnd.resize(f, c);
+  // Now we move through the data buffer and store the stereo audio channels
+  // in the two columns of the Eigen Array
+  // TODO alexbrandmeyer: this assumes stereo and will give problems otherwise.
   int j = 0;
-  for(int i = 0; i < num_items; i = i + 2){
-    mysnd(j,0) = buf[i];
-    mysnd(j,1) = buf[i+1];
+  for(int i = 0; i < num_items; i = i + 2) {
+    mysnd(j, 0) = buf[i];
+    mysnd(j, 1) = buf[i + 1];
     j++;
   }
+  // Now it's time to close the audio file and return the Eigen Array as output.
   sf_close(sf);
   return mysnd;
 };
@@ -72,46 +90,46 @@
 // ReadSoundInfo takes a character array (filename string) as an argument and
 // returns an SF_INFO structure containing the sample rate and channel info
 // needed during the call to CARFAC::Design
-SF_INFO ReadSoundInfo(const char * filename){
+SF_INFO ReadSoundInfo(const char * filename) {
   SNDFILE *sf;
   SF_INFO info;
   sf = sf_open(filename,SFM_READ,&info);
-  if (sf == NULL)
-  {
+  if (sf == NULL) {
     std::cout << "Failed to open the file" << std::endl;
     return info;
   }
   return info;
 };
 
-
-
-
 // This 'main' function serves as the primary testbed for this C++ CARFAC
 // implementation. It currently uses a hardcoded filename to obtain sound file
 // info and sound data, and designs a CARFAC on the basis of the header data
 // using the default parameters.
-int main()
-{
-  //Here we specify a path to a test file
+int main(int argc, char **argv) {
+  // This initializes the GoogleTest unit testing framework.
+  ::testing::InitGoogleTest(&argc, argv);
+  // Here we specify a path to a test file.
   const char * filename = "/Users/alexbrandmeyer/aimc/carfac/test_signal.wav";
-  
-  //Now we load the header info and sound data
+  // This loads the header info and sound data.
   SF_INFO info = ReadSoundInfo(filename);
   FloatArray2d mysnd = ReadSound(filename);
-  //These initialze the default parameter objects needed for the CARFAC design
+  // These initialze the default parameter objects needed for the CARFAC design.
   CARParams *car_params = new CARParams();
   IHCParams *ihc_params = new IHCParams();
   AGCParams *agc_params = new AGCParams();
-  
-  //This initializes the CARFAC object and runs the design method
+  // This initializes the CARFAC object and runs the design method.
   CARFAC *mycf = new CARFAC();
-  mycf->Design(info.channels,info.samplerate, *car_params, *ihc_params,
+  mycf->Design(info.channels, info.samplerate, *car_params, *ihc_params,
                *agc_params);
-  std::cout << "CARFAC Object Created" << std::endl;
-  
-  //Now we run the model on the test data
-  CARFACOutput output = mycf->Run(mysnd);
+  // This clears the parameters which are no longer needed.
+  delete car_params;
+  delete ihc_params;
+  delete agc_params;
+  // Now we run the model on the test data (using a closed loop for now).
+  CARFACOutput output = mycf->Run(mysnd, false);
+  // Finally we clear the CARFAC object when we're done.
+  delete mycf;
+  //return RUN_ALL_TESTS();;
   return 0;
 }