Mercurial > hg > aimc
changeset 523:2b96cb7ea4f7
Major AGC improvements mostly
| author | dicklyon@google.com |
|---|---|
| date | Thu, 01 Mar 2012 19:49:24 +0000 |
| parents | c3c85000f804 |
| children | 58d7d67bd138 |
| files | trunk/matlab/bmm/carfac/CARFAC_AGCStep.m trunk/matlab/bmm/carfac/CARFAC_Design.m trunk/matlab/bmm/carfac/CARFAC_FilterStep.m trunk/matlab/bmm/carfac/CARFAC_IHCStep.m trunk/matlab/bmm/carfac/CARFAC_Init.m trunk/matlab/bmm/carfac/CARFAC_Run.m trunk/matlab/bmm/carfac/CARFAC_Spatial_Smooth.m trunk/matlab/bmm/carfac/MultiScaleSmooth.m trunk/matlab/bmm/carfac/SmoothDoubleExponential.m |
| diffstat | 9 files changed, 329 insertions(+), 168 deletions(-) [+] |
line wrap: on
line diff
--- a/trunk/matlab/bmm/carfac/CARFAC_AGCStep.m Mon Feb 27 21:50:20 2012 +0000 +++ b/trunk/matlab/bmm/carfac/CARFAC_AGCStep.m Thu Mar 01 19:49:24 2012 +0000 @@ -17,66 +17,99 @@ % 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) +function [state, updated] = CARFAC_AGCStep(AGC_coeffs, detects, state) +% function [state, updated] = CARFAC_AGCStep(AGC_coeffs, 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 +[n_ch, n_AGC_stages] = size(state(1).AGC_memory); % number of channels optimize_for_mono = n_mics == 1; % mono optimization -if ~optimize_for_mono - stage_sum = zeros(n_ch, 1); + +stage = 1; +ins = AGC_coeffs.detect_scale * detects; +[state, updated] = CARFAC_AGC_Recurse(AGC_coeffs, ins, n_AGC_stages, ... + n_mics, n_ch, optimize_for_mono, stage, state); + + + + + +function [state, updated] = CARFAC_AGC_Recurse(coeffs, ins, n_stages, ... + n_mics, n_ch, mono, stage, state) +% function [state, updated = CARFAC_AGC_Recurse(coeffs, ins, n_stages, ... +% n_mics, n_ch, mono, stage, state) + +decim = coeffs.decimation(stage); % decim phase for this stage +decim_phase = mod(state(1).decim_phase(stage) + 1, decim); +state(1).decim_phase(stage) = decim_phase; + +% accumulate input for this stage from detect or previous stage: +for mic = 1:n_mics + state(mic).input_accum(:, stage) = ... + state(mic).input_accum(:, stage) + ins(:, mic); end -for stage = 1:n_AGC_stages - if ~optimize_for_mono % skip if mono - if stage > 1 - prev_stage_mean = stage_sum / n_mics; +% nothing else to do if it's not the right decim_phase +if decim_phase == 0 + % do lots of work, at decimated rate + + % decimated inputs for this stage, and to be decimated more for next: + for mic = 1:n_mics + ins(:,mic) = state(mic).input_accum(:, stage) / decim; + state(mic).input_accum(:, stage) = 0; % reset accumulator + end + + if stage < n_stages % recurse to evaluate next stage(s) + state = CARFAC_AGC_Recurse(coeffs, ins, n_stages, ... + n_mics, n_ch, mono, stage+1, state); + end + + epsilon = coeffs.AGC_epsilon(stage); % for this stage's LPF pole + stage_gain = coeffs.AGC_stage_gain; + + for mic = 1:n_mics + AGC_in = ins(:,mic); % the newly decimated input for this mic + +% AGC_in = max(0, AGC_in); % don't let neg inputs in + + % add the latest output (state) of next stage... + if stage < n_stages + AGC_in = AGC_in + stage_gain * state(mic).AGC_memory(:, stage+1); 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); + + AGC_stage_state = state(mic).AGC_memory(:, stage); + % first-order recursive smoothing filter update, in time: + AGC_stage_state = AGC_stage_state + ... + epsilon * (AGC_in - AGC_stage_state); + % spatial smooth: + AGC_stage_state = ... + CARFAC_Spatial_Smooth(coeffs, stage, AGC_stage_state); + % and store the state back (in C++, do it all in place?) + state(mic).AGC_memory(:, stage) = AGC_stage_state; + + if ~mono + if mic == 1 + this_stage_sum = AGC_stage_state; 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)); + this_stage_sum = this_stage_sum + AGC_stage_state; 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 + if ~mono + mix_coeff = coeffs.AGC_mix_coeffs(stage); + if mix_coeff > 0 + this_stage_mean = this_stage_sum / n_mics; + for mic = 1:n_mics + state(mic).AGC_memory(:, stage) = ... + state(mic).AGC_memory(:, stage) + ... + mix_coeff * ... + (this_stage_mean - state(mic).AGC_memory(:, stage)); + end end end + updated = 1; % bool to say we have new state +else + updated = 0; end - -
--- a/trunk/matlab/bmm/carfac/CARFAC_Design.m Mon Feb 27 21:50:20 2012 +0000 +++ b/trunk/matlab/bmm/carfac/CARFAC_Design.m Thu Mar 01 19:49:24 2012 +0000 @@ -51,14 +51,14 @@ else if one_cap CF_IHC_params = struct( ... - 'just_hwr', 0, ... % not just a simple HWR + 'just_hwr', just_hwr, ... % 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 + 'just_hwr', just_hwr, ... % 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 @@ -83,16 +83,19 @@ '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 + 'decimation', [8, 2, 2, 2], ... % how often to update the AGC states + 'AGC1_scales', [1, 2, 4, 8]*1, ... % in units of channels + 'AGC2_scales', [1, 2, 4, 8]*1.5, ... % spread more toward base 'detect_scale', 0.15, ... % the desired damping range - 'AGC_mix_coeff', 0.25); + 'AGC_mix_coeff', 0.5); end if nargin < 2 CF_filter_params = struct( ... - 'velocity_scale', 0.2, ... % for the cubic nonlinearity + 'velocity_scale', 0.2, ... % for the "cubic" velocity nonlinearity + '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, ... 'first_pole_theta', 0.78*pi, ... 'zero_ratio', sqrt(2), ... @@ -147,7 +150,11 @@ % the filter design coeffs: -filter_coeffs = struct('velocity_scale', filter_params.velocity_scale); +filter_coeffs = struct('velocity_scale', filter_params.velocity_scale, ... + 'v_offset', filter_params.v_offset, ... + 'v2_corner', filter_params.v2_corner, ... + 'v_damp_max', filter_params.v_damp_max ... + ); filter_coeffs.r_coeffs = zeros(n_ch, 1); filter_coeffs.a_coeffs = zeros(n_ch, 1); @@ -187,9 +194,7 @@ %% 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); - +AGC_coeffs = struct('AGC_stage_gain', AGC_params.AGC_stage_gain); % AGC1 pass is smoothing from base toward apex; % AGC2 pass is back, which is done first now @@ -198,23 +203,95 @@ 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; +decim = 1; +AGC_coeffs.decimation = AGC_params.decimation; + +total_DC_gain = 0; for stage = 1:n_AGC_stages tau = AGC_params.time_constants(stage); + decim = decim * AGC_params.decimation(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 + + n_iterations = 1; % how many times to apply smoothing filter in a stage + % effective number of smoothings in a time constant: + ntimes = n_iterations * tau * (fs / decim); % 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); + % t is the variance of the distribution (impulse response width squared) + t1 = (AGC1_scales(stage)^2) / ntimes; + t2 = (AGC2_scales(stage)^2) / ntimes; + % the back-and-forth IIR method coefficients: + polez1 = 1 + 1/t1 - sqrt((1+1/t1)^2 - 1); + polez2 = 1 + 1/t2 - sqrt((1+1/t2)^2 - 1); + AGC_coeffs.AGC_polez1(stage) = polez1; + AGC_coeffs.AGC_polez2(stage) = polez2; + + % above method has spatial "delay" near sqrt(t2) - sqrt(t1) per time, + % or net delay that is not independent of ntimes. A problem??? + + % try a 3-tap FIR as an alternative: + n_taps = 3; + % from geometric distribution mean and variance from wikipedia: + delay = polez2/(1-polez2) - polez1/(1-polez1); + total_delay = delay * ntimes + spread_sq = polez1/(1-polez1)^2 + polez2/(1-polez2)^2; + total_spread = sqrt(spread_sq * ntimes) + + delay1 = delay; % keep this as 1-iteration delay reference... + a = (spread_sq + delay*delay - delay) / 2; + b = (spread_sq + delay*delay + delay) / 2; + AGC_spatial_FIR = [a, 1 - a - b, b]; % stored as 5 taps + done = AGC_spatial_FIR(2) > 0.1; % not OK if center tap is too low + % if 1 iteration is not good with 3 taps go to 5 taps, then more + % iterations if needed, and then fall back to double-exponential IIR: + while ~done % smoothing condition, middle value + if n_taps == 3 + % first time through, go wider but stick to 1 iteration + n_taps = 5; + n_iterations = 1; + else + % already at 5 taps, so just increase iterations + n_iterations = n_iterations + 1; % number of times to apply spatial + end + ntimes = n_iterations * tau * (fs / decim); % effective number of smoothings + % divide the spatial variance by effective number of smoothings: + % t is the variance of the distribution (impulse response width squared) + t1 = (AGC1_scales(stage)^2) / ntimes; + t2 = (AGC2_scales(stage)^2) / ntimes; + spread_sq = t1 + t2; + delay = delay1 / n_iterations; % maybe better than sqrt(t2) - sqrt(t1); + % 5-tap design duplicates the a and b coeffs; stores just 3 coeffs: + % a and b from their sum and diff as before: (sum \pm diff) / 2: + a = ((spread_sq + delay*delay)*2/5 - delay*2/3) / 2; + b = ((spread_sq + delay*delay)*2/5 + delay*2/3) / 2; + AGC_spatial_FIR = [a/2, 1 - a - b, b/2]; % implicit dup of a and b + done = AGC_spatial_FIR(2) > 0.1; + end + % store the resulting FIR design in coeffs: + AGC_coeffs.AGC_spatial_iterations(stage) = n_iterations; + AGC_coeffs.AGC_spatial_FIR(:,stage) = AGC_spatial_FIR; + AGC_coeffs.AGC_n_taps(stage) = n_taps; + + total_DC_gain = total_DC_gain + AGC_params.AGC_stage_gain^(stage-1); + + % TODO (dicklyon) -- is this what we want? + if stage == 1 + AGC_coeffs.AGC_mix_coeffs(stage) = 0; + else + AGC_coeffs.AGC_mix_coeffs(stage) = AGC_params.AGC_mix_coeff / ... + (tau * (fs / decim)); + end end -AGC_coeffs.AGC_gain = gain; +AGC_coeffs.AGC_spatial_FIR +AGC_coeffs.AGC_spatial_iterations +AGC_coeffs.AGC_n_taps +AGC_coeffs.AGC_polez1 +AGC_coeffs.AGC_polez2 + +AGC_coeffs.AGC_gain = total_DC_gain; + %% the IHC design coeffs: function IHC_coeffs = CARFAC_DesignIHC(IHC_params, fs) @@ -276,6 +353,7 @@ % 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 @@ -296,6 +374,9 @@ % 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 @@ -305,27 +386,33 @@ % 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] +% decimation: [8 2 2 2] +% AGC1_scales: [1 2 4 8] +% AGC2_scales: [1.5000 3 6 12] % detect_scale: 0.1500 -% AGC_mix_coeff: 0.2500 +% AGC_mix_coeff: 0.3500 % 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 = -% 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 +% 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_spatial_FIR: [3x4 double] +% AGC_n_taps: [3 5 5 5] +% AGC_mix_coeffs: [0 0.0317 0.0159 0.0079] +% AGC_gain: 15 +% detect_scale: 0.0664 % ans = +% just_hwr: 0 % lpf_coeff: 0.4327 % out1_rate: 0.0023 % in1_rate: 0.0023 @@ -336,7 +423,7 @@ % rest_cap: 0 % rest_cap1: 0.9635 % rest_cap2: 0.9269 -% saturation_output: 0.1587 +% saturation_output: 0.1507
--- a/trunk/matlab/bmm/carfac/CARFAC_FilterStep.m Mon Feb 27 21:50:20 2012 +0000 +++ b/trunk/matlab/bmm/carfac/CARFAC_FilterStep.m Thu Mar 01 19:49:24 2012 +0000 @@ -29,6 +29,10 @@ zA = state.zA_memory; zB = state.zB_memory + state.dzB_memory; % AGC interpolation r0 = filter_coeffs.r_coeffs; +v_offset = filter_coeffs.v_offset; +v2_corner = filter_coeffs.v2_corner; +v_damp_max = filter_coeffs.v_damp_max; + r = r0 - filter_coeffs.c_coeffs .* (zA + zB); % now reduce state by r and rotate with the fixed cos/sin coeffs: @@ -39,10 +43,9 @@ 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 +zA = (((state.z2_memory - z2) .* filter_coeffs.velocity_scale) + ... + v_offset) .^ 2; +zA = v_damp_max * zA ./ (v2_corner + zA); % make it more like an "essential" nonlinearity % Get outputs from inputs and new z2 values: zY = filter_coeffs.g_coeffs .* (inputs + filter_coeffs.h_coeffs .* z2);
--- a/trunk/matlab/bmm/carfac/CARFAC_IHCStep.m Mon Feb 27 21:50:20 2012 +0000 +++ b/trunk/matlab/bmm/carfac/CARFAC_IHCStep.m Thu Mar 01 19:49:24 2012 +0000 @@ -27,11 +27,11 @@ if just_hwr ihc_out = max(0, filters_out); - state.ihc_accum = state.ihc_accum + max(0, ihc_out); + state.ihc_accum = state.ihc_accum + ihc_out; else one_cap = coeffs.one_cap; - detect = CARFAC_Detect(filters_out/2); % detect with HWR or so + detect = CARFAC_Detect(filters_out); % detect with HWR or so if one_cap ihc_out = detect .* state.cap_voltage; @@ -56,7 +56,7 @@ state.lpf2_state = state.lpf2_state + coeffs.lpf_coeff * ... (state.lpf1_state - state.lpf2_state); - ihc_out = state.lpf2_state; + ihc_out = state.lpf2_state - coeffs.rest_output; - state.ihc_accum = state.ihc_accum + max(0, ihc_out - coeffs.rest_output); + state.ihc_accum = state.ihc_accum + max(0, ihc_out); end
--- a/trunk/matlab/bmm/carfac/CARFAC_Init.m Mon Feb 27 21:50:20 2012 +0000 +++ b/trunk/matlab/bmm/carfac/CARFAC_Init.m Thu Mar 01 19:49:24 2012 +0000 @@ -49,33 +49,37 @@ 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: +CF_struct.AGC_state(1).decim_phase = zeros(n_AGC_stages, 1); % ints % 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); - CF_struct.filter_state(mic).zA_memory = zeros(CF_struct.n_ch, 1); % cubic loop - CF_struct.filter_state(mic).zB_memory = zeros(CF_struct.n_ch, 1); % AGC interp - CF_struct.filter_state(mic).dzB_memory = zeros(CF_struct.n_ch, 1); % AGC incr - 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); + CF_struct.filter_state(mic).z1_memory = zeros(n_ch, 1); + CF_struct.filter_state(mic).z2_memory = zeros(n_ch, 1); + CF_struct.filter_state(mic).zA_memory = zeros(n_ch, 1); % cubic loop + CF_struct.filter_state(mic).zB_memory = zeros(n_ch, 1); % AGC interp + 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); % AGC loop filters' state: - CF_struct.AGC_state(mic).AGC_memory = zeros(CF_struct.n_ch, n_AGC_stages); % HACK init - CF_struct.AGC_state(mic).AGC_sum = zeros(CF_struct.n_ch, 1); + 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 % IHC state: if CF_struct.IHC_coeffs.just_hwr - CF_struct.IHC_state(mic).ihc_accum = zeros(CF_struct.n_ch, 1); + CF_struct.IHC_state(mic).ihc_accum = zeros(n_ch, 1); else CF_struct.IHC_state(mic).cap_voltage = ... - CF_struct.IHC_coeffs.rest_cap * ones(CF_struct.n_ch, 1); + CF_struct.IHC_coeffs.rest_cap * ones(n_ch, 1); CF_struct.IHC_state(mic).cap1_voltage = ... - CF_struct.IHC_coeffs.rest_cap1 * ones(CF_struct.n_ch, 1); + CF_struct.IHC_coeffs.rest_cap1 * ones(n_ch, 1); CF_struct.IHC_state(mic).cap2_voltage = ... - CF_struct.IHC_coeffs.rest_cap2 * ones(CF_struct.n_ch, 1); + CF_struct.IHC_coeffs.rest_cap2 * ones(n_ch, 1); CF_struct.IHC_state(mic).lpf1_state = ... - CF_struct.IHC_coeffs.rest_output * zeros(CF_struct.n_ch, 1); + CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1); CF_struct.IHC_state(mic).lpf2_state = ... - CF_struct.IHC_coeffs.rest_output * zeros(CF_struct.n_ch, 1); - CF_struct.IHC_state(mic).ihc_accum = zeros(CF_struct.n_ch, 1); + CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1); + CF_struct.IHC_state(mic).ihc_accum = zeros(n_ch, 1); end end
--- a/trunk/matlab/bmm/carfac/CARFAC_Run.m Mon Feb 27 21:50:20 2012 +0000 +++ b/trunk/matlab/bmm/carfac/CARFAC_Run.m Thu Mar 01 19:49:24 2012 +0000 @@ -1,5 +1,5 @@ % Copyright 2012, Google, Inc. -% Author: Richard F. Lyon +% 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" @@ -19,8 +19,8 @@ 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) +% function [naps, CF, decim_naps] = CARFAC_Run ... +% (CF, input_waves, 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). % @@ -58,91 +58,89 @@ 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; +% 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; +NAP_decim = 8; if nargout > 2 % make decimated detect output: - decim_naps = zeros(ceil(n_samp/decim), CF.n_ch, CF.n_mics); + decim_naps = zeros(ceil(n_samp/NAP_decim), CF.n_ch, CF.n_mics); else decim_naps = []; end -decim_k = 0; -sum_abs_response = 0; +k_AGC = 0; +AGC_plot_decim = 16; % how often to plot AGC state; TODO: use segments + +detects = zeros(n_ch, n_mics); for k = 1:n_samp - CF.k_mod_decim = mod(CF.k_mod_decim + 1, decim); % 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 [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); - + + detects(:, mic) = ihc_out; % for input to AGC, and out to SAI + 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 + if ~isempty(decim_naps) && (k_NAP_decim == 0) + decim_k = decim_k + 1; % index of decimated NAP + for mic = 1:n_mics + decim_naps(decim_k, :, mic) = CF.IHC_state(mic).ihc_accum / ... + NAP_decim; % for cochleagram + CF.IHC_state(mic).ihc_accum = zeros(n_ch,1); end - - % get the avg_detects to connect filter_state to AGC_state: - avg_detects = zeros(n_ch, n_mics); + end + % run the AGC update step, taking input from IHC_state, decimating + % internally, all mics at once due to mixing across them: + [CF.AGC_state, updated] = ... + CARFAC_AGCStep(CF.AGC_coeffs, detects, CF.AGC_state); + + % connect the feedback from AGC_state to filter_state when it updates + if updated 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); + 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) ... - / 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, max(maxsum))]); - drawnow + / decim1; end end + + k_AGC = mod(k_AGC + 1, AGC_plot_decim); + if AGC_plot_fig_num && k_AGC == 0 + figure(AGC_plot_fig_num); hold off; clf + set(gca, 'Position', [.25, .25, .5, .5]) + + maxsum = 0; + for mic = 1:n_mics + plot(CF.AGC_state(mic).AGC_memory(:, 1), 'k-', 'LineWidth', 1) + maxes(mic) = max(CF.AGC_state(mic).AGC_memory(:)); + hold on + stage1 = 4; % as opposed to stage + for stage = 1:3; + plot(2^(stage1-1) * (CF.AGC_state(mic).AGC_memory(:, stage) - ... + 2 * CF.AGC_state(mic).AGC_memory(:, stage+1))); + end + stage = 4; + plot(2^(stage1-1) * CF.AGC_state(mic).AGC_memory(:, stage)); + end + axis([0, CF.n_ch+1, -0.01, max(maxes) + 0.01]); + drawnow + end + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/trunk/matlab/bmm/carfac/CARFAC_Spatial_Smooth.m Thu Mar 01 19:49:24 2012 +0000 @@ -0,0 +1,35 @@ +function stage_state = CARFAC_Spatial_Smooth(coeffs, stage, stage_state) +% function AGC_state = CARFAC_Spatial_Smooth( ... +% n_taps, n_iterations, FIR_coeffs, AGC_state) + +n_iterations = coeffs.AGC_spatial_iterations(stage); + +use_FIR = n_iterations < 4; % or whatever condition we want to try + +if use_FIR + FIR_coeffs = coeffs.AGC_spatial_FIR(:,stage); + switch coeffs.AGC_n_taps(stage) + case 3 + for iter = 1:n_iterations + stage_state = ... + FIR_coeffs(1) * stage_state([1, 1:(end-1)], :) + ... + FIR_coeffs(2) * stage_state + ... + FIR_coeffs(3) * stage_state([2:end, end], :); + end + case 5 % 5-tap smoother duplicates first and last coeffs: + for iter = 1:n_iterations + stage_state = ... + FIR_coeffs(1) * (stage_state([1, 1, 1:(end-2)], :) + ... + stage_state([1, 1:(end-1)], :)) + ... + FIR_coeffs(2) * stage_state + ... + FIR_coeffs(3) * (stage_state([2:end, end], :) + ... + stage_state([3:end, end, end], :)); + end + otherwise + error('Bad n_taps in CARFAC_Spatial_Smooth'); + end +else + % use IIR method, back-and-forth firt-order smoothers: + stage_state = SmoothDoubleExponential(stage_state, ... + coeffs.AGC_polez1(stage), coeffs.AGC_polez2(stage)); +end
--- a/trunk/matlab/bmm/carfac/MultiScaleSmooth.m Mon Feb 27 21:50:20 2012 +0000 +++ b/trunk/matlab/bmm/carfac/MultiScaleSmooth.m Thu Mar 01 19:49:24 2012 +0000 @@ -51,10 +51,10 @@ end figure(scale_no + fig_offset1) - imagesc(smoothed') + imagesc(squeeze(smoothed(:,:,1))') figure(scale_no + fig_offset2) - plot(mean(smoothed, 2)); + plot(squeeze(mean(smoothed, 2))); drawnow pause(1)