Mercurial > hg > aimc
changeset 473:b4da807f4318
Add CARFAC_Design_Doc.txt, CARFAC_Run_Segment.m, and some renames; rename various variables to be more parallel; clean up init code and such.
| author | dicklyon@google.com |
|---|---|
| date | Fri, 16 Mar 2012 04:19:24 +0000 |
| parents | 7ce064380002 |
| children | 3a873d04a7fe |
| files | matlab/bmm/carfac/CARFAC_AGCStep.m 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_FilterStep.m matlab/bmm/carfac/CARFAC_IHCStep.m matlab/bmm/carfac/CARFAC_IHC_Step.m matlab/bmm/carfac/CARFAC_Init.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_Spatial_Smooth.m matlab/bmm/carfac/CARFAC_Stage_g.m matlab/bmm/carfac/CARFAC_Transfer_Functions.m matlab/bmm/carfac/CARFAC_hacking.m |
| diffstat | 17 files changed, 716 insertions(+), 465 deletions(-) [+] |
line wrap: on
line diff
--- a/matlab/bmm/carfac/CARFAC_AGCStep.m Mon Mar 12 06:14:53 2012 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,115 +0,0 @@ -% 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, 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_mics = length(state); -[n_ch, n_AGC_stages] = size(state(1).AGC_memory); % number of channels - -optimize_for_mono = n_mics == 1; % mono optimization - -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 - -% 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 - - 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 - this_stage_sum = this_stage_sum + AGC_stage_state; - end - end - 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/matlab/bmm/carfac/CARFAC_AGC_Step.m Fri Mar 16 04:19:24 2012 +0000 @@ -0,0 +1,115 @@ +% 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, updated] = CARFAC_AGC_Step(AGC_coeffs, detects, state) +% function [state, updated] = CARFAC_AGC_Step(AGC_coeffs, detects, state) +% +% one time step (at decimated low AGC rate) of the AGC state update + +n_ears = length(state); +[n_ch, n_AGC_stages] = size(state(1).AGC_memory); % number of channels + +optimize_for_mono = n_ears == 1; % mono optimization + +stage = 1; +ins = AGC_coeffs.detect_scale * detects; +[state, updated] = CARFAC_AGC_Recurse(AGC_coeffs, ins, n_AGC_stages, ... + n_ears, n_ch, optimize_for_mono, stage, state); + + + + + +function [state, updated] = CARFAC_AGC_Recurse(coeffs, ins, n_stages, ... + n_ears, n_ch, mono, stage, state) +% function [state, updated = CARFAC_AGC_Recurse(coeffs, ins, n_stages, ... +% n_ears, 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 ear = 1:n_ears + state(ear).input_accum(:, stage) = ... + state(ear).input_accum(:, stage) + ins(:, ear); +end + +% 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 ear = 1:n_ears + ins(:,ear) = state(ear).input_accum(:, stage) / decim; + state(ear).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_ears, 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 ear = 1:n_ears + AGC_in = ins(:,ear); % the newly decimated input for this ear + +% 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(ear).AGC_memory(:, stage+1); + end + + AGC_stage_state = state(ear).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(ear).AGC_memory(:, stage) = AGC_stage_state; + + if ~mono + if ear == 1 + this_stage_sum = AGC_stage_state; + else + this_stage_sum = this_stage_sum + AGC_stage_state; + end + end + end + if ~mono + mix_coeff = coeffs.AGC_mix_coeffs(stage); + if mix_coeff > 0 + this_stage_mean = this_stage_sum / n_ears; + for ear = 1:n_ears + state(ear).AGC_memory(:, stage) = ... + state(ear).AGC_memory(:, stage) + ... + mix_coeff * ... + (this_stage_mean - state(ear).AGC_memory(:, stage)); + end + end + end + updated = 1; % bool to say we have new state +else + updated = 0; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/matlab/bmm/carfac/CARFAC_CAR_Step.m Fri Mar 16 04:19:24 2012 +0000 @@ -0,0 +1,77 @@ +% 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_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: +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_offset = CAR_coeffs.v_offset; +v2_corner = CAR_coeffs.v2_corner; +v_damp_max = CAR_coeffs.v_damp_max; + +% 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 - ... + CAR_coeffs.c0_coeffs .* state.z2_memory); +% z1 = z1 + inputs; +z2 = r .* (CAR_coeffs.c0_coeffs .* state.z1_memory + ... + CAR_coeffs.a0_coeffs .* state.z2_memory); + +% update the "velocity" for cubic nonlinearity, into zA: +zA = (((state.z2_memory - z2) .* CAR_coeffs.velocity_scale) + ... + v_offset) .^ 2; +% soft saturation to make it more like an "essential" nonlinearity: +zA = v_damp_max * zA ./ (v2_corner + zA); + +zY = CAR_coeffs.h_coeffs .* z2; % partial output + +% 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 final channel outputs in zY + in_out = g(ch) * (in_out + zY(ch)); + zY(ch) = in_out; +end + +% 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 +state.detect_accum = state.detect_accum + hwr_detect; +
--- a/matlab/bmm/carfac/CARFAC_Close_AGC_Loop.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Close_AGC_Loop.m Fri Mar 16 04:19:24 2012 +0000 @@ -23,13 +23,13 @@ % 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 +for ear = 1:CF.n_ears + extra_damping = CF.AGC_state(ear).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); + new_g = CARFAC_Stage_g(CF.CAR_coeffs, 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; + CF.CAR_state(ear).dzB_memory = ... + (extra_damping - CF.CAR_state(ear).zB_memory) / decim1; + CF.CAR_state(ear).dg_memory = ... + (new_g - CF.CAR_state(ear).g_memory) / decim1; end
--- a/matlab/bmm/carfac/CARFAC_Design.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Design.m Fri Mar 16 04:19:24 2012 +0000 @@ -17,9 +17,9 @@ % See the License for the specific language governing permissions and % limitations under the License. -function CF = CARFAC_Design(fs, CF_filter_params, ... +function CF = CARFAC_Design(fs, CF_CAR_params, ... CF_AGC_params, ERB_break_freq, ERB_Q, CF_IHC_params) -% function CF = CARFAC_Design(fs, CF_filter_params, ... +% function CF = CARFAC_Design(fs, CF_CAR_params, ... % CF_AGC_params, ERB_break_freq, ERB_Q, CF_IHC_params) % % This function designs the CARFAC (Cascade of Asymmetric Resonators with @@ -27,7 +27,7 @@ % 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_CAR_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 % @@ -91,7 +91,7 @@ end if nargin < 2 - CF_filter_params = struct( ... + CF_CAR_params = struct( ... '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 @@ -109,20 +109,20 @@ end % first figure out how many filter stages (PZFC/CARFAC channels): -pole_Hz = CF_filter_params.first_pole_theta * fs / (2*pi); +pole_Hz = CF_CAR_params.first_pole_theta * fs / (2*pi); n_ch = 0; -while pole_Hz > CF_filter_params.min_pole_Hz +while pole_Hz > CF_CAR_params.min_pole_Hz n_ch = n_ch + 1; - pole_Hz = pole_Hz - CF_filter_params.ERB_per_step * ... + pole_Hz = pole_Hz - CF_CAR_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); +pole_Hz = CF_CAR_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 * ... + pole_Hz = pole_Hz - CF_CAR_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 @@ -132,45 +132,48 @@ CF = struct( ... 'fs', fs, ... 'max_channels_per_octave', max_channels_per_octave, ... - 'filter_params', CF_filter_params, ... + 'CAR_params', CF_CAR_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 ); + 'CAR_coeffs', CARFAC_DesignFilters(CF_CAR_params, fs, pole_freqs), ... + 'AGC_coeffs', CARFAC_DesignAGC(CF_AGC_params, fs, n_ch), ... + 'IHC_coeffs', CARFAC_DesignIHC(CF_IHC_params, fs, n_ch), ... + 'n_ears', 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) +function CAR_coeffs = CARFAC_DesignFilters(CAR_params, fs, pole_freqs) n_ch = length(pole_freqs); % the filter design coeffs: -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 ... +CAR_coeffs = struct( ... + '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 ... ); -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.g0_coeffs = zeros(n_ch, 1); +CAR_coeffs.r1_coeffs = zeros(n_ch, 1); +CAR_coeffs.a0_coeffs = zeros(n_ch, 1); +CAR_coeffs.c0_coeffs = zeros(n_ch, 1); +CAR_coeffs.h_coeffs = zeros(n_ch, 1); +CAR_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. % 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 +f = CAR_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: @@ -180,47 +183,50 @@ a0 = cos(theta); % different possible interpretations for min-damping r: -% r = exp(-theta * CF_filter_params.min_zeta). +% r = exp(-theta * CF_CAR_params.min_zeta). % Compress theta to give somewhat higher Q at highest thetas: -ff = filter_params.high_f_damping_compression; % 0 to 1; typ. 0.5 +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. -filter_coeffs.zr_coeffs = zr_coeffs; % how r relates to zeta +CAR_coeffs.zr_coeffs = zr_coeffs; % how r relates to zeta -min_zeta = filter_params.min_zeta; +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) ./ pole_freqs - min_zeta); r1 = (1 - zr_coeffs .* min_zeta); % "1" for the min-damping condition -filter_coeffs.r1_coeffs = r1; +CAR_coeffs.r1_coeffs = r1; % undamped coupled-form coefficients: -filter_coeffs.a0_coeffs = a0; -filter_coeffs.c0_coeffs = c0; +CAR_coeffs.a0_coeffs = a0; +CAR_coeffs.c0_coeffs = c0; % the zeros follow via the h_coeffs h = c0 .* f; -filter_coeffs.h_coeffs = h; +CAR_coeffs.h_coeffs = h; % for unity gain at min damping, radius r; only used in CARFAC_Init: extra_damping = zeros(size(r1)); -% this function needs to take filter_coeffs even if we haven't finished +% this function needs to take CAR_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); +CAR_coeffs.g0_coeffs = CARFAC_Stage_g(CAR_coeffs, extra_damping); %% the AGC design coeffs: -function AGC_coeffs = CARFAC_DesignAGC(AGC_params, fs) +function AGC_coeffs = CARFAC_DesignAGC(AGC_params, fs, n_ch) -AGC_coeffs = struct('AGC_stage_gain', AGC_params.AGC_stage_gain); +n_AGC_stages = AGC_params.n_stages; +AGC_coeffs = struct( ... + 'n_ch', n_ch, ... + 'n_AGC_stages', n_AGC_stages, ... + '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 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 = 1; AGC_coeffs.decimation = AGC_params.decimation; @@ -327,14 +333,15 @@ %% the IHC design coeffs: -function IHC_coeffs = CARFAC_DesignIHC(IHC_params, fs) +function IHC_coeffs = CARFAC_DesignIHC(IHC_params, fs, n_ch) 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(... + IHC_coeffs = struct( ... + 'n_ch', n_ch, ... 'just_hwr', 0, ... 'lpf_coeff', 1 - exp(-1/(IHC_params.tau_lpf * fs)), ... 'out_rate', 1 / (IHC_params.tau_out * fs), ... @@ -342,6 +349,7 @@ 'one_cap', IHC_params.one_cap); else IHC_coeffs = struct(... + 'n_ch', n_ch, ... 'just_hwr', 0, ... 'lpf_coeff', 1 - exp(-1/(IHC_params.tau_lpf * fs)), ... 'out1_rate', 1 / (IHC_params.tau1_out * fs), ... @@ -361,9 +369,9 @@ '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); + IHC_in = 0; % the get the IHC output rest level + for k = 1:20000 + [IHC_out, IHC_state] = CARFAC_IHC_Step(IHC_in, IHC_coeffs, IHC_state); end IHC_coeffs.rest_output = IHC_out; @@ -373,8 +381,8 @@ 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); + for k = 1:20000 + [IHC_out, IHC_state] = CARFAC_IHC_Step(IHC_in, IHC_coeffs, IHC_state); prev_IHC_out = IHC_out; IHC_in = -IHC_in; end @@ -388,24 +396,24 @@ % % % CF = CARFAC_Design -% CF.filter_params +% CF.CAR_params % CF.AGC_params -% CF.filter_coeffs +% CF.CAR_coeffs % CF.AGC_coeffs % CF.IHC_coeffs % % CF = % fs: 22050 % max_channels_per_octave: 12.1873 -% filter_params: [1x1 struct] +% CAR_params: [1x1 struct] % AGC_params: [1x1 struct] % IHC_params: [1x1 struct] % n_ch: 66 % pole_freqs: [66x1 double] -% filter_coeffs: [1x1 struct] +% CAR_coeffs: [1x1 struct] % AGC_coeffs: [1x1 struct] % IHC_coeffs: [1x1 struct] -% n_mics: 0 +% n_ears: 0 % ans = % velocity_scale: 0.2000 % v_offset: 0.0100
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/matlab/bmm/carfac/CARFAC_Design_Doc.txt Fri Mar 16 04:19:24 2012 +0000 @@ -0,0 +1,142 @@ +CARFAC Design Doc +by "Richard F. Lyon" <dicklyon@google.com> + + +The CAR-FAC (cascade of asymmetric resonators with fast-acting +compression) is a cochlear model implemented as an efficient sound +processor, for mono, stereo, or multi-channel sound inputs. This file +describes aspects of the software design. + +The implementation will include equivalent Maltab, C++, Python, and +perhaps other versions. They should all be based on the same set of +classes (or structs) and functions, with roughly equivalent +functionality and names where possible. The examples here are Matlab, +in which the structs are typeless, but in other languages each will be +class or type. + +The top-level class is CARFAC. A CARFAC object knows how to design +its details from a modest set of parameters, and knows how to process +sound signals to produce "neural activity patterns" (NAPs). It +includes various sub-objects representing the parameters, designs, and +states of different parts of the CAR-FAC model. + +The three main subsystems of the CARFAC 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: + + CAR_params, CAR_coeffs, CAR_state + IHC_params, IHC_coeffs, IHC_state + AGC_params, AGC_coeffs, AGC_state + +These names can be used both for the classes, and slightly modified +for the member variables, arguments, temps, etc. + + -- (presently we have "filters" where I say "CAR"; I'll change it) -- + +The params are inputs that specify the system; the coeffs are things +like filter coefficients, things computed once and used at run time; +the state is whatever internal state is needed between running the +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: + + CF = CARFAC_Design(CAR_params, IHC_params, AGC_params) + + +On "channels" and "ears": + +Stages of the cascade are called "channels", and there are n_ch of +them (n_ch is determined from other parameters, but we could also make +a function to pick params to get a desired n_ch). Since we already +used "channels" for that, sound input channels (one for monaural, two +for binaural, or more) are called "ears", and there are n_ears of them. +The parameter n_ears is set when state is initialized, but +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. + + +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 +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 +will all fit and have good hit rate in a typical 32 KB L1 d-cache). +In Matlab we use doubles (8-byte), but in C++ we intend to use floats +(4-byte) and SSE (via Eigen) for higher performance. Alternate +implementations are OK. + + +Run-time strategy: + +To support real-time applications, sound can be processed in short +segments, producing segments of NAPs: + + [NAPs, CF] = CARFAC_Run_Segment(CF, input_waves) + +Here the plurals ("NAPs" and "input_waves") suggest multiple ears. +These can be vectors, with length 1 for mono, or Matlab +multi-D arrays with a singleton last dimension. + +The CF's states are updated such that there will be no glitch when a next +segment is processed. Segments are of arbitrary positive-integer +length, not necessarily equal. It is not inefficient to use a segment +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: + + [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) + +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 +boolean "updated" indicates whether the AGC actually did some work and +has a new output that needs to be used to "close the loop" and modify +the CAR state. When it's true, there's one more step, involving both AGC +and CAR states, so it's a function on the CARFAC: + + CF = CARFAC_Close_AGC_Loop(CF) + +In Matlab, these functions return a modified copy of the state or +CARFAC; that's why we have the multiple returned values. In languages +that do more than just call-by-value, a reference will be passed so +the state can be modified in place instead. + + +C++ Eigen strategy: + +The mapping from Matlab's "parallel" operations on vectors of values +into C++ code will be done by using Eigen +(http://eigen.tuxfamily.org/) arrays, which support similar +source-level parallel arithmetic, so we don't have to write loops over +channels. Eigen also allows fairly efficient compilation to SSE or +Arm/Neon instructions, which can do four float operations per cycle. +So we should be able to easily get to efficient code on various +platforms. + +If there are similar strategies available in other languages, we should +use them. Python may have a route to using Eigen +(http://eigen.tuxfamily.org/index.php?title=PythonInterop).
--- a/matlab/bmm/carfac/CARFAC_FilterStep.m Mon Mar 12 06:14:53 2012 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,77 +0,0 @@ -% 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. - -% Most of the update is parallel; finally we ripple inputs at the end. - -% Local nonlinearity zA and AGC feedback zB reduce pole radius: -zA = state.zA_memory; -zB = state.zB_memory + state.dzB_memory; % AGC interpolation -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 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.a0_coeffs .* state.z1_memory - ... - filter_coeffs.c0_coeffs .* state.z2_memory); -% z1 = z1 + inputs; -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) + ... - v_offset) .^ 2; -% soft saturation to make it more like an "essential" nonlinearity: -zA = v_damp_max * zA ./ (v2_corner + zA); - -zY = filter_coeffs.h_coeffs .* z2; % partial output - -% 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 final channel outputs in zY - in_out = g(ch) * (in_out + zY(ch)); - zY(ch) = in_out; -end - -% 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 -state.detect_accum = state.detect_accum + hwr_detect; -
--- a/matlab/bmm/carfac/CARFAC_IHCStep.m Mon Mar 12 06:14:53 2012 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,62 +0,0 @@ -% 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 + ihc_out; -else - one_cap = coeffs.one_cap; - - detect = CARFAC_Detect(filters_out); % 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 - coeffs.rest_output; - - state.ihc_accum = state.ihc_accum + max(0, ihc_out); -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/matlab/bmm/carfac/CARFAC_IHC_Step.m Fri Mar 16 04:19:24 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_IHC_Step(filters_out, coeffs, state); +% function [ihc_out, state] = CARFAC_IHC_Step(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 + ihc_out; +else + one_cap = coeffs.one_cap; + + detect = CARFAC_Detect(filters_out); % 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 - coeffs.rest_output; + + state.ihc_accum = state.ihc_accum + max(0, ihc_out); +end
--- a/matlab/bmm/carfac/CARFAC_Init.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Init.m Fri Mar 16 04:19:24 2012 +0000 @@ -17,10 +17,10 @@ % 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) +function CF_struct = CARFAC_Init(CF_struct, n_ears) +% function CF_struct = CARFAC_Init(CF_struct, n_ears) % -% Initialize state for n_mics channels (default 1). +% Initialize state for n_ears 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 @@ -28,60 +28,99 @@ % level of object. I like fewer structs and class types. if nargin < 2 - n_mics = 1; % monaural + n_ears = 1; % monaural end % % this is probably what I'd do in the C++ version: -% if CF_struct.n_mics ~= n_mics; +% if CF_struct.n_ears ~= n_ears; % % 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); +% % make a struct arrray, one element per ear channel, numbered: +% for k = 1:n_ears +% CF_struct.state(k) = struct('ear_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. +% For now I don't ever free anything if n_ears is reduced; +% so be sure to respect n_ears, 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_ears = n_ears; -CF_struct.n_mics = n_mics; -n_ch = CF_struct.n_ch; +% These inits grow the struct arrays as needed: +for ear = 1:n_ears + % for now there's only one coeffs, not one per ear + CF_struct.CAR_state(ear) = CAR_Init_State(CF_struct.CAR_coeffs); + CF_struct.IHC_state(ear) = IHC_Init_State(CF_struct.IHC_coeffs); + CF_struct.AGC_state(ear) = AGC_Init_State(CF_struct.AGC_coeffs); +end -% keep all the decimator phase info in mic 1 state only: -CF_struct.AGC_state(1).decim_phase = zeros(n_AGC_stages, 1); % ints +% for ear = 1:n_ears +% CF_struct.CAR_state(ear).z1_memory = zeros(n_ch, 1); +% CF_struct.CAR_state(ear).z2_memory = zeros(n_ch, 1); +% CF_struct.CAR_state(ear).zA_memory = zeros(n_ch, 1); % cubic loop +% CF_struct.CAR_state(ear).zB_memory = zeros(n_ch, 1); % AGC interp +% CF_struct.CAR_state(ear).dzB_memory = zeros(n_ch, 1); % AGC incr +% CF_struct.CAR_state(ear).zY_memory = zeros(n_ch, 1); +% CF_struct.CAR_state(ear).detect_accum = zeros(n_ch, 1); +% CF_struct.CAR_state(ear).g_memory = ... +% CF_struct.CAR_coeffs(ear).g0_coeffs; % initial g for min_zeta +% CF_struct.CAR_state(ear).dg_memory = zeros(n_ch, 1); % g interp +% % AGC loop filters' state: +% CF_struct.AGC_state(ear).AGC_memory = zeros(n_ch, n_AGC_stages); % HACK init +% CF_struct.AGC_state(ear).input_accum = zeros(n_ch, n_AGC_stages); % HACK init +% % IHC state: +% if CF_struct.IHC_coeffs.just_hwr +% CF_struct.IHC_state(ear).ihc_accum = zeros(n_ch, 1); +% else +% CF_struct.IHC_state(ear).cap_voltage = ... +% CF_struct.IHC_coeffs.rest_cap * ones(n_ch, 1); +% CF_struct.IHC_state(ear).cap1_voltage = ... +% CF_struct.IHC_coeffs.rest_cap1 * ones(n_ch, 1); +% CF_struct.IHC_state(ear).cap2_voltage = ... +% CF_struct.IHC_coeffs.rest_cap2 * ones(n_ch, 1); +% CF_struct.IHC_state(ear).lpf1_state = ... +% CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1); +% CF_struct.IHC_state(ear).lpf2_state = ... +% CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1); +% CF_struct.IHC_state(ear).ihc_accum = zeros(n_ch, 1); +% end +% end -% This zeroing grows the struct array as needed: -for mic = 1:n_mics - 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); - 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 - % IHC state: - if CF_struct.IHC_coeffs.just_hwr - 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(n_ch, 1); - CF_struct.IHC_state(mic).cap1_voltage = ... - CF_struct.IHC_coeffs.rest_cap1 * ones(n_ch, 1); - CF_struct.IHC_state(mic).cap2_voltage = ... - CF_struct.IHC_coeffs.rest_cap2 * ones(n_ch, 1); - CF_struct.IHC_state(mic).lpf1_state = ... - CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1); - CF_struct.IHC_state(mic).lpf2_state = ... - CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1); - CF_struct.IHC_state(mic).ihc_accum = zeros(n_ch, 1); - end -end + +function state = CAR_Init_State(coeffs) +n_ch = coeffs.n_ch; +state = struct( ... + 'z1_memory', zeros(n_ch, 1), ... + 'z2_memory', zeros(n_ch, 1), ... + 'zA_memory', zeros(n_ch, 1), ... + 'zB_memory', zeros(n_ch, 1), ... + 'dzB_memory', zeros(n_ch, 1), ... + 'zY_memory', zeros(n_ch, 1), ... + 'detect_accum', zeros(n_ch, 1), ... + 'g_memory', coeffs.g0_coeffs, ... + 'dg_memory', zeros(n_ch, 1) ... + ); + + +function state = AGC_Init_State(coeffs) +n_ch = coeffs.n_ch; +n_AGC_stages = coeffs.n_AGC_stages; +state = struct( ... + 'AGC_memory', zeros(n_ch, n_AGC_stages), ... + 'input_accum', zeros(n_ch, n_AGC_stages), ... + 'decim_phase', zeros(n_AGC_stages, 1) ... % integer decimator phase + ); + + +function state = IHC_Init_State(coeffs) +n_ch = coeffs.n_ch; +state = struct( ... + 'ihc_accum', zeros(n_ch, 1), ... + 'cap_voltage', coeffs.rest_cap * ones(n_ch, 1), ... + '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) ... + ); + +
--- a/matlab/bmm/carfac/CARFAC_Run.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Run.m Fri Mar 16 04:19:24 2012 +0000 @@ -17,9 +17,9 @@ % See the License for the specific language governing permissions and % limitations under the License. -function [naps, CF, decim_naps] = CARFAC_Run ... +function [CF, decim_naps, naps] = CARFAC_Run ... (CF, input_waves, AGC_plot_fig_num) -% function [naps, CF, decim_naps] = CARFAC_Run ... +% function [CF, decim_naps, 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). @@ -42,94 +42,82 @@ % 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_samp, n_ears] = size(input_waves); n_ch = CF.n_ch; if nargin < 3 AGC_plot_fig_num = 0; end -if n_mics ~= CF.n_mics +if n_ears ~= CF.n_ears error('bad number of input_waves channels passed to CARFAC_Run') end -naps = zeros(n_samp, n_ch, n_mics); -decim_k = 0; -k_NAP_decim = 0; -NAP_decim = 8; -if nargout > 2 + +naps = zeros(n_samp, n_ch, n_ears); + +seglen = 16; +n_segs = ceil(n_samp / seglen); + +if nargout > 1 % make decimated detect output: - decim_naps = zeros(ceil(n_samp/NAP_decim), CF.n_ch, CF.n_mics); + decim_naps = zeros(n_segs, CF.n_ch, CF.n_ears); else decim_naps = []; end +if nargout > 2 + % make decimated detect output: + naps = zeros(n_samp, CF.n_ch, CF.n_ears); +else + naps = []; +end -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, 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)); - - detects(:, mic) = ihc_out; % for input to AGC, and out to SAI - - naps(k, :, mic) = ihc_out; % output to neural activity pattern - +for seg_num = 1:n_segs + if seg_num == n_segs + % The last segement may be short of seglen, but do it anyway: + k_range = (seglen*(seg_num - 1) + 1):n_samp; + else + k_range = seglen*(seg_num - 1) + (1:seglen); 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); + % Process a segment to get a slice of decim_naps, and plot AGC state: + [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(k_range, :)); + + if ~isempty(naps) + for ear = 1:n_ears + % Accumulate segment naps to make full naps + naps(k_range, :, ear) = seg_naps(:, :, ear); end 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 - CF = CARFAC_Close_AGC_Loop(CF); + if ~isempty(decim_naps) + for ear = 1:n_ears + decim_naps(seg_num, :, ear) = CF.IHC_state(ear).ihc_accum / seglen; + CF.IHC_state(ear).ihc_accum = zeros(n_ch,1); + end end - k_AGC = mod(k_AGC + 1, AGC_plot_decim); - if AGC_plot_fig_num && k_AGC == 0 + if AGC_plot_fig_num 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(:)); + for ear = 1:n_ears + plot(CF.AGC_state(ear).AGC_memory(:, 1), 'k-', 'LineWidth', 1) + maxes(ear) = max(CF.AGC_state(ear).AGC_memory(:)); hold on for stage = 1:3; - plot(2^(stage-1) * (CF.AGC_state(mic).AGC_memory(:, stage) - ... - 2 * CF.AGC_state(mic).AGC_memory(:, stage+1))); + plot(2^(stage-1) * (CF.AGC_state(ear).AGC_memory(:, stage) - ... + 2 * CF.AGC_state(ear).AGC_memory(:, stage+1))); end stage = 4; - plot(2^(stage-1) * CF.AGC_state(mic).AGC_memory(:, stage)); + plot(2^(stage-1) * CF.AGC_state(ear).AGC_memory(:, stage)); end - axis([0, CF.n_ch+1, -0.01, max(maxes) + 0.01]); + axis([0, CF.n_ch+1, 0.0, max(maxes) + 0.01]); drawnow end - + end + +
--- a/matlab/bmm/carfac/CARFAC_Run_Linear.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Run_Linear.m Fri Mar 16 04:19:24 2012 +0000 @@ -25,36 +25,36 @@ % however, unlike CARFAC_Run, it forces it to be linear, and gives a % linear (not detected) output. -saved_v_damp_max = CF.filter_coeffs.v_damp_max; -CF.filter_coeffs.v_damp_max = 0.00; % make it linear for now +saved_v_damp_max = CF.CAR_coeffs.v_damp_max; +CF.CAR_coeffs.v_damp_max = 0.00; % make it linear for now -[n_samp, n_mics] = size(input_waves); +[n_samp, n_ears] = size(input_waves); n_ch = CF.n_ch; -if n_mics ~= CF.n_mics +if n_ears ~= CF.n_ears error('bad number of input_waves channels passed to CARFAC_Run') end -for mic = 1:CF.n_mics +for ear = 1:CF.n_ears % 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 + CF.CAR_state(ear).zB_memory(:) = extra_damping; % interpolator state + CF.CAR_state(ear).dzB_memory(:) = 0; % interpolator slope + CF.CAR_state(ear).g_memory = CARFAC_Stage_g( ... + CF.CAR_coeffs(ear), extra_damping); + CF.CAR_state(ear).dg_memory(:) = 0; % interpolator slope end -naps = zeros(n_samp, n_ch, n_mics); +naps = zeros(n_samp, n_ch, n_ears); for k = 1:n_samp % 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)); - naps(k, :, mic) = filters_out; % linear + for ear = 1:n_ears + [filters_out, CF.CAR_state(ear)] = CARFAC_FilterStep( ... + input_waves(k, ear), CF.CAR_coeffs, CF.CAR_state(ear)); + naps(k, :, ear) = filters_out; % linear end % skip IHC and AGC updates end -CF.filter_coeffs.v_damp_max = saved_v_damp_max; +CF.CAR_coeffs.v_damp_max = saved_v_damp_max;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/matlab/bmm/carfac/CARFAC_Run_Segment.m Fri Mar 16 04:19:24 2012 +0000 @@ -0,0 +1,81 @@ +% 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] = CARFAC_Run_Segment(CF, input_waves) +% function [naps, CF, decim_naps] = CARFAC_Run_Segment(CF, input_waves) +% +% This function runs the CARFAC; that is, filters a 1 or more channel +% sound input segment to make one or more neural activity patterns (naps); +% it can be called multiple times for successive segments of any length, +% as long as the returned CF with modified state is passed back in each +% time. +% +% 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_CAR_params, CF_AGC_params, n_ears) +% transfns = CARFAC_Transfer_Functions(CF, to_chans, from_chans) + +[n_samp, n_ears] = size(input_waves); + +if n_ears ~= CF.n_ears + error('bad number of input_waves channels passed to CARFAC_Run') +end + +n_ch = CF.n_ch; +naps = zeros(n_samp, n_ch, n_ears); % allocate space for result + +detects = zeros(n_ch, n_ears); +for k = 1:n_samp + % at each time step, possibly handle multiple channels + for ear = 1:n_ears + [car_out, CF.CAR_state(ear)] = CARFAC_CAR_Step( ... + input_waves(k, ear), CF.CAR_coeffs, CF.CAR_state(ear)); + + % update IHC state & output on every time step, too + [ihc_out, CF.IHC_state(ear)] = CARFAC_IHC_Step( ... + car_out, CF.IHC_coeffs, CF.IHC_state(ear)); + + detects(:, ear) = ihc_out; % for input to AGC, and out to SAI + naps(k, :, ear) = ihc_out; % output to neural activity pattern + end + % run the AGC update step, taking input from IHC_state, + % decimating internally, all ears at once due to mixing across them: + [CF.AGC_state, updated] = CARFAC_AGC_Step( ... + CF.AGC_coeffs, detects, CF.AGC_state); + + % connect the feedback from AGC_state to CAR_state when it updates + if updated + CF = CARFAC_Close_AGC_Loop(CF); + end +end +
--- a/matlab/bmm/carfac/CARFAC_Spatial_Smooth.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Spatial_Smooth.m Fri Mar 16 04:19:24 2012 +0000 @@ -38,11 +38,11 @@ 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)], :) + ... + FIR_coeffs(1) * (stage_state([1, 2, 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], :)); + stage_state([3:end, end, end-1], :)); end otherwise error('Bad n_taps in CARFAC_Spatial_Smooth');
--- a/matlab/bmm/carfac/CARFAC_Stage_g.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Stage_g.m Fri Mar 16 04:19:24 2012 +0000 @@ -17,15 +17,14 @@ % 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) +function g = CARFAC_Stage_g(CAR_coeffs, extra_damping) +% function g = CARFAC_Stage_g(CAR_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; +r1 = CAR_coeffs.r1_coeffs; % not zero damping, but min 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; % 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 Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_Transfer_Functions.m Fri Mar 16 04:19:24 2012 +0000 @@ -102,8 +102,8 @@ % Return transfer functions of all stages as rational functions. n_ch = CF.n_ch; -coeffs = CF.filter_coeffs; -min_zeta = CF.filter_params.min_zeta; +coeffs = CF.CAR_coeffs; +min_zeta = CF.CAR_params.min_zeta; a0 = coeffs.a0_coeffs; c0 = coeffs.c0_coeffs; @@ -111,8 +111,8 @@ % get r, adapted if we have state: r = coeffs.r1_coeffs; -if isfield(CF, 'filter_state') - state = CF.filter_state; +if isfield(CF, 'CAR_state') + state = CF.CAR_state; zB = state.zB_memory; % current extra damping r = r - zr .* zB; else
--- a/matlab/bmm/carfac/CARFAC_hacking.m Mon Mar 12 06:14:53 2012 +0000 +++ b/matlab/bmm/carfac/CARFAC_hacking.m Fri Mar 16 04:19:24 2012 +0000 @@ -43,29 +43,25 @@ agc_plot_fig_num = 6; -for n_mics = 1:2 - CF_struct = CARFAC_Init(CF_struct, n_mics); +for n_ears = 1:2 + CF_struct = CARFAC_Init(CF_struct, n_ears); - [nap, CF_struct, nap_decim] = CARFAC_Run(CF_struct, ... - test_signal, agc_plot_fig_num); + [CF_struct, nap_decim, nap] = CARFAC_Run(CF_struct, test_signal, ... + agc_plot_fig_num); % nap = deskew(nap); % deskew doesn't make much difference - if n_mics == 1 % because this hack doesn't work for binarual yet + if n_ears == 1 % because this hack doesn't work for binarual yet MultiScaleSmooth(nap_decim, 4); end -% 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 + % Display results for 1 or 2 ears: + for ear = 1:n_ears + smooth_nap = nap_decim(:, :, ear); + if n_ears == 1 mono_max = max(smooth_nap(:)); end - figure(3 + mic + n_mics) % Makes figures 5, ... + figure(3 + ear + n_ears) % Makes figures 5, ... image(63 * ((max(0, smooth_nap)/mono_max)' .^ 0.5)) title('smooth nap from nap decim') colormap(1 - gray); @@ -73,10 +69,8 @@ % Show resulting data, even though M-Lint complains: CF_struct - CF_struct.k_mod_decim - CF_struct.filter_state + CF_struct.CAR_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: