Mercurial > hg > aimc

changeset 469:bc0618485ad4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reparameterize stage gain g and compressed damping with theta; interpolate g
author dicklyon@google.com
date Sun, 11 Mar 2012 00:31:57 +0000
parents 6fed17027b28
children 8e8381785b2b
files matlab/bmm/carfac/CARFAC_Close_AGC_Loop.m matlab/bmm/carfac/CARFAC_Design.m matlab/bmm/carfac/CARFAC_FilterStep.m matlab/bmm/carfac/CARFAC_Init.m matlab/bmm/carfac/CARFAC_Run.m matlab/bmm/carfac/CARFAC_Run_Linear.m matlab/bmm/carfac/CARFAC_Stage_g.m matlab/bmm/carfac/CARFAC_Transfer_Functions.m
diffstat 8 files changed, 162 insertions(+), 95 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/matlab/bmm/carfac/CARFAC_Close_AGC_Loop.m	Sun Mar 11 00:31:57 2012 +0000
@@ -0,0 +1,35 @@
+% 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_Close_AGC_Loop(CF)
+% function CF = CARFAC_Close_AGC_Loop(CF)
+
+% fastest decimated rate determines interp needed:
+decim1 = CF.AGC_params.decimation(1);
+
+for mic = 1:CF.n_mics
+  extra_damping = CF.AGC_state(mic).AGC_memory(:, 1);  % stage 1 result
+  % Update the target stage gain for the new damping:
+  new_g = CARFAC_Stage_g(CF.filter_coeffs(mic), extra_damping);
+  % set the deltas needed to get to the new damping:
+  CF.filter_state(mic).dzB_memory = ...
+    (extra_damping - CF.filter_state(mic).zB_memory) / decim1;
+  CF.filter_state(mic).dg_memory = ...
+    (new_g - CF.filter_state(mic).g_memory) / decim1;
+end
--- a/matlab/bmm/carfac/CARFAC_Design.m	Sat Mar 10 06:22:56 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Design.m	Sun Mar 11 00:31:57 2012 +0000
@@ -96,9 +96,10 @@
     'v_offset', 0.01, ...  % 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.12, ...
+    'min_zeta', 0.10, ...
     '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
     'ERB_per_step', 0.5, ... % assume G&M's ERB formula
     'min_pole_Hz', 30 );
 end
@@ -159,11 +160,11 @@
   'v_damp_max', filter_params.v_damp_max ...
   );
 
-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.r1_coeffs = zeros(n_ch, 1);
+filter_coeffs.a0_coeffs = zeros(n_ch, 1);
+filter_coeffs.c0_coeffs = zeros(n_ch, 1);
 filter_coeffs.h_coeffs = zeros(n_ch, 1);
-filter_coeffs.g_coeffs = zeros(n_ch, 1);
+filter_coeffs.g0_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.
@@ -175,31 +176,35 @@
 % which mostly depend on the pole angle theta:
 theta = pole_freqs .* (2 * pi / fs);
 
+c0 = sin(theta);
+a0 = cos(theta);
+
 % different possible interpretations for min-damping r:
 % r = exp(-theta * CF_filter_params.min_zeta).
-% Using sin gives somewhat higher Q at highest thetas.
-ff = 5;  % fudge factor for theta distortion; at least 1.0
-r = (1 - ff*sin(theta/ff) * filter_params.min_zeta);
-filter_coeffs.r_coeffs = r;
+% Compress theta to give somewhat higher Q at highest thetas:
+ff = filter_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.
+filter_coeffs.zr_coeffs = zr_coeffs;  % how r relates to zeta
+
+r = (1 - zr_coeffs * filter_params.min_zeta);
+filter_coeffs.r1_coeffs = r;
 
 % undamped coupled-form coefficients:
-filter_coeffs.a_coeffs = cos(theta);
-filter_coeffs.c_coeffs = sin(theta);
+filter_coeffs.a0_coeffs = a0;
+filter_coeffs.c0_coeffs = c0;
 
 % the zeros follow via the h_coeffs
-h = sin(theta) .* f;
+h = c0 .* f;
 filter_coeffs.h_coeffs = h;
 
-% % unity gain at min damping, radius r:
-g = (1 - 2*r.*cos(theta) + r.^2) ./ ...
-  (1 - 2*r .* cos(theta) + h .* r .* sin(theta) + r.^2);
-% or assume r is 1, for the zero-damping gain g0:
-g0 = (2 - 2*cos(theta)) ./ ...
-  (2 - 2 * cos(theta) + h .* sin(theta));
+% for unity gain at min damping, radius r; only used in CARFAC_Init:
+extra_damping = zeros(size(r));
+% this function needs to take filter_coeffs even if we haven't finished
+% constucting it by putting in the g0_coeffs:
+filter_coeffs.g0_coeffs = CARFAC_Stage_g(filter_coeffs, extra_damping);
 
-filter_coeffs.g_coeffs = g0;
-% make coeffs that can correct g0 to make g based on (1 - r).^2:
-filter_coeffs.gr_coeffs = ((g ./ g0) - 1) ./ ((1 - r).^2);
 
 %% the AGC design coeffs:
 function AGC_coeffs = CARFAC_DesignAGC(AGC_params, fs)
@@ -385,57 +390,62 @@
 % 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
-%             v_offset: 0.0100
-%            v2_corner: 0.2000
-%           v_damp_max: 0.0100
-%             min_zeta: 0.1200
-%     first_pole_theta: 2.4504
-%           zero_ratio: 1.4142
-%         ERB_per_step: 0.3333
-%          min_pole_Hz: 40
-% ans =
+% CF = 
+%                          fs: 22050
+%     max_channels_per_octave: 12.1873
+%               filter_params: [1x1 struct]
+%                  AGC_params: [1x1 struct]
+%                  IHC_params: [1x1 struct]
+%                        n_ch: 66
+%                  pole_freqs: [66x1 double]
+%               filter_coeffs: [1x1 struct]
+%                  AGC_coeffs: [1x1 struct]
+%                  IHC_coeffs: [1x1 struct]
+%                      n_mics: 0
+% ans = 
+%                 velocity_scale: 0.2000
+%                       v_offset: 0.0100
+%                      v2_corner: 0.2000
+%                     v_damp_max: 0.0100
+%                       min_zeta: 0.1200
+%               first_pole_theta: 2.6704
+%                     zero_ratio: 1.4142
+%     high_f_damping_compression: 0.5000
+%                   ERB_per_step: 0.5000
+%                    min_pole_Hz: 30
+% ans = 
 %           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 2 4 8]
-%        AGC2_scales: [1.5000 3 6 12]
+%        AGC1_scales: [1 2 4 6]
+%        AGC2_scales: [1.5000 3 6 9]
 %       detect_scale: 0.1500
-%      AGC_mix_coeff: 0.3500
-% ans =
+%      AGC_mix_coeff: 0.5000
+% ans = 
 %     velocity_scale: 0.2000
 %           v_offset: 0.0100
 %          v2_corner: 0.2000
 %         v_damp_max: 0.0100
-%           r_coeffs: [96x1 double]
-%           a_coeffs: [96x1 double]
-%           c_coeffs: [96x1 double]
-%           h_coeffs: [96x1 double]
-%           g_coeffs: [96x1 double]
-% ans =
+%          r1_coeffs: [66x1 double]
+%          a0_coeffs: [66x1 double]
+%          c0_coeffs: [66x1 double]
+%           h_coeffs: [66x1 double]
+%          g0_coeffs: [66x1 double]
+%          zr_coeffs: [66x1 double]
+% ans = 
 %             AGC_stage_gain: 2
 %                AGC_epsilon: [0.1659 0.0867 0.0443 0.0224]
 %                 decimation: [8 2 2 2]
-%     AGC_spatial_iterations: [1 1 2 3]
+%                 AGC_polez1: [0.1627 0.2713 0.3944 0.4194]
+%                 AGC_polez2: [0.2219 0.3165 0.4260 0.4414]
+%     AGC_spatial_iterations: [1 1 2 2]
 %            AGC_spatial_FIR: [3x4 double]
 %                 AGC_n_taps: [3 5 5 5]
-%             AGC_mix_coeffs: [0 0.0317 0.0159 0.0079]
+%             AGC_mix_coeffs: [0 0.0454 0.0227 0.0113]
 %                   AGC_gain: 15
 %               detect_scale: 0.0664
-% ans =
+% ans = 
 %              just_hwr: 0
 %             lpf_coeff: 0.4327
 %             out1_rate: 0.0023
@@ -449,4 +459,3 @@
 %             rest_cap2: 0.9269
 %     saturation_output: 0.1507
 
-
--- a/matlab/bmm/carfac/CARFAC_FilterStep.m	Sat Mar 10 06:22:56 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_FilterStep.m	Sun Mar 11 00:31:57 2012 +0000
@@ -27,20 +27,21 @@
 % 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;
+r1 = filter_coeffs.r1_coeffs;
+g = state.g_memory + state.dg_memory;  % interp g
 v_offset  = filter_coeffs.v_offset;
 v2_corner = filter_coeffs.v2_corner;
 v_damp_max = filter_coeffs.v_damp_max;
 
-% zB and zA are "extra damping", and multiply c or sin(theta):
-r = r0 - filter_coeffs.c_coeffs .* (zA + zB); 
+% zB and zA are "extra damping", and multiply zr (compressed theta):
+r = r1 - filter_coeffs.zr_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 = r .* (filter_coeffs.a0_coeffs .* state.z1_memory - ...
+  filter_coeffs.c0_coeffs .* state.z2_memory);
 % z1 = z1 + inputs;
-z2 = r .* (filter_coeffs.c_coeffs .* state.z1_memory + ...
-  filter_coeffs.a_coeffs .* state.z2_memory);
+z2 = r .* (filter_coeffs.c0_coeffs .* state.z1_memory + ...
+  filter_coeffs.a0_coeffs .* state.z2_memory);
 
 % update the "velocity" for cubic nonlinearity, into zA:
 zA = (((state.z2_memory - z2) .* filter_coeffs.velocity_scale) + ...
@@ -48,31 +49,27 @@
 % soft saturation to make it more like an "essential" nonlinearity:
 zA = v_damp_max * zA ./ (v2_corner + zA);
 
-% Adjust gain for r variation:
-g = filter_coeffs.g_coeffs;
-g = g .* (1 + filter_coeffs.gr_coeffs .* (1 - r).^2);
+zY = filter_coeffs.h_coeffs .* z2;  % partial output
 
-gh = g .* filter_coeffs.h_coeffs;
-zY = gh .* z2;  % partial output; still need to ripple in_out
-% ripples input-output path instead of parallel, to avoid delay...
-% this is the only path that doesn't get computed "in parallel":
+% Ripple input-output path, instead of parallel, to avoid delay...
+% this is the only part that doesn't get computed "in parallel":
 in_out = x_in;
 for ch = 1:length(zY)
   % could do this here, or later in parallel:
   z1(ch) = z1(ch) + in_out;
-  % ripple, saving output in zY
-  in_out = g(ch) * in_out + zY(ch);
+  % ripple, saving final channel outputs in zY
+  in_out = g(ch) * (in_out + zY(ch));
   zY(ch) = in_out;
 end
-% % final parallel step is the effect of inputs on state z1:
-% z1 = z1 + [x_in; zY(1:(end-1))];
 
 % put new state back in place of old
+% (z1 and z2 are genuine temps; the others can update by reference in C)
 state.z1_memory = z1;
 state.z2_memory = z2;
 state.zA_memory = zA;
 state.zB_memory = zB;
 state.zY_memory = zY;
+state.g_memory = g;
 
 % accum the straight hwr version in case someone wants it:
 hwr_detect = max(0, zY);  % detect with HWR
--- a/matlab/bmm/carfac/CARFAC_Init.m	Sat Mar 10 06:22:56 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Init.m	Sun Mar 11 00:31:57 2012 +0000
@@ -48,7 +48,6 @@
 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
 n_ch = CF_struct.n_ch;
 
 % keep all the decimator phase info in mic 1 state only:
@@ -63,6 +62,9 @@
   CF_struct.filter_state(mic).dzB_memory = zeros(n_ch, 1);  % AGC incr
   CF_struct.filter_state(mic).zY_memory = zeros(n_ch, 1);
   CF_struct.filter_state(mic).detect_accum = zeros(n_ch, 1);
+  CF_struct.filter_state(mic).g_memory = ...
+    CF_struct.filter_coeffs(mic).g0_coeffs;  % initial g for min_zeta
+  CF_struct.filter_state(mic).dg_memory = zeros(n_ch, 1);    % g interp
   % AGC loop filters' state:
   CF_struct.AGC_state(mic).AGC_memory = zeros(n_ch, n_AGC_stages);  % HACK init
   CF_struct.AGC_state(mic).input_accum = zeros(n_ch, n_AGC_stages);  % HACK init
--- a/matlab/bmm/carfac/CARFAC_Run.m	Sat Mar 10 06:22:56 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Run.m	Sun Mar 11 00:31:57 2012 +0000
@@ -58,9 +58,6 @@
   error('bad number of input_waves channels passed to CARFAC_Run')
 end
 
-% fastest decimated rate determines some interp needed:
-decim1 = CF.AGC_params.decimation(1);
-
 naps = zeros(n_samp, n_ch, n_mics);
 decim_k = 0;
 k_NAP_decim = 0;
@@ -79,7 +76,7 @@
 
 detects = zeros(n_ch, n_mics);
 for k = 1:n_samp
-  CF.k_mod_decim = mod(CF.k_mod_decim + 1, decim1);  % global time phase
+%   CF.k_mod_decim = mod(CF.k_mod_decim + 1, decim1);  % global time phase
   k_NAP_decim = mod(k_NAP_decim + 1, NAP_decim);  % phase of decimated nap
   % at each time step, possibly handle multiple channels
   for mic = 1:n_mics
@@ -110,14 +107,7 @@
   
   % connect the feedback from AGC_state to filter_state when it updates
   if updated
-    for mic = 1:n_mics
-      new_damping = CF.AGC_state(mic).AGC_memory(:, 1);  % stage 1 result
-      % set the delta needed to get to new_damping:
-      % TODO: update this to use da and dc instead of dr maybe?
-      CF.filter_state(mic).dzB_memory = ...
-        (new_damping - CF.filter_state(mic).zB_memory) ...
-        / decim1;
-    end
+    CF = CARFAC_Close_AGC_Loop(CF);
   end
   
   k_AGC = mod(k_AGC + 1, AGC_plot_decim);
--- a/matlab/bmm/carfac/CARFAC_Run_Linear.m	Sat Mar 10 06:22:56 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Run_Linear.m	Sun Mar 11 00:31:57 2012 +0000
@@ -36,9 +36,12 @@
 end
 
 for mic = 1:CF.n_mics
-  % for the state of the AGC interpolator:
+  % Set the state of damping, and prevent interpolation from there:
   CF.filter_state(mic).zB_memory(:) = extra_damping;  % interpolator state
   CF.filter_state(mic).dzB_memory(:) = 0;  % interpolator slope
+  CF.filter_state(mic).g_memory = CARFAC_Stage_g( ...
+    CF.filter_coeffs(mic), extra_damping);
+  CF.filter_state(mic).dg_memory(:) = 0;  % interpolator slope
 end
 
 naps = zeros(n_samp, n_ch, n_mics);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/matlab/bmm/carfac/CARFAC_Stage_g.m	Sun Mar 11 00:31:57 2012 +0000
@@ -0,0 +1,31 @@
+% 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 g = CARFAC_Stage_g(filter_coeffs, extra_damping)
+% function g = CARFAC_Stage_g(filter_coeffs, extra_damping)
+% Return the stage gain g needed to get unity gain at DC
+
+r1 = filter_coeffs.r1_coeffs;  % not zero damping, but min damping
+a0 = filter_coeffs.a0_coeffs;
+c0 = filter_coeffs.c0_coeffs;
+h  = filter_coeffs.h_coeffs;
+zr = filter_coeffs.zr_coeffs;
+r  = r1 - zr.*extra_damping;  % HACK??? or use sin(ff*theta) instead of c?
+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 Mar 10 06:22:56 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Transfer_Functions.m	Sun Mar 11 00:31:57 2012 +0000
@@ -97,25 +97,25 @@
 coeffs = CF.filter_coeffs;
 min_zeta = CF.filter_params.min_zeta;
 
-a0 = coeffs.a_coeffs;
-c0 = coeffs.c_coeffs;
+a0 = coeffs.a0_coeffs;
+c0 = coeffs.c0_coeffs;
+zr = coeffs.zr_coeffs;
 
 % get r, adapted if we have state:
-r =  coeffs.r_coeffs;
+r =  coeffs.r1_coeffs;
 if isfield(CF, 'filter_state')
   state = CF.filter_state;
   zB = state.zB_memory; % current extra damping
-  r = r - c0 .* zB;
+  r = r - zr .* zB;
 else
   zB = 0;
 end
 
+g = CARFAC_Stage_g(coeffs, zB);
 a = a0 .* r;
 c = c0 .* r;
 r2 = r .* r;
 h = coeffs.h_coeffs;
-g0 = coeffs.g_coeffs;
-g = g0 .* (1 + coeffs.gr_coeffs .* (1 - r).^2);
 
 stage_denominators = [ones(n_ch, 1), -2 * a, r2];
 stage_numerators = [g .* ones(n_ch, 1), g .* (-2 * a + h .* c), g .* r2];