Mercurial > hg > aimc

changeset 504:a0869cb1c99b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Major update to how the DOHC works; like in recent book OHC chapter; Design Doc update (a bit)
author dicklyon@google.com
date Thu, 24 May 2012 22:26:56 +0000
parents 37c007925536
children db0e5e86fddd
files matlab/bmm/carfac/CARFAC_AGC_Step.m matlab/bmm/carfac/CARFAC_CAR_Step.m matlab/bmm/carfac/CARFAC_Close_AGC_Loop.m matlab/bmm/carfac/CARFAC_Design.m matlab/bmm/carfac/CARFAC_Design_Doc.txt matlab/bmm/carfac/CARFAC_IHC_Step.m matlab/bmm/carfac/CARFAC_Init.m matlab/bmm/carfac/CARFAC_OHC_NLF.m matlab/bmm/carfac/CARFAC_Run.m matlab/bmm/carfac/CARFAC_Run_Linear.m matlab/bmm/carfac/CARFAC_Run_Segment.m matlab/bmm/carfac/CARFAC_Stage_g.m matlab/bmm/carfac/CARFAC_Transfer_Functions.m matlab/bmm/carfac/CARFAC_hacking.m matlab/bmm/carfac/MultiScaleSmooth.m
diffstat 15 files changed, 188 insertions(+), 199 deletions(-) [+]
line wrap: on
line diff
--- a/matlab/bmm/carfac/CARFAC_AGC_Step.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_AGC_Step.m	Thu May 24 22:26:56 2012 +0000
@@ -17,8 +17,8 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.
 
-function [state, updated] = CARFAC_AGC_Step(coeffs, detects, state)
-% function [state, updated] = CARFAC_AGC_Step(coeffs, detects, state)
+function [state, updated] = CARFAC_AGC_Step(detects, coeffs, state)
+% function [state, updated] = CARFAC_AGC_Step(detects, coeffs, state)
 %
 % one time step of the AGC state update; decimates internally
 
--- a/matlab/bmm/carfac/CARFAC_CAR_Step.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_CAR_Step.m	Thu May 24 22:26:56 2012 +0000
@@ -17,21 +17,27 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.
 
-function [zY, state] = CARFAC_CAR_Step(x_in, CAR_coeffs, state)
+function [car_out, state] = CARFAC_CAR_Step(x_in, CAR_coeffs, state)
 % function [zY, state] = CARFAC_CAR_Step(x_in, CAR_coeffs, state)
 %
 % One sample-time update step for the filter part of the CARFAC.
 
 % Most of the update is parallel; finally we ripple inputs at the end.
 
-% Local nonlinearity zA and AGC feedback zB reduce pole radius:
+
+% do the DOHC stuff:
+
+g = state.g_memory + state.dg_memory;  % interp g
+zB = state.zB_memory + state.dzB_memory; % AGC interpolation state
+% update the nonlinear function of "velocity", and zA (delay of z2):
 zA = state.zA_memory;
-zB = state.zB_memory + state.dzB_memory; % AGC interpolation
-r1 = CAR_coeffs.r1_coeffs;
-g = state.g_memory + state.dg_memory;  % interp g
+v = state.z2_memory - zA;
+% nlf = CARFAC_OHC_NLF(v .* widen, CAR_coeffs);  % widen v with feedback
+nlf = CARFAC_OHC_NLF(v, CAR_coeffs);
+% zB * nfl is "undamping" delta r:
+r = CAR_coeffs.r1_coeffs + zB .* nlf; 
+zA = state.z2_memory;
 
-% zB and zA are "extra damping", and multiply zr (compressed theta):
-r = r1 - CAR_coeffs.zr_coeffs .* (zA + zB); 
 
 % now reduce state by r and rotate with the fixed cos/sin coeffs:
 z1 = r .* (CAR_coeffs.a0_coeffs .* state.z1_memory - ...
@@ -40,9 +46,6 @@
 z2 = r .* (CAR_coeffs.c0_coeffs .* state.z1_memory + ...
   CAR_coeffs.a0_coeffs .* state.z2_memory);
 
-% update the nonlinear function of "velocity", into zA:
-zA = CARFAC_OHC_NLF(state.z2_memory - z2, CAR_coeffs);
-
 zY = CAR_coeffs.h_coeffs .* z2;  % partial output
 
 % Ripple input-output path, instead of parallel, to avoid delay...
@@ -57,7 +60,7 @@
 end
 
 % put new state back in place of old
-% (z1 and z2 are genuine temps; the others can update by reference in C)
+% (z1 is a genuine temp; the others can update by reference in C)
 state.z1_memory = z1;
 state.z2_memory = z2;
 state.zA_memory = zA;
@@ -65,3 +68,6 @@
 state.zY_memory = zY;
 state.g_memory = g;
 
+car_out = zY;
+
+
--- a/matlab/bmm/carfac/CARFAC_Close_AGC_Loop.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Close_AGC_Loop.m	Thu May 24 22:26:56 2012 +0000
@@ -24,12 +24,13 @@
 decim1 = CF.AGC_params.decimation(1);
 
 for ear = 1:CF.n_ears
-  extra_damping = CF.ears(ear).AGC_state.AGC_memory(:, 1);  % stage 1 result
+  undamping = 1 - CF.ears(ear).AGC_state.AGC_memory(:, 1); % stage 1 result
   % Update the target stage gain for the new damping:
-  new_g = CARFAC_Stage_g(CF.ears(ear).CAR_coeffs, extra_damping);
+  new_g = CARFAC_Stage_g(CF.ears(ear).CAR_coeffs, undamping);
   % set the deltas needed to get to the new damping:
   CF.ears(ear).CAR_state.dzB_memory = ...
-    (extra_damping - CF.ears(ear).CAR_state.zB_memory) / decim1;
+    (CF.ears(ear).CAR_coeffs.zr_coeffs .* undamping - ...
+    CF.ears(ear).CAR_state.zB_memory) / decim1;
   CF.ears(ear).CAR_state.dg_memory = ...
     (new_g - CF.ears(ear).CAR_state.g_memory) / decim1;
 end
--- a/matlab/bmm/carfac/CARFAC_Design.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Design.m	Thu May 24 22:26:56 2012 +0000
@@ -52,11 +52,11 @@
 
 if nargin < 3
   CF_CAR_params = struct( ...
-    'velocity_scale', 0.2, ...  % for the "cubic" velocity nonlinearity
+    'velocity_scale', 0.05, ...  % for the velocity nonlinearity
     'v_offset', 0.04, ...  % offset gives a quadratic part
     'v2_corner', 0.2, ...  % corner for essential nonlin
-    'v_damp_max', 0.01, ... % damping delta damping from velocity nonlin
     'min_zeta', 0.10, ... % minimum damping factor in mid-freq channels
+    'max_zeta', 0.35, ... % maximum damping factor in mid-freq channels
     'first_pole_theta', 0.85*pi, ...
     'zero_ratio', sqrt(2), ... % how far zero is above pole
     'high_f_damping_compression', 0.5, ... % 0 to 1 to compress zeta
@@ -71,10 +71,9 @@
     'n_stages', 4, ...
     'time_constants', [1, 4, 16, 64]*0.002, ...
     'AGC_stage_gain', 2, ...  % gain from each stage to next slower stage
-    'decimation', [4, 2, 2, 2], ...  % how often to update the AGC states
+    'decimation', [8, 2, 2, 2], ...  % how often to update the AGC states
     'AGC1_scales', [1.0, 1.4,  2.0, 2.8], ...   % in units of channels
     'AGC2_scales', [1.6, 2.25, 3.2, 4.5], ... % spread more toward base
-    'detect_scale', 0.25, ...  % the desired damping range
     'AGC_mix_coeff', 0.5);
 end
 
@@ -83,7 +82,8 @@
   one_cap = 0;         % bool; 0 for new two-cap hack
   just_hwr = 0;        % book; 0 for normal/fancy IHC; 1 for HWR
   if just_hwr
-    CF_IHC_params = struct('just_hwr', 1);  % just a simple HWR
+    CF_IHC_params = struct('just_hwr', 1, ...  % just a simple HWR
+        'ac_corner_Hz', 20);
   else
     if one_cap
       CF_IHC_params = struct( ...
@@ -91,7 +91,8 @@
         'one_cap', one_cap, ...   % bool; 0 for new two-cap hack
         'tau_lpf', 0.000080, ...  % 80 microseconds smoothing twice
         'tau_out', 0.0005, ...    % depletion tau is pretty fast
-        'tau_in', 0.010 );        % recovery tau is slower
+        'tau_in', 0.010, ...        % recovery tau is slower
+        'ac_corner_Hz', 20);
     else
       CF_IHC_params = struct( ...
         'just_hwr', just_hwr, ...        % not just a simple HWR
@@ -100,7 +101,8 @@
         'tau1_out', 0.010, ...    % depletion tau is pretty fast
         'tau1_in', 0.020, ...     % recovery tau is slower
         'tau2_out', 0.0025, ...   % depletion tau is pretty fast
-        'tau2_in', 0.005 );        % recovery tau is slower
+        'tau2_in', 0.005, ...        % recovery tau is slower
+        'ac_corner_Hz', 20);
     end
   end
 end
@@ -164,8 +166,7 @@
   'n_ch', n_ch, ...
   'velocity_scale', CAR_params.velocity_scale, ...
   'v_offset', CAR_params.v_offset, ...
-  'v2_corner', CAR_params.v2_corner, ...
-  'v_damp_max', CAR_params.v_damp_max ...
+  'v2_corner', CAR_params.v2_corner ...
   );
 
 % don't really need these zero arrays, but it's a clue to what fields
@@ -194,18 +195,19 @@
 % Compress theta to give somewhat higher Q at highest thetas:
 ff = CAR_params.high_f_damping_compression;  % 0 to 1; typ. 0.5
 x = theta/pi;
+
 zr_coeffs = pi * (x - ff * x.^3);  % when ff is 0, this is just theta,
 %                       and when ff is 1 it goes to zero at theta = pi.
-CAR_coeffs.zr_coeffs = zr_coeffs;  % how r relates to zeta
+max_zeta = CAR_params.max_zeta;
+CAR_coeffs.r1_coeffs = (1 - zr_coeffs .* max_zeta);  % "r1" for the max-damping condition
 
 min_zeta = CAR_params.min_zeta;
-% increase the min damping where channels are spaced out more:
-
-min_zeta = min_zeta + 0.25*(ERB_Hz(pole_freqs, ...
+% Increase the min damping where channels are spaced out more, by pulling 
+% 25% of the way toward ERB_Hz/pole_freqs (close to 0.1 at high f)
+min_zetas = min_zeta + 0.25*(ERB_Hz(pole_freqs, ...
   CAR_params.ERB_break_freq, CAR_params.ERB_Q) ./ pole_freqs - min_zeta);
-r1 = (1 - zr_coeffs .* min_zeta);  % "1" for the min-damping condition
-
-CAR_coeffs.r1_coeffs = r1;
+CAR_coeffs.zr_coeffs = zr_coeffs .* ...
+  (max_zeta - min_zetas);  % how r relates to undamping
 
 % undamped coupled-form coefficients:
 CAR_coeffs.a0_coeffs = a0;
@@ -216,10 +218,10 @@
 CAR_coeffs.h_coeffs = h;
 
 % for unity gain at min damping, radius r; only used in CARFAC_Init:
-extra_damping = zeros(size(r1));
+relative_undamping = ones(n_ch, 1);  % max undamping to start
 % this function needs to take CAR_coeffs even if we haven't finished
 % constucting it by putting in the g0_coeffs:
-CAR_coeffs.g0_coeffs = CARFAC_Stage_g(CAR_coeffs, extra_damping);
+CAR_coeffs.g0_coeffs = CARFAC_Stage_g(CAR_coeffs, relative_undamping);
 
 
 %% the AGC design coeffs:
@@ -310,14 +312,8 @@
 
 AGC_coeffs.AGC_gain = total_DC_gain;
 
-% adjust the detect_scale by the total DC gain of the AGC filters:
-AGC_coeffs.detect_scale = AGC_params.detect_scale / total_DC_gain;
-
-% % print some results
-AGC_coeffs
-AGC_spatial_FIR = AGC_coeffs.AGC_spatial_FIR
-AGC_spatial_iterations = AGC_coeffs.AGC_spatial_iterations
-AGC_spatial_n_taps = AGC_coeffs.AGC_spatial_n_taps
+% adjust the detect_scale to be the reciprocal DC gain of the AGC filters:
+AGC_coeffs.detect_scale = 1 / total_DC_gain;
 
 
 %%
@@ -355,7 +351,7 @@
     'just_hwr', 1);
 else
   if IHC_params.one_cap
-    ro = 1 / CARFAC_Detect(2);  % output resistance
+    ro = 1 / CARFAC_Detect(10);  % output resistance at a very high level
     c = IHC_params.tau_out / ro;
     ri = IHC_params.tau_in / c;
     % to get steady-state average, double ro for 50% duty cycle
@@ -381,7 +377,7 @@
       'lpf2_state', 0, ...
       'ihc_accum', 0);
   else
-    ro = 1 / CARFAC_Detect(2);  % output resistance
+    ro = 1 / CARFAC_Detect(10);  % output resistance at a very high level
     c2 = IHC_params.tau2_out / ro;
     r2 = IHC_params.tau2_in / c2;
     c1 = IHC_params.tau1_out / r2;
@@ -415,6 +411,8 @@
       'ihc_accum', 0);
   end
 end
+% one more late addition that applies to all cases:
+IHC_coeffs.ac_coeff = 2 * pi * IHC_params.ac_corner_Hz / fs;
 
 %%
 % default design result, running this function with no args, should look
@@ -422,12 +420,15 @@
 %
 %
 % CF = CARFAC_Design
-% CF.CAR_params
-% CF.AGC_params
-% CF.CAR_coeffs
-% CF.AGC_coeffs
-% CF.IHC_coeffs
-% CF =
+% CAR_params = CF.CAR_params
+% AGC_params = CF.AGC_params
+% IHC_params = CF.IHC_params
+% CAR_coeffs = CF.ears(1).CAR_coeffs
+% AGC_coeffs = CF.ears(1).AGC_coeffs
+% AGC_spatial_FIR = AGC_coeffs.AGC_spatial_FIR
+% IHC_coeffs = CF.ears(1).IHC_coeffs
+
+% CF = 
 %                          fs: 22050
 %     max_channels_per_octave: 12.2709
 %                  CAR_params: [1x1 struct]
@@ -435,16 +436,14 @@
 %                  IHC_params: [1x1 struct]
 %                        n_ch: 71
 %                  pole_freqs: [71x1 double]
-%                  CAR_coeffs: [1x1 struct]
-%                  AGC_coeffs: [1x1 struct]
-%                  IHC_coeffs: [1x1 struct]
-%                      n_ears: 0
-% ans =
-%                 velocity_scale: 0.2000
-%                       v_offset: 0.0100
+%                        ears: [1x1 struct]
+%                      n_ears: 1
+% CAR_params = 
+%                 velocity_scale: 0.0500
+%                       v_offset: 0.0400
 %                      v2_corner: 0.2000
-%                     v_damp_max: 0.0100
 %                       min_zeta: 0.1000
+%                       max_zeta: 0.3500
 %               first_pole_theta: 2.6704
 %                     zero_ratio: 1.4142
 %     high_f_damping_compression: 0.5000
@@ -452,28 +451,35 @@
 %                    min_pole_Hz: 30
 %                 ERB_break_freq: 165.3000
 %                          ERB_Q: 9.2645
-% ans =
+% AGC_params = 
 %           n_stages: 4
 %     time_constants: [0.0020 0.0080 0.0320 0.1280]
 %     AGC_stage_gain: 2
 %         decimation: [8 2 2 2]
 %        AGC1_scales: [1 1.4000 2 2.8000]
 %        AGC2_scales: [1.6000 2.2500 3.2000 4.5000]
-%       detect_scale: 0.2500
 %      AGC_mix_coeff: 0.5000
-% ans =
+% IHC_params = 
+%         just_hwr: 0
+%          one_cap: 0
+%          tau_lpf: 8.0000e-05
+%         tau1_out: 0.0100
+%          tau1_in: 0.0200
+%         tau2_out: 0.0025
+%          tau2_in: 0.0050
+%     ac_corner_Hz: 20
+% CAR_coeffs = 
 %               n_ch: 71
-%     velocity_scale: 0.2000
-%           v_offset: 0.0100
+%     velocity_scale: 0.0500
+%           v_offset: 0.0400
 %          v2_corner: 0.2000
-%         v_damp_max: 0.0100
 %          r1_coeffs: [71x1 double]
 %          a0_coeffs: [71x1 double]
 %          c0_coeffs: [71x1 double]
 %           h_coeffs: [71x1 double]
 %          g0_coeffs: [71x1 double]
 %          zr_coeffs: [71x1 double]
-% ans =
+% AGC_coeffs = 
 %                       n_ch: 71
 %               n_AGC_stages: 4
 %             AGC_stage_gain: 2
@@ -486,17 +492,25 @@
 %         AGC_spatial_n_taps: [3 3 3 3]
 %             AGC_mix_coeffs: [0 0.0454 0.0227 0.0113]
 %                   AGC_gain: 15
-%               detect_scale: 0.0167
-% ans =
+%               detect_scale: 0.0667
+% AGC_spatial_FIR =
+%     0.2744    0.2829    0.2972    0.2999
+%     0.3423    0.3571    0.3512    0.3616
+%     0.3832    0.3600    0.3516    0.3385
+% IHC_coeffs = 
 %            n_ch: 71
 %        just_hwr: 0
 %       lpf_coeff: 0.4327
 %       out1_rate: 0.0045
 %        in1_rate: 0.0023
-%       out2_rate: 0.0267
+%       out2_rate: 0.0199
 %        in2_rate: 0.0091
 %         one_cap: 0
-%     output_gain: 17.9162
-%     rest_output: 0.5240
-%       rest_cap2: 0.7421
-%       rest_cap1: 0.8281
+%     output_gain: 12.1185
+%     rest_output: 0.3791
+%       rest_cap2: 0.7938
+%       rest_cap1: 0.8625
+%        ac_coeff: 0.0057
+
+
+
--- a/matlab/bmm/carfac/CARFAC_Design_Doc.txt	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Design_Doc.txt	Thu May 24 22:26:56 2012 +0000
@@ -1,6 +1,6 @@
 CARFAC Design Doc
 by "Richard F. Lyon" <dicklyon@google.com>
-
+updated 24 May 2012
 
 The CAR-FAC (cascade of asymmetric resonators with fast-acting
 compression) is a cochlear model implemented as an efficient sound
@@ -25,11 +25,10 @@
 
 The three main subsystems of the EAR are the cascade of asymmetric
 resonators (CAR), the inner hair cell (IHC), and the automatic gain
