Mercurial > hg > aimc

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_AGCStep.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,82 @@
+% 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 state = CARFAC_AGCStep(AGC_coeffs, avg_detects, state)
+% function state = CARFAC_AGCStep(AGC_coeffs, avg_detects, state)
+%
+% one time step (at decimated low AGC rate) of the AGC state update
+
+n_AGC_stages = length(AGC_coeffs.AGC_epsilon);
+n_mics = length(state);
+n_ch = size(state(1).AGC_sum, 1);  % number of channels
+
+optimize_for_mono = n_mics == 1;  % mono optimization
+if ~optimize_for_mono
+  stage_sum = zeros(n_ch, 1);
+end
+
+for stage = 1:n_AGC_stages
+  if ~optimize_for_mono  % skip if mono
+    if stage > 1
+      prev_stage_mean = stage_sum / n_mics;
+    end
+    stage_sum(:) = 0;  % sum accumulating over mics at this stage
+  end
+  epsilon = AGC_coeffs.AGC_epsilon(stage);  % for this stage's LPF pole
+  polez1 = AGC_coeffs.AGC1_polez(stage);
+  polez2 = AGC_coeffs.AGC2_polez(stage);
+  for mic = 1:n_mics
+    if stage == 1
+      AGC_in = AGC_coeffs.detect_scale * avg_detects(:,mic);
+      AGC_in = max(0, AGC_in);  % don't let neg inputs in
+    else
+      % prev. stage mixed with prev_stage_sum
+      if optimize_for_mono
+        % Mono optimization ignores AGC_mix_coeff,
+        % assuming all(prev_stage_mean == AGC_memory(:, stage - 1));
+        % but we also don't even allocate or compute the sum or mean.
+        AGC_in = AGC_coeffs.AGC_stage_gain * ...
+          state(mic).AGC_memory(:, stage - 1);
+      else
+        AGC_in = AGC_coeffs.AGC_stage_gain * ...
+          (AGC_coeffs.AGC_mix_coeff * prev_stage_mean + ...
+            (1 - AGC_coeffs.AGC_mix_coeff) * ...
+              state(mic).AGC_memory(:, stage - 1));
+      end
+    end
+    AGC_stage = state(mic).AGC_memory(:, stage);
+    % first-order recursive smooting filter update:
+    AGC_stage = AGC_stage + epsilon * (AGC_in - AGC_stage);
+
+    % spatially spread it; using diffusion coeffs like in smooth1d
+    AGC_stage = SmoothDoubleExponential(AGC_stage, polez1, polez2);
+
+    state(mic).AGC_memory(:, stage) = AGC_stage;
+    if stage == 1
+      state(mic).sum_AGC = AGC_stage;
+    else
+      state(mic).sum_AGC = state(mic).sum_AGC + AGC_stage;
+    end
+    if ~optimize_for_mono
+      stage_sum = stage_sum + AGC_stage;
+    end
+  end
+end
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Design.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,342 @@
+% 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 CF = CARFAC_Design(fs, CF_filter_params, ...
+  CF_AGC_params, ERB_break_freq, ERB_Q, CF_IHC_params)
+% function CF = CARFAC_Design(fs, CF_filter_params, ...
+%   CF_AGC_params, ERB_break_freq, ERB_Q, CF_IHC_params)
+%
+% This function designs the CARFAC (Cascade of Asymmetric Resonators with
+% Fast-Acting Compression); that is, it take bundles of parameters and
+% computes all the filter coefficients needed to run it.
+%
+% fs is sample rate (per second)
+% CF_filter_params bundles all the pole-zero filter cascade parameters
+% CF_AGC_params bundles all the automatic gain control parameters
+% CF_IHC_params bundles all the inner hair cell parameters
+%
+% See other functions for designing and characterizing the CARFAC:
+% [naps, CF] = CARFAC_Run(CF, input_waves)
+% transfns = CARFAC_Transfer_Functions(CF, to_channels, from_channels)
+%
+% Defaults to Glasberg & Moore's ERB curve:
+% ERB_break_freq = 1000/4.37;  % 228.833
+% ERB_Q = 1000/(24.7*4.37);    % 9.2645
+%
+% All args are defaultable; for sample/default args see the code; they
+% make 96 channels at default fs = 22050, 114 channels at 44100.
+
+if nargin < 6
+  % HACK: these constant control the defaults
+  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
+  else
+    if one_cap
+      CF_IHC_params = struct( ...
+        'just_hwr', 0, ...        % not just a simple HWR
+        '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
+    else
+      CF_IHC_params = struct( ...
+        'just_hwr', 0, ...        % not just a simple HWR
+        'one_cap', one_cap, ...   % bool; 0 for new two-cap hack
+        'tau_lpf', 0.000080, ...  % 80 microseconds smoothing twice
+        'tau1_out', 0.020, ...    % depletion tau is pretty fast
+        'tau1_in', 0.020, ...     % recovery tau is slower
+        'tau2_out', 0.005, ...   % depletion tau is pretty fast
+        'tau2_in', 0.005 );        % recovery tau is slower
+    end
+  end
+end
+
+if nargin < 5
+  %	Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138
+  % ERB = 24.7 * (1 + 4.37 * CF_Hz / 1000);
+  ERB_Q = 1000/(24.7*4.37);  % 9.2645
+  if nargin < 4
+    ERB_break_freq = 1000/4.37;  % 228.833
+  end
+end
+
+if nargin < 3
+  CF_AGC_params = struct( ...
+    'n_stages', 4, ...
+    'time_constants', [1, 4, 16, 64]*0.002, ...
+    'AGC_stage_gain', 2, ...  % gain from each stage to next slower stage
+    'decimation', 16, ...  % how often to update the AGC states
+    'AGC1_scales', [1, 2, 3, 4]*1, ...   % in units of channels
+    'AGC2_scales', [1, 2, 3, 4]*1.25, ... % spread more toward base
+    'detect_scale', 0.15, ...  % the desired damping range
+    'AGC_mix_coeff', 0.25);
+end
+
+if nargin < 2
+  CF_filter_params = struct( ...
+    'velocity_scale', 0.2, ...  % for the cubic nonlinearity
+    'min_zeta', 0.12, ...
+    'first_pole_theta', 0.78*pi, ...
+    'zero_ratio', sqrt(2), ...
+    'ERB_per_step', 0.3333, ... % assume G&M's ERB formula
+    'min_pole_Hz', 40 );
+end
+
+if nargin < 1
+  fs = 22050;
+end
+
+% first figure out how many filter stages (PZFC/CARFAC channels):
+pole_Hz = CF_filter_params.first_pole_theta * fs / (2*pi);
+n_ch = 0;
+while pole_Hz > CF_filter_params.min_pole_Hz
+  n_ch = n_ch + 1;
+  pole_Hz = pole_Hz - CF_filter_params.ERB_per_step * ...
+    ERB_Hz(pole_Hz, ERB_break_freq, ERB_Q);
+end
+% Now we have n_ch, the number of channels, so can make the array
+% and compute all the frequencies again to put into it:
+pole_freqs = zeros(n_ch, 1);
+pole_Hz = CF_filter_params.first_pole_theta * fs / (2*pi);
+for ch = 1:n_ch
+  pole_freqs(ch) = pole_Hz;
+  pole_Hz = pole_Hz - CF_filter_params.ERB_per_step * ...
+    ERB_Hz(pole_Hz, ERB_break_freq, ERB_Q);
+end
+% now we have n_ch, the number of channels, and pole_freqs array
+
+CF = struct( ...
+  'fs', fs, ...
+  'filter_params', CF_filter_params, ...
+  'AGC_params', CF_AGC_params, ...
+  'IHC_params', CF_IHC_params, ...
+  'n_ch', n_ch, ...
+  'pole_freqs', pole_freqs, ...
+  'filter_coeffs', CARFAC_DesignFilters(CF_filter_params, fs, pole_freqs), ...
+  'AGC_coeffs', CARFAC_DesignAGC(CF_AGC_params, fs), ...
+  'IHC_coeffs', CARFAC_DesignIHC(CF_IHC_params, fs), ...
+  'n_mics', 0 );
+
+% adjust the AGC_coeffs to account for IHC saturation level to get right
+% damping change as specified in CF.AGC_params.detect_scale
+CF.AGC_coeffs.detect_scale = CF.AGC_params.detect_scale / ...
+  (CF.IHC_coeffs.saturation_output * CF.AGC_coeffs.AGC_gain);
+
+%% Design the filter coeffs:
+function filter_coeffs = CARFAC_DesignFilters(filter_params, fs, pole_freqs)
+
+n_ch = length(pole_freqs);
+
+% the filter design coeffs:
+
+filter_coeffs = struct('velocity_scale', filter_params.velocity_scale);
+
+filter_coeffs.r_coeffs = zeros(n_ch, 1);
+filter_coeffs.a_coeffs = zeros(n_ch, 1);
+filter_coeffs.c_coeffs = zeros(n_ch, 1);
+filter_coeffs.h_coeffs = zeros(n_ch, 1);
+filter_coeffs.g_coeffs = zeros(n_ch, 1);
+
+% zero_ratio comes in via h.  In book's circuit D, zero_ratio is 1/sqrt(a),
+% and that a is here 1 / (1+f) where h = f*c.
+% solve for f:  1/zero_ratio^2 = 1 / (1+f)
+% zero_ratio^2 = 1+f => f = zero_ratio^2 - 1
+f = filter_params.zero_ratio^2 - 1;  % nominally 1 for half-octave
+
+% Make pole positions, s and c coeffs, h and g coeffs, etc.,
+% which mostly depend on the pole angle theta:
+theta = pole_freqs .* (2 * pi / fs);
+
+% different possible interpretations for min-damping r:
+% r = exp(-theta * CF_filter_params.min_zeta).
+% Using sin gives somewhat higher Q at highest thetas.
+r = (1 - sin(theta) * filter_params.min_zeta);
+filter_coeffs.r_coeffs = r;
+
+% undamped coupled-form coefficients:
+filter_coeffs.a_coeffs = cos(theta);
+filter_coeffs.c_coeffs = sin(theta);
+
+% the zeros follow via the h_coeffs
+h = sin(theta) .* f;
+filter_coeffs.h_coeffs = h;
+
+r2 = r;  % aim for unity DC gain at min damping, here; or could try r^2
+filter_coeffs.g_coeffs = 1 ./ (1 + h .* r2 .* sin(theta) ./ ...
+  (1 - 2 * r2 .* cos(theta) + r2 .^ 2));
+
+
+%% the AGC design coeffs:
+function AGC_coeffs = CARFAC_DesignAGC(AGC_params, fs)
+
+AGC_coeffs = struct('AGC_stage_gain', AGC_params.AGC_stage_gain, ...
+  'AGC_mix_coeff', AGC_params.AGC_mix_coeff);
+
+
+% AGC1 pass is smoothing from base toward apex;
+% AGC2 pass is back, which is done first now
+AGC1_scales = AGC_params.AGC1_scales;
+AGC2_scales = AGC_params.AGC2_scales;
+
+n_AGC_stages = AGC_params.n_stages;
+AGC_coeffs.AGC_epsilon = zeros(1, n_AGC_stages);  % the 1/(tau*fs) roughly
+decim = AGC_params.decimation;
+gain = 0;
+for stage = 1:n_AGC_stages
+  tau = AGC_params.time_constants(stage);
+  % epsilon is how much new input to take at each update step:
+  AGC_coeffs.AGC_epsilon(stage) = 1 - exp(-decim / (tau * fs));
+  % and these are the smoothing scales and poles for decimated rate:
+  ntimes = tau * (fs / decim);  % effective number of smoothings
+  % divide the spatial variance by effective number of smoothings:
+  t = (AGC1_scales(stage)^2) / ntimes;  % adjust scale for diffusion
+  AGC_coeffs.AGC1_polez(stage) = 1 + 1/t - sqrt((1+1/t)^2 - 1);
+  t = (AGC2_scales(stage)^2) / ntimes;  % adjust scale for diffusion
+  AGC_coeffs.AGC2_polez(stage) = 1 + 1/t - sqrt((1+1/t)^2 - 1);
+  gain = gain + AGC_params.AGC_stage_gain^(stage-1);
+end
+
+AGC_coeffs.AGC_gain = gain;
+
+%% the IHC design coeffs:
+function IHC_coeffs = CARFAC_DesignIHC(IHC_params, fs)
+
+if IHC_params.just_hwr
+  IHC_coeffs = struct('just_hwr', 1);
+  IHC_coeffs.saturation_output = 10;  % HACK: assume some max out
+else
+  if IHC_params.one_cap
+    IHC_coeffs = struct(...
+      'just_hwr', 0, ...
+      'lpf_coeff', 1 - exp(-1/(IHC_params.tau_lpf * fs)), ...
+      'out_rate', 1 / (IHC_params.tau_out * fs), ...
+      'in_rate', 1 / (IHC_params.tau_in * fs), ...
+      'one_cap', IHC_params.one_cap);
+  else
+    IHC_coeffs = struct(...
+      'just_hwr', 0, ...
+      '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', 1 / (IHC_params.tau2_out * fs), ...
+      'in2_rate', 1 / (IHC_params.tau2_in * fs), ...
+      'one_cap', IHC_params.one_cap);
+  end
+
+  % run one channel to convergence to get rest state:
+  IHC_coeffs.rest_output = 0;
+  IHC_state = struct( ...
+    'cap_voltage', 0, ...
+    'cap1_voltage', 0, ...
+    'cap2_voltage', 0, ...
+    'lpf1_state', 0, ...
+    'lpf2_state', 0, ...
+    'ihc_accum', 0);
+
+  IHC_in = 0;
+  for k = 1:30000
+    [IHC_out, IHC_state] = CARFAC_IHCStep(IHC_in, IHC_coeffs, IHC_state);
+  end
+
+  IHC_coeffs.rest_output = IHC_out;
+  IHC_coeffs.rest_cap = IHC_state.cap_voltage;
+  IHC_coeffs.rest_cap1 = IHC_state.cap1_voltage;
+  IHC_coeffs.rest_cap2 = IHC_state.cap2_voltage;
+
+  LARGE = 2;
+  IHC_in = LARGE;  % "Large" saturating input to IHC; make it alternate
+  for k = 1:30000
+    [IHC_out, IHC_state] = CARFAC_IHCStep(IHC_in, IHC_coeffs, IHC_state);
+    prev_IHC_out = IHC_out;
+    IHC_in = -IHC_in;
+  end
+
+  IHC_coeffs.saturation_output = (IHC_out + prev_IHC_out) / 2;
+end
+
+%%
+% default design result, running this function with no args, should look
+% like this, before CARFAC_Init puts state storage into it:
+%
+% CF = CARFAC_Design
+% CF.filter_params
+% CF.AGC_params
+% CF.filter_coeffs
+% CF.AGC_coeffs
+% CF.IHC_coeffs
+%
+% CF =
+%                fs: 22050
+%     filter_params: [1x1 struct]
+%        AGC_params: [1x1 struct]
+%        IHC_params: [1x1 struct]
+%              n_ch: 96
+%        pole_freqs: [96x1 double]
+%     filter_coeffs: [1x1 struct]
+%        AGC_coeffs: [1x1 struct]
+%        IHC_coeffs: [1x1 struct]
+%            n_mics: 0
+% ans =
+%       velocity_scale: 0.2000
+%             min_zeta: 0.1200
+%     first_pole_theta: 2.4504
+%           zero_ratio: 1.4142
+%         ERB_per_step: 0.3333
+%          min_pole_Hz: 40
+% ans =
+%           n_stages: 4
+%     time_constants: [0.0020 0.0080 0.0320 0.1280]
+%     AGC_stage_gain: 2
+%         decimation: 16
+%        AGC1_scales: [1 2 3 4]
+%        AGC2_scales: [1.2500 2.5000 3.7500 5]
+%       detect_scale: 0.1500
+%      AGC_mix_coeff: 0.2500
+% ans =
+%     velocity_scale: 0.2000
+%           r_coeffs: [96x1 double]
+%           a_coeffs: [96x1 double]
+%           c_coeffs: [96x1 double]
+%           h_coeffs: [96x1 double]
+%           g_coeffs: [96x1 double]
+% ans =
+%     AGC_stage_gain: 2
+%      AGC_mix_coeff: 0.2500
+%        AGC_epsilon: [0.3043 0.0867 0.0224 0.0057]
+%         AGC1_polez: [0.1356 0.1356 0.0854 0.0417]
+%         AGC2_polez: [0.1872 0.1872 0.1227 0.0623]
+%           AGC_gain: 15
+%       detect_scale: 0.0630
+% ans =
+%             lpf_coeff: 0.4327
+%             out1_rate: 0.0023
+%              in1_rate: 0.0023
+%             out2_rate: 0.0091
+%              in2_rate: 0.0091
+%               one_cap: 0
+%           rest_output: 0.0365
+%              rest_cap: 0
+%             rest_cap1: 0.9635
+%             rest_cap2: 0.9269
+%     saturation_output: 0.1587
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Detect.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,92 @@
+% 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 conductance = CARFAC_Detect(x_in)
+% function conductance = CARFAC_detect(x_in)
+% An IHC-like sigmoidal detection nonlinearity for the CARFAC.
+% Resulting conductance is in about [0...1.3405]
+
+
+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
+
+set = x_in > -a;
+z = x_in(set) + a;
+
+% zero is the final answer for many points:
+conductance = zeros(size(x_in));
+conductance(set) = z.^3 ./ (z.^3 + z.^2 + 0.1);
+
+
+%% other things I tried:
+%
+% % zero is the final answer for many points:
+% conductance = zeros(size(x_in));
+%
+% order = 4;  % 3 is a little cheaper; 4 has continuous second deriv.
+%
+% % thresholds and terms involving just a, b, s are scalar ops; x are vectors
+%
+% switch order
+%   case 3
+%     a = 0.15;  % offset of low-end tail into neg x territory
+%     b = 1; % 0.44;   % width of poly segment
+%     slope = 0.7;
+%
+%     threshold1 = -a;
+%     threshold2 = b - a;
+%
+%     set2 = x_in > threshold2;
+%     set1 = x_in > threshold1 & ~set2;
+%
+%     s = slope/(2*b - 3/2*b^2);  % factor to make slope at breakpoint
+%     t = s * (b^2 - (b^3) / 2);
+%
+%     x = x_in(set1) - threshold1;
+%     conductance(set1) = s * x .* (x - x .* x / 2);  % x.^2 - 0.5x.^3
+%
+%     x = x_in(set2) - threshold2;
+%     conductance(set2) = t + slope * x ./ (1 + x);
+%
+%
+%   case 4
+%     a = 0.24;  % offset of low-end tail into neg x territory
+%     b = 0.57;   % width of poly segment; 0.5 to end at zero curvature,
+%     a = 0.18;  % offset of low-end tail into neg x territory
+%     b = 0.57;   % width of poly segment; 0.5 to end at zero curvature,
+%     % 0.57 to approx. match curvature of the upper segment.
+%     threshold1 = -a;
+%     threshold2 = b - a;
+%
+%
+%     set2 = x_in > threshold2;
+%     set1 = x_in > threshold1 & ~set2;
+%
+%     s = 1/(3*b^2 - 4*b^3);  % factor to make slope 1 at breakpoint
+%     t = s * (b^3 - b^4);
+%
+%     x = x_in(set1) - threshold1;
+%     conductance(set1) = s * x .* x .* (x - x .* x);  % x.^3 - x.^4
+%
+%     x = x_in(set2) - threshold2;
+%     conductance(set2) = t + x ./ (1 + x);
+%
+% end
+%
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_FilterStep.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,60 @@
+% 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 [zY, state] = CARFAC_FilterStep(x_in, filter_coeffs, state)
+% function [zY, state] = CARFAC_FilterStep(x_in, filter_coeffs, state)
+%
+% One sample-time update step for the filter part of the CARFAC.
+
+% Use each stage previous Y as input to next:
+inputs = [x_in; state.zY_memory(1:(end-1))];
+
+% Local nonlinearity zA and AGC feedback zB reduce pole radius:
+zA = state.zA_memory;
+zB = state.zB_memory + state.dzB_memory; % AGC interpolation
+r0 = filter_coeffs.r_coeffs;
+r = r0 - filter_coeffs.c_coeffs .* (zA + zB);
+
+% now reduce state by r and rotate with the fixed cos/sin coeffs:
+z1 = r .* (filter_coeffs.a_coeffs .* state.z1_memory - ...
+  filter_coeffs.c_coeffs .* state.z2_memory);
+z1 = z1 + inputs;
+z2 = r .* (filter_coeffs.c_coeffs .* state.z1_memory + ...
+  filter_coeffs.a_coeffs .* state.z2_memory);
+
+% update the "velocity" for cubic nonlinearity, into zA:
+zA = (((state.z2_memory - z2) .* filter_coeffs.velocity_scale) - 0.2) .^ 2;
+
+velocity_damp_max = 1/16;
+zA = velocity_damp_max * zA ./ (1 + zA);  % soft max at velocity_damp_max
+
+% Get outputs from inputs and new z2 values:
+zY = filter_coeffs.g_coeffs .* (inputs + filter_coeffs.h_coeffs .* z2);
+
+% put new state back in place of old
+state.z1_memory = z1;
+state.z2_memory = z2;
+state.zA_memory = zA;
+state.zB_memory = zB;
+state.zY_memory = zY;
+
+% accum the straight hwr version for the sake of AGC range:
+hwr_detect = max(0, zY);  % detect with HWR
+state.detect_accum = state.detect_accum + hwr_detect;
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_IHCStep.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,62 @@
+% 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 [ihc_out, state] = CARFAC_IHCStep(filters_out, coeffs, state);
+% function [ihc_out, state] = CARFAC_IHCStep(filters_out, coeffs, state);
+%
+% 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;
+
+if just_hwr
+  ihc_out = max(0, filters_out);
+  state.ihc_accum = state.ihc_accum + max(0, ihc_out);
+else
+  one_cap = coeffs.one_cap;
+
+  detect = CARFAC_Detect(filters_out/2);  % detect with HWR or so
+
+  if one_cap
+    ihc_out = detect .* state.cap_voltage;
+    state.cap_voltage = state.cap_voltage - ihc_out .* coeffs.out_rate + ...
+      (1 - state.cap_voltage) .* coeffs.in_rate;
+  else
+    % change to 2-cap version more like Meddis's:
+    ihc_out = detect .* state.cap2_voltage;
+    state.cap1_voltage = state.cap1_voltage - ...
+      (state.cap1_voltage - state.cap2_voltage) .* coeffs.out1_rate + ...
+      (1 - state.cap1_voltage) .* coeffs.in1_rate;
+
+    state.cap2_voltage = state.cap2_voltage - ihc_out .* coeffs.out2_rate + ...
+      (state.cap1_voltage - state.cap2_voltage) .* coeffs.in2_rate;
+  end
+
+  % smooth it twice with LPF:
+
+  state.lpf1_state = state.lpf1_state + coeffs.lpf_coeff * ...
+    (ihc_out - state.lpf1_state);
+
+  state.lpf2_state = state.lpf2_state + coeffs.lpf_coeff * ...
+    (state.lpf1_state - state.lpf2_state);
+
+  ihc_out = state.lpf2_state;
+
+  state.ihc_accum = state.ihc_accum + max(0, ihc_out - coeffs.rest_output);
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Init.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,87 @@
+% 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 CF_struct = CARFAC_Init(CF_struct, n_mics)
+% function CF_struct = CARFAC_Init(CF_struct, n_mics)
+%
+% Initialize state for n_mics channels (default 1).
+% This allocates and zeros all the state vector storage in the CF_struct.
+
+% TODO (dicklyon): Review whether storing state in the same struct as
+% the design is a good thing, or whether we want another
+% level of object.  I like fewer structs and class types.
+
+if nargin < 2
+  n_mics = 1;  % monaural
+end
+
+% % this is probably what I'd do in the C++ version:
+% if CF_struct.n_mics ~= n_mics;
+%   % free the state and make new number of channels
+%   % make a struct arrray, one element per mic channel, numbered:
+%   for k = 1:n_mics
+%     CF_struct.state(k)  = struct('mic_number', k);
+%   end
+% end
+% But this code doesn't work because I don't understand struct arrays.
+
+% For now I don't ever free anything if n_mics is reduced;
+% so be sure to respect n_mics, not the size of the state struct array.
+
+AGC_time_constants = CF_struct.AGC_params.time_constants;
+n_AGC_stages = length(AGC_time_constants);
+
+CF_struct.n_mics = n_mics;
+CF_struct.k_mod_decim = 0;  % time index phase, cumulative over segments
+
+% This zeroing grows the struct array as needed:
+for mic = 1:n_mics
+  CF_struct.filter_state(mic).z1_memory = zeros(CF_struct.n_ch, 1);
+  CF_struct.filter_state(mic).z2_memory = zeros(CF_struct.n_ch, 1);
+  % cubic loop
+  CF_struct.filter_state(mic).zA_memory = zeros(CF_struct.n_ch, 1);
+  % AGC interp
+  CF_struct.filter_state(mic).zB_memory = zeros(CF_struct.n_ch, 1);
+  % AGC incr
+  CF_struct.filter_state(mic).dzB_memory = zeros(CF_struct.n_ch, 1);
+  CF_struct.filter_state(mic).zY_memory = zeros(CF_struct.n_ch, 1);
+  CF_struct.filter_state(mic).detect_accum = zeros(CF_struct.n_ch, 1);
+  % AGC loop filters' state:
+  % HACK init
+  CF_struct.AGC_state(mic).AGC_memory = zeros(CF_struct.n_ch, n_AGC_stages);
+  CF_struct.AGC_state(mic).AGC_sum = zeros(CF_struct.n_ch, 1);
+  % IHC state:
+  if CF_struct.IHC_coeffs.just_hwr
+    CF_struct.IHC_state(mic).ihc_accum = zeros(CF_struct.n_ch, 1);
+  else
+    CF_struct.IHC_state(mic).cap_voltage = ...
+      CF_struct.IHC_coeffs(mic).rest_cap * ones(CF_struct.n_ch, 1);
+    CF_struct.IHC_state(mic).cap1_voltage = ...
+      CF_struct.IHC_coeffs(mic).rest_cap1 * ones(CF_struct.n_ch, 1);
+    CF_struct.IHC_state(mic).cap2_voltage = ...
+      CF_struct.IHC_coeffs(mic).rest_cap2 * ones(CF_struct.n_ch, 1);
+    CF_struct.IHC_state(mic).lpf1_state = ...
+      CF_struct.IHC_coeffs(mic).rest_output * zeros(CF_struct.n_ch, 1);
+    CF_struct.IHC_state(mic).lpf2_state = ...
+      CF_struct.IHC_coeffs(mic).rest_output * zeros(CF_struct.n_ch, 1);
+    CF_struct.IHC_state(mic).ihc_accum = zeros(CF_struct.n_ch, 1);
+  end
+end
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Run.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,148 @@
+% 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 [naps, CF, decim_naps] = CARFAC_Run ...
+  (CF, input_waves, AGC_plot_fig_num)
+% function [naps, CF, CF.cum_k, decim_naps] = CARFAC_Run ...
+%    (CF, input_waves, CF.cum_k, AGC_plot_fig_num)
+% This function runs the CARFAC; that is, filters a 1 or more channel
+% sound input to make one or more neural activity patterns (naps).
+%
+% The CF struct holds the filterbank design and state; if you want to
+% break the input up into segments, you need to use the updated CF
+% to keep the state between segments.
+%
+% input_waves is a column vector if there's just one audio channel;
+% more generally, it has a row per time sample, a column per audio channel.
+%
+% naps has a row per time sample, a column per filterbank channel, and
+% a layer per audio channel if more than 1.
+% decim_naps is like naps but time-decimated by the int CF.decimation.
+%
+% the input_waves are assumed to be sampled at the same rate as the
+% CARFAC is designed for; a resampling may be needed before calling this.
+%
+% The function works as an outer iteration on time, updating all the
+% filters and AGC states concurrently, so that the different channels can
+% interact easily.  The inner loops are over filterbank channels, and
+% this level should be kept efficient.
+%
+% See other functions for designing and characterizing the CARFAC:
+% CF = CARFAC_Design(fs, CF_filter_params, CF_AGC_params, n_mics)
+% transfns = CARFAC_Transfer_Functions(CF, to_chans, from_chans)
+
+[n_samp, n_mics] = size(input_waves);
+n_ch = CF.n_ch;
+
+if nargin < 3
+  AGC_plot_fig_num = 0;
+end
+
+if n_mics ~= CF.n_mics
+  error('bad number of input_waves channels passed to CARFAC_Run')
+end
+
+% pull coeffs out of struct first, into local vars for convenience
+decim = CF.AGC_params.decimation;
+
+naps = zeros(n_samp, n_ch, n_mics);
+if nargout > 2
+  % make decimated detect output:
+  decim_naps = zeros(ceil(n_samp/decim), CF.n_ch, CF.n_mics);
+else
+  decim_naps = [];
+end
+
+decim_k = 0;
+
+sum_abs_response = 0;
+
+for k = 1:n_samp
+  CF.k_mod_decim = mod(CF.k_mod_decim + 1, decim);  % global time phase
+  % at each time step, possibly handle multiple channels
+  for mic = 1:n_mics
+    [filters_out, CF.filter_state(mic)] = CARFAC_FilterStep( ...
+      input_waves(k, mic), CF.filter_coeffs, CF.filter_state(mic));
+
+    % update IHC state & output on every time step, too
+    [ihc_out, CF.IHC_state(mic)] = CARFAC_IHCStep( ...
+      filters_out, CF.IHC_coeffs, CF.IHC_state(mic));
+
+%     sum_abs_response = sum_abs_response + abs(filters_out);
+
+    naps(k, :, mic) = ihc_out;  % output to neural activity pattern
+  end
+
+  % conditionally update all the AGC stages and channels now:
+  if CF.k_mod_decim == 0
+    % just for the plotting option:
+    decim_k = decim_k + 1;   % index of decimated signal for display
+    if ~isempty(decim_naps)
+      for mic = 1:n_mics
+        % this is HWR out of filters, not IHCs
+        avg_detect = CF.filter_state(mic).detect_accum / decim;
+        % This HACK is the IHC version:
+        avg_detect = CF.IHC_state(mic).ihc_accum / decim;  % for cochleagram
+        decim_naps(decim_k, :, mic) = avg_detect;  % for cochleagram
+%         decim_naps(decim_k, :, mic) = sum_abs_response / decim;  % HACK for mechanical out ABS
+%         sum_abs_response(:) = 0;
+      end
+    end
+
+    % get the avg_detects to connect filter_state to AGC_state:
+    avg_detects = zeros(n_ch, n_mics);
+    for mic = 1:n_mics
+%       % mechanical response from filter output through HWR as AGC in:
+%       avg_detects(:, mic) = CF.filter_state(mic).detect_accum / decim;
+      CF.filter_state(mic).detect_accum(:) = 0;  % zero the detect accumulator
+      % New HACK, IHC output relative to rest as input to AGC:
+      avg_detects(:, mic) = CF.IHC_state(mic).ihc_accum / decim;
+      CF.IHC_state(mic).ihc_accum(:) = 0;  % zero the detect accumulator
+    end
+
+    % run the AGC update step:
+    CF.AGC_state = CARFAC_AGCStep(CF.AGC_coeffs, avg_detects, CF.AGC_state);
+
+    % connect the feedback from AGC_state to filter_state:
+    for mic = 1:n_mics
+      new_damping = CF.AGC_state(mic).sum_AGC;
+%       max_damping = 0.15;  % HACK
+%       new_damping = min(new_damping, max_damping);
+      % set the delta needed to get to new_damping:
+      CF.filter_state(mic).dzB_memory = ...
+        (new_damping - CF.filter_state(mic).zB_memory) ...
+          / decim;
+    end
+
+    if AGC_plot_fig_num
+      figure(AGC_plot_fig_num); hold off
+      maxsum = 0;
+      for mic = 1:n_mics
+        plot(CF.AGC_state(mic).AGC_memory)
+        agcsum = sum(CF.AGC_state(mic).AGC_memory, 2);
+        maxsum(mic) = max(maxsum, max(agcsum));
+        hold on
+        plot(agcsum, 'k-')
+      end
+      axis([0, CF.n_ch, 0, max(0.001, maxsum)]);
+      drawnow
+    end
+  end
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_SAI.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,70 @@
+% 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 [CF, sai] = CARFAC_SAI(CF, k, n_mics, naps, sai)
+% function sai = CARFAC_SAI(CF_struct, n_mics, naps, sai)
+%
+% Calculate the Stabilized Auditory Image from naps
+%
+
+  threshold_alpha = CF.sai_params.threshold_alpha;
+  threshold_jump = CF.sai_params.threshold_jump_factor;
+  threshold_offset = CF.sai_params.threshold_jump_offset;
+
+  sai2 = reshape(sai,CF.sai_params.sai_width * CF.n_ch,n_mics);
+  naps2 = reshape(naps,CF.n_samp * CF.n_ch,n_mics);
+
+  for mic = 1:n_mics
+    data = naps(k, :, mic)';
+    above_threshold = (CF.sai_state(mic).lastdata > ...
+                       CF.sai_state(mic).thresholds) & ...
+                      (CF.sai_state(mic).lastdata > data);
+    CF.sai_state(mic).thresholds(above_threshold) = ...
+        data(above_threshold) * threshold_jump + threshold_offset;
+    CF.sai_state(mic).thresholds(~above_threshold) = ...
+        CF.sai_state(mic).thresholds(~above_threshold) * threshold_alpha;
+    CF.sai_state(mic).lastdata = data;
+
+    % Update SAI image with strobe data.
+    othermic = 3 - mic;
+
+    % Channels that are above the threhsold
+    above_ch = find(above_threshold);
+
+    % If we are above the threshold, set the trigger index and reset the
+    % sai_index
+    CF.sai_state(mic).trigger_index(above_ch) = k;
+    CF.sai_state(mic).sai_index(above_ch) = 1;
+
+    % Copy the right data from the nap to the sai
+    chans = (1:CF.n_ch)';
+    fromindices = CF.sai_state(mic).trigger_index() + (chans - 1) * CF.n_samp;
+    toindices = min((CF.sai_state(mic).sai_index() + (chans - 1) * ...
+                     CF.sai_params.sai_width), ...
+                     CF.sai_params.sai_width * CF.n_ch);
+    sai2(toindices,mic) = naps2(fromindices,othermic);
+
+    CF.sai_state(mic).trigger_index(:) = CF.sai_state(mic).trigger_index(:) + 1;
+    CF.sai_state(mic).sai_index(:) = CF.sai_state(mic).sai_index(:) + 1;
+
+  end
+
+  sai = reshape(sai2,CF.sai_params.sai_width,CF.n_ch,n_mics);
+  naps = reshape(naps2,CF.n_samp, CF.n_ch,n_mics);
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_binaural.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,60 @@
+% 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.
+
+% Test/demo hacking for audio/hearing/filterbanks/carfac/matlab/ stuff
+
+clear variables
+
+agc_plot_fig_num = 1;
+
+tic
+
+%test_signal = wavread('pbav4p-22050.wav');
+test_signal = wavread('plan-twochannel-1ms.wav')
+%test_signal = wavread('binaural.wav')
+%test_signal = test_signal(20000:2);  % trim for a faster test
+
+CF_struct = CARFAC_Design;  % default design
+cum_k = 0;  % not doing segments, so just count from 0
+
+% Run mono, then stereo test:
+n_mics = 2
+CF_struct = CARFAC_Init(CF_struct, n_mics);
+
+[nap, CF_struct, cum_k, nap_decim] = ...
+    CARFAC_Run(CF_struct, test_signal, cum_k, agc_plot_fig_num);
+
+% Display results for 1 or 2 mics:
+for mic = 1:n_mics
+  smooth_nap = nap_decim(:, :, mic);
+  figure(mic + n_mics)  % Makes figures 2, 3, and 4
+  image(63 * ((smooth_nap)' .^ 0.5))
+
+  colormap(1 - gray);
+end
+
+toc
+
+% Show resulting data, even though M-Lint complains:
+cum_k
+CF_struct
+CF_struct.filter_state
+CF_struct.AGC_state
+min_max = [min(nap(:)), max(nap(:))]
+min_max_decim = [min(nap_decim(:)), max(nap_decim(:))]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_hacking.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,85 @@
+% 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.
+
+%% Test/demo hacking for CARFAC Matlab stuff:
+
+clear variables
+
+%%
+file_signal = wavread('plan.wav');
+
+% file_signal = file_signal(9000+(1:10000));  % trim for a faster test
+file_signal = file_signal(9300+(1:5000));  % trim for a faster test
+
+% repeat with negated signal to compare responses:
+file_signal = [file_signal; -file_signal];
+
+% make a long test signal by repeating at different levels:
+test_signal = [];
+for dB = -60:20:40  % -80:20:60
+  test_signal = [test_signal; file_signal * 10^(dB/20)];
+end
+
+%%
+CF_struct = CARFAC_Design;  % default design
+
+%% Run mono, then stereo test:
+
+agc_plot_fig_num = 6;
+
+for n_mics = 1  % :2
+  CF_struct = CARFAC_Init(CF_struct, n_mics);
+
+  [nap, CF_struct, nap_decim] = CARFAC_Run(CF_struct, ...
+    test_signal, agc_plot_fig_num);
+
+%   nap = deskew(nap);  % deskew doesn't make much difference
+
+  MultiScaleSmooth(nap_decim, 10);
+
+%   nap_decim = nap;
+%   nap_decim = smooth1d(nap_decim, 1);
+%   nap_decim = nap_decim(1:8:size(nap_decim, 1), :);
+
+  % Display results for 1 or 2 mics:
+  for mic = 1:n_mics
+    smooth_nap = nap_decim(:, :, mic);
+    if n_mics == 1
+      mono_max = max(smooth_nap(:));
+    end
+    figure(3 + mic + n_mics)  % Makes figures 5, ...
+    image(63 * ((max(0, smooth_nap)/mono_max)' .^ 0.5))
+    title('smooth nap from nap decim')
+    colormap(1 - gray);
+  end
+
+  % Show resulting data, even though M-Lint complains:
+  CF_struct
+  CF_struct.k_mod_decim
+  CF_struct.filter_state
+  CF_struct.AGC_state
+  min_max = [min(nap(:)), max(nap(:))]
+  min_max_decim = [min(nap_decim(:)), max(nap_decim(:))]
+
+  % For the 2-channel pass, add a silent second channel:
+  test_signal = [test_signal, zeros(size(test_signal))];
+end
+
+% Expected result:  Figure 3 looks like figure 2, a tiny bit darker.
+% and figure 4 is empty (all zero)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/ERB_Hz.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,15 @@
+function ERB = ERB_Hz(CF_Hz, ERB_break_freq, ERB_Q)
+% function ERB = ERB_Hz(CF_Hz, ERB_break_freq, ERB_Q)
+%
+% Auditory filter nominal Equivalent Rectangular Bandwidth
+%	Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138
+% ERB = 24.7 * (1 + 4.37 * CF_Hz / 1000);
+
+if nargin < 3
+  ERB_Q = 1000/(24.7*4.37);  % 9.2645
+  if nargin < 2
+    ERB_break_freq = 1000/4.37;  % 228.833
+  end
+end
+
+ERB = (ERB_break_freq + CF_Hz) / ERB_Q;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/MultiScaleSmooth.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,70 @@
+% 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 MultiScaleSmooth(waves, n_scales)
+% function MultiScaleSmooth(waves, n_scales)
+%
+% Let's take columns as waveforms, and smooth them to different scales;
+% these inputs can be carfac NAPs, for example, and the peaks of the
+% smoothed versions can be used as trigger events, even tracking back
+% to less-smoothed versions.
+% And we'll deciamte 2:1 at every other smoothing.
+%
+% 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;
+end
+
+smoothed = waves;
+
+for scale_no = 1:n_scales
+
+  if mod(scale_no, 2) == 1
+    newsmoothed = filter([1, 1]/2, 1, smoothed);
+    diffsmoothed = max(0, smoothed - newsmoothed);
+    smoothed = newsmoothed;
+  else
+    newsmoothed = filter([1, 2, 1]/4, 1, smoothed);
+    diffsmoothed = max(0, smoothed - newsmoothed);
+    smoothed = newsmoothed(1:2:end, :);
+  end
+
+  figure(scale_no + fig_offset1)
+  imagesc(smoothed')
+
+  figure(scale_no + fig_offset2)
+  plot(mean(smoothed, 2));
+
+  drawnow
+  pause(1)
+
+end
+
+
+function waves = deskew(waves)
+
+for col = 1:size(waves, 2)
+  waves(1:(end-col+1), col) = waves(col:end, col);
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/SmoothDoubleExponential.m	Wed Feb 15 21:26:40 2012 +0000
@@ -0,0 +1,64 @@
+% 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 signal_vecs = SmoothDoubleExponential(signal_vecs, ...
+  polez1, polez2, fast_matlab_way)
+% function signal_vecs = SmoothDoubleExponential(signal_vecs, ...
+%   polez1, polez2, fast_matlab_way)
+%
+% Smooth the input column vectors in signal_vecs using forward
+% and backwards one-pole smoothing filters, backwards first, with
+% approximately reflecting edge conditions.
+%
+% It will be done with Matlab's filter function if "fast_matlab_way"
+% is nonzero or defaulted; use 0 to test the algorithm for how to do it
+% in sequential c code.
+
+if nargin < 4
+  fast_matlab_way = 1;
+  % can also use the slow way with explicit loop like we'll do in C++
+end
+
+if fast_matlab_way
+  [junk, Z_state] = filter(1-polez1, [1, -polez1], ...
+    signal_vecs((end-10):end, :));  % initialize state from 10 points
+  [signal_vecs(end:-1:1), Z_state] = filter(1-polez2, [1, -polez2], ...
+    signal_vecs(end:-1:1), Z_state*polez2/polez1);
+  signal_vecs = filter(1-polez1, [1, -polez1], signal_vecs, ...
+    Z_state*polez1/polez2);
+else
+  npts = size(signal_vecs, 1);
+  state = zeros(size(signal_vecs, 2));
+  for index = npts-10:npts
+    input = signal_vecs(index, :);
+    state = state + (1 - polez1) * (input - state);
+  end
+  % smooth backward with polez2, starting with state from above:
+  for index = npts:-1:1
+    input = signal_vecs(index, :);
+    state = state + (1 - polez2) * (input - state);
+    signal_vecs(index, :) = state;
+  end
+  % smooth forward with polez1, starting with state from above:
+  for index = 1:npts
+    input = signal_vecs(index, :);
+    state = state + (1 - polez1) * (input - state);
+    signal_vecs(index, :) = state;
+  end
+end