-control (AGC).  These are not intended to work independently, so may
-not have corresponding classes to combine all their aspects (in the
-Matlab so far, they don't).  Rather, each part has three classes
-associated with it, and the CARFAC stores instances of all of these as
-members:
+control (AGC).  These are not intended to work independently, but
+each part has three groups of data associated with it.  The CARFAC 
+stores instances of the "params" that drive the design, and the
+"coeffs" and "state" are stored per ear:
 
   CAR_params, CAR_coeffs, CAR_state
   IHC_params, IHC_coeffs, IHC_state
@@ -44,10 +43,9 @@
 model on samples or segments of sound waveforms, such as the state
 variables of the digital filters.
 
-At construction ("design") time, params are provided, and coeffs are
-created.  The created CARFAC stores the params that it was designed
-from, and the coeffs that will be used at run time; there is no state
-yet:
+At construction ("design") time, params are provided, and ears with coeffs 
+are created.  The created CARFAC stores the params that it was designed
+from; there is no state yet:
 
   CF = CARFAC_Design(CAR_params, IHC_params, AGC_params)
 
@@ -63,19 +61,22 @@
 not earlier at design time; it is not in params since it doesn't
 affect coefficient design).
 
-Multi-ear designs always share the same coeffs, but need separate state.
-The only  place the ears interact is in the AGC.  The AGC_state is
-one object that holds the state of all the ears, so they can be
-concurrently updated.  For the CAR and IHC, though, the state object
-represents a single ear's state, and the CARFAC stores arrays (vectors) 
-of state objects for the different ears.  The functions that update the
-CAR and IHC states are then simple single-ear functions.
+Multi-ear designs usually have the same coeffs across several ears, but 
+always need separate state.  
+The coeffs are kept separate per ear in case someone wants to
+modify one or both to simulate asymmetric hearing loss or such.
+
+The only  place the ears interact is in the AGC.  The function that
+closes the AGC loop can "mix" their states, making an inter-ear
+coupling.  Other than that, the functions that update the CAR and IHC 
+and AGC states are simple single-ear functions (methods of the ear
+class?).
 
 
 Data size and performance:
 
 The coeffs and states are several kilobytes each, since they store a
-handful (5 to 10) of floating-point values (at 4 or 8 bytes each) for
+handful (10 or so) of floating-point values (at 4 or 8 bytes each) for
 each channel (typically 60 to 100 channels); that 240 to 800 bytes per
 coefficient and per state variable.  That's the entire data memory
 footprint; most of it is accessed at every sample time (hopefully it
@@ -102,15 +103,11 @@
 length of 1 sample if results are needed with very low latency.
 
 Internally, CARFAC_Run updates each part of the state, one sample at a
-time.  First the CAR and IHC are updated ("stepped") for all ears:
+time.  First the CAR, IHC, and AGC are updated ("stepped") for all ears:
 
   [car_out, CAR_state] = CARFAC_CAR_Step(sample,  CAR_coeffs, CAR_state)
   [ihc_out, IHC_state] = CARFAC_IHC_Step(car_out, IHC_coeffs, IHC_state)
-
-Then the AGC is updated, for all ears at once (note the plural
-"ihc_outs", representing the collection of n_ears):
-
-  [AGC_state, updated] = CARFAC_AGC_Step(AGC_coeffs, ihc_outs, AGC_state)
+  [AGC_state, updated] = CARFAC_AGC_Step(ihc_out, AGC_coeffs, AGC_state)
 
 The AGC filter mostly runs at a lower sample rate (by a factor of 8 by 
 default). Usually it just accumulates its input and returns quickly. The 
--- a/matlab/bmm/carfac/CARFAC_IHC_Step.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_IHC_Step.m	Thu May 24 22:26:56 2012 +0000
@@ -23,14 +23,15 @@
 % One sample-time update of inner-hair-cell (IHC) model, including the
 % detection nonlinearity and one or two capacitor state variables.
 
-just_hwr = coeffs.just_hwr;
+% AC couple the filters_out, with 20 Hz corner
+ac_diff = filters_out - state.ac_coupler;
+state.ac_coupler = state.ac_coupler + coeffs.ac_coeff * ac_diff;
 
-if just_hwr
-  ihc_out = min(2, max(0, filters_out));  % limit it for stability
-  state.ihc_accum = state.ihc_accum + ihc_out;
+if coeffs.just_hwr
+  ihc_out = min(2, max(0, ac_diff));  % limit it for stability
 else
-  conductance = CARFAC_Detect(filters_out);  % detect with HWR or so
-  
+  conductance = CARFAC_Detect(ac_diff);  % rectifying nonlinearity
+
   if coeffs.one_cap;
     ihc_out = conductance .* state.cap_voltage;
     state.cap_voltage = state.cap_voltage - ihc_out .* coeffs.out_rate + ...
@@ -54,5 +55,7 @@
   state.lpf2_state = state.lpf2_state + coeffs.lpf_coeff * ...
     (state.lpf1_state - state.lpf2_state);
   ihc_out = state.lpf2_state - coeffs.rest_output;
-  state.ihc_accum = state.ihc_accum + ihc_out;  % for where decimated output is useful
 end
+
+state.ihc_accum = state.ihc_accum + ihc_out;  % for where decimated output is useful
+
--- a/matlab/bmm/carfac/CARFAC_Init.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Init.m	Thu May 24 22:26:56 2012 +0000
@@ -67,7 +67,8 @@
       'ihc_accum', zeros(n_ch, 1), ...
       'cap_voltage', coeffs.rest_cap * ones(n_ch, 1), ...
       'lpf1_state', coeffs.rest_output * ones(n_ch, 1), ...
-      'lpf2_state', coeffs.rest_output * ones(n_ch, 1) ...
+      'lpf2_state', coeffs.rest_output * ones(n_ch, 1), ...
+      'ac_coupler', zeros(n_ch, 1) ...
       );
   else
     state = struct( ...
@@ -75,7 +76,8 @@
       'cap1_voltage', coeffs.rest_cap1 * ones(n_ch, 1), ...
       'cap2_voltage', coeffs.rest_cap2* ones(n_ch, 1), ...
       'lpf1_state', coeffs.rest_output * ones(n_ch, 1), ...
-      'lpf2_state', coeffs.rest_output * ones(n_ch, 1) ...
+      'lpf2_state', coeffs.rest_output * ones(n_ch, 1), ...
+      'ac_coupler', zeros(n_ch, 1) ...
       );
   end
 end
--- a/matlab/bmm/carfac/CARFAC_OHC_NLF.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_OHC_NLF.m	Thu May 24 22:26:56 2012 +0000
@@ -1,7 +1,28 @@
+% Copyright 2012, Google, Inc.
+% Author: Richard F. Lyon
+%
+% This Matlab 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.
+
 function nlf = CARFAC_OHC_NLF(velocities, CAR_coeffs)
+% function nlf = CARFAC_OHC_NLF(velocities, CAR_coeffs)
+% start with a quadratic nonlinear function, and limit it via a
+% rational function; make the result go to zero a high
+% absolute velocities, so it will do nothing there.
 
-nlf = ((velocities .* CAR_coeffs.velocity_scale) + ...
+qnlf = ((velocities .* CAR_coeffs.velocity_scale) + ...
   CAR_coeffs.v_offset) .^ 2;
-% soft saturation to make it more like an "essential" nonlinearity:
-nlf = CAR_coeffs.v_damp_max * nlf ./ (CAR_coeffs.v2_corner + nlf);
-
+nlf = 1 - qnlf ./ (CAR_coeffs.v2_corner + qnlf);
--- a/matlab/bmm/carfac/CARFAC_Run.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Run.m	Thu May 24 22:26:56 2012 +0000
@@ -72,7 +72,7 @@
 
 naps = zeros(n_samp, n_ch, n_ears);
 
-seglen = 256;
+seglen = 441;  % anything should work; this is 20 ms at default fs
 n_segs = ceil(n_samp / seglen);
 
 if nargout > 1
--- a/matlab/bmm/carfac/CARFAC_Run_Linear.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Run_Linear.m	Thu May 24 22:26:56 2012 +0000
@@ -17,8 +17,8 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.
 
-function [naps, CF] = CARFAC_Run_Linear(CF, input_waves, extra_damping)
-% function [naps, CF] = CARFAC_Run_Linear(CF, input_waves, extra_damping)
+function [naps, CF] = CARFAC_Run_Linear(CF, input_waves, relative_undamping)
+% function [naps, CF] = CARFAC_Run_Linear(CF, input_waves, relative_undamping)
 %
 % This function runs the CARFAC; that is, filters a 1 or more channel
 % sound input to make one or more neural activity patterns (naps);
@@ -26,16 +26,17 @@
 % linear (not detected) output.
 
 % only saving one of these, really:
-saved_v_damp_max = CF.ears(1).CAR_coeffs.v_damp_max;
+velocity_scale = CF.ears(1).CAR_coeffs.velocity_scale;
 for ear = 1:CF.n_ears
-  CF.ears(ear).CAR_coeffs.v_damp_max = 0.00;  % make it linear for now
+  % make it effectively linear for now
+  CF.ears(ear).CAR_coeffs.velocity_scale = 0;
 end
 
 [n_samp, n_ears] = size(input_waves);
 n_ch = CF.n_ch;
 
 if nargin < 3
-  extra_damping = 0;
+  relative_undamping = 1;  % default to min-damping condition
 end
 
 if n_ears ~= CF.n_ears
@@ -43,11 +44,11 @@
 end
 
 for ear = 1:CF.n_ears
+  coeffs = CF.ears(ear).CAR_coeffs;
   % Set the state of damping, and prevent interpolation from there:
-  CF.ears(ear).CAR_state.zB_memory(:) = extra_damping;  % interpolator state
+  CF.ears(ear).CAR_state.zB_memory(:) = coeffs.zr_coeffs .* relative_undamping;  % interpolator state
   CF.ears(ear).CAR_state.dzB_memory(:) = 0;  % interpolator slope
-  CF.ears(ear).CAR_state.g_memory = CARFAC_Stage_g( ...
-    CF.ears(ear).CAR_coeffs(ear), extra_damping);
+  CF.ears(ear).CAR_state.g_memory = CARFAC_Stage_g(coeffs, relative_undamping);
   CF.ears(ear).CAR_state.dg_memory(:) = 0;  % interpolator slope
 end
 
@@ -64,6 +65,6 @@
 end
 
 for ear = 1:CF.n_ears
-  CF.ears(ear).CAR_coeffs.v_damp_max = saved_v_damp_max;
+  CF.ears(ear).CAR_coeffs.velocity_scale = velocity_scale;
 end
 
--- a/matlab/bmm/carfac/CARFAC_Run_Segment.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Run_Segment.m	Thu May 24 22:26:56 2012 +0000
@@ -87,7 +87,7 @@
     
     % run the AGC update step, decimating internally,
     [CF.ears(ear).AGC_state, updated] = CARFAC_AGC_Step( ...
-      CF.ears(ear).AGC_coeffs, ihc_out, CF.ears(ear).AGC_state);
+      ihc_out, CF.ears(ear).AGC_coeffs, CF.ears(ear).AGC_state);
     
     % save some output data:
     naps(k, :, ear) = ihc_out;  % output to neural activity pattern
--- a/matlab/bmm/carfac/CARFAC_Stage_g.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Stage_g.m	Thu May 24 22:26:56 2012 +0000
@@ -17,14 +17,14 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.
 
-function g = CARFAC_Stage_g(CAR_coeffs, extra_damping)
-% function g = CARFAC_Stage_g(CAR_coeffs, extra_damping)
+function g = CARFAC_Stage_g(CAR_coeffs, relative_undamping)
+% function g = CARFAC_Stage_g(CAR_coeffs, relative_undamping)
 % Return the stage gain g needed to get unity gain at DC
 
-r1 = CAR_coeffs.r1_coeffs;  % not zero damping, but min damping
+r1 = CAR_coeffs.r1_coeffs;  % at max damping
 a0 = CAR_coeffs.a0_coeffs;
 c0 = CAR_coeffs.c0_coeffs;
 h  = CAR_coeffs.h_coeffs;
 zr = CAR_coeffs.zr_coeffs;
-r  = r1 - zr.*extra_damping;
+r  = r1 + zr .* relative_undamping;
 g  = (1 - 2*r.*a0 + r.^2) ./ (1 - 2*r.*a0 + h.*r.*c0 + r.^2);
--- a/matlab/bmm/carfac/CARFAC_Transfer_Functions.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Transfer_Functions.m	Thu May 24 22:26:56 2012 +0000
@@ -50,7 +50,7 @@
   
   % Now multiply gains from input to output places; use logs?
   log_gains = log(gains);
-  cum_log_gains = cumsum(log_gains);  % accum across cascaded stages  
+  cum_log_gains = cumsum(log_gains);  % accum across cascaded stages
   
   % And figure out which cascade products we want:
   n_ch = CF.n_ch;
@@ -95,6 +95,7 @@
 gains = (numerators * zz) ./ (denominators * zz);
 
 
+
 function [stage_numerators, stage_denominators] = ...
   CARFAC_Rational_Functions(CF, ear)
 % function [stage_z_numerators, stage_z_denominators] = ...
@@ -107,23 +108,23 @@
 
 n_ch = CF.n_ch;
 coeffs = CF.ears(ear).CAR_coeffs;
-min_zeta = CF.CAR_params.min_zeta;
 
 a0 = coeffs.a0_coeffs;
 c0 = coeffs.c0_coeffs;
 zr = coeffs.zr_coeffs;
 
 % get r, adapted if we have state:
-r =  coeffs.r1_coeffs;
+r1 =  coeffs.r1_coeffs;  % max-damping condition
 if isfield(CF.ears(ear), 'CAR_state')
   state = CF.ears(ear).CAR_state;
-  zB = state.zB_memory; % current extra damping
-  r = r - zr .* zB;
+  zB = state.zB_memory; % current delta-r from undamping
+  r = r1 + zB;
 else
-  zB = 0;
+  zB = 0;  % HIGH-level linear condition by default
 end
 
-g = CARFAC_Stage_g(coeffs, zB);
+relative_undamping = zB ./ zr;
+g = CARFAC_Stage_g(coeffs, relative_undamping);
 a = a0 .* r;
 c = c0 .* r;
 r2 = r .* r;
--- a/matlab/bmm/carfac/CARFAC_hacking.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_hacking.m	Thu May 24 22:26:56 2012 +0000
@@ -22,7 +22,10 @@
 clear variables
 
 %%
+
 use_plan_file = 1;
+dB_list = -60:20:40
+
 if use_plan_file
   
   file_signal = wavread('plan.wav');
@@ -30,12 +33,10 @@
   file_signal = file_signal(10000+(1:10000));  % trim for a faster test
   
 else
-  flist = [1000];
-  alist = [1];
-  flist = 1000;
-  alist = 1;
+  flist = 1400 + (1:4)*200;
+  alist = [1, 1, 1, 1];
   sine_signal = 0;
-  times = (0:19999)' / 22050;
+  times = (0:9999)' / 22050;
   for fno = 1:length(flist)
     sine_signal = sine_signal + alist(fno)*sin(flist(fno)*2*pi*times);
   end
@@ -43,16 +44,7 @@
   file_signal = 1.0 * (sine_signal .* (times/max(times)).^growth_power);
 end
 
-% repeat with negated signal to compare responses:
-% file_signal = [file_signal; -file_signal];
-
 % make a long test signal by repeating at different levels:
-do_distortion_figs = 0;  % use 1 for distortion figure hack....
-if do_distortion_figs
-  dB_list = -40:20:0
-else
-  dB_list = -80:20:60
-end
 % dB = dB_list(1);
 % test_signal =  10^(dB/20)* file_signal(1:4000) % lead-in [];
 test_signal = [];
@@ -65,7 +57,7 @@
 
 agc_plot_fig_num = 6;
 
-for n_ears = 1    %  1:2
+for n_ears = 1:2
   
   CF_struct = CARFAC_Design(n_ears);  % default design
   
@@ -82,11 +74,13 @@
   %   nap = deskew(nap);  % deskew doesn't make much difference
   
   %   dB_BM = 10/log(10) * log(filter(1, [1, -0.995], BM(:, 20:50, :).^2));
-  dB_BM = 10/log(10) * log(filter(1, [1, -0.995], BM(:, :, :).^2));
+  sm_BM = filter(1, [1, -0.995], BM(:, :, :).^2);
   
   % only ear 1:
-  MultiScaleSmooth(dB_BM(5000:200:end, :, 1), 1);
+  smoothed = sm_BM(100:100:end, :, 1);
+  MultiScaleSmooth(10/log(10) * log(smoothed), 1);
   
+ 
   % Display results for 1 or 2 ears:
   for ear = 1:n_ears
     smooth_nap = nap_decim(:, :, ear);
@@ -105,52 +99,6 @@
   CF_struct.ears(1).AGC_state
   min_max_decim = [min(nap_decim(:)), max(nap_decim(:))]
   
-  %%
-  if do_distortion_figs
-    channels = [38, 39, 40];
-    times = 0000 + (3200:3599);
-    smoothed_ohc = ohc;
-    epsi = 0.4;
-    %   smoothed_ohc = filter(epsi, [1, -(1-epsi)], smoothed_ohc);
-    %   smoothed_ohc = filter(epsi, [1, -(1-epsi)], smoothed_ohc);
-    
-    figure(101); hold off
-    plot(smoothed_ohc(times, channels))
-    hold on
-    plot(smoothed_ohc(times, channels(2)), 'k-', 'LineWidth', 1.5)
-    title('OHC')
-    
-    figure(105); hold off
-    plot(agc(times, channels))
-    hold on
-    plot(agc(times, channels(2)), 'k-', 'LineWidth', 1.5)
-    title('AGC')
-    
-    figure(103); hold off
-    plot(BM(times, channels))
-    hold on
-    plot(BM(times, channels(2)), 'k-', 'LineWidth', 1.5)
-    title('BM')
-    
-    extra_damping = smoothed_ohc + agc;
-    figure(102); hold off
-    plot(extra_damping(times, channels))
-    hold on
-    plot(extra_damping(times, channels(2)), 'k-', 'LineWidth', 1.5)
-    title('extra damping')
-    
-    distortion = -(extra_damping - smooth1d(extra_damping, 10)) .* BM;
-    distortion = filter(epsi, [1, -(1-epsi)], distortion);
-    distortion = filter(epsi, [1, -(1-epsi)], distortion);
-    %   distortion = filter(epsi, [1, -(1-epsi)], distortion);
-    figure(104); hold off
-    plot(distortion(times, channels))
-    hold on
-    plot(distortion(times, channels(2)), 'k-', 'LineWidth', 1.5)
-    title('distortion')
-    
-  end
-  %%
 end
 
 % Expected result:  Figure 3 looks like figure 2, a tiny bit darker.
--- a/matlab/bmm/carfac/MultiScaleSmooth.m	Sat May 12 04:31:59 2012 +0000
+++ b/matlab/bmm/carfac/MultiScaleSmooth.m	Thu May 24 22:26:56 2012 +0000
@@ -29,8 +29,6 @@
 % Until we decide what we want, we'll just plot things, one plot per scale.
 
 fig_offset1 = 10;
-fig_offset2 = 30;
-fig_offset3 = 50;
 
 if nargin < 2
   n_scales = 20;
@@ -53,9 +51,6 @@
   figure(scale_no + fig_offset1)
   imagesc(squeeze(smoothed(:,:,1))')
 
-  figure(scale_no + fig_offset2)
-  plot(squeeze(mean(smoothed, 2)));
-
   drawnow
   pause(1)