tom@455: % Copyright 2012, Google, Inc. tom@455: % Author: Richard F. Lyon tom@455: % tom@455: % This Matlab file is part of an implementation of Lyon's cochlear model: tom@455: % "Cascade of Asymmetric Resonators with Fast-Acting Compression" tom@455: % to supplement Lyon's upcoming book "Human and Machine Hearing" tom@455: % tom@455: % Licensed under the Apache License, Version 2.0 (the "License"); tom@455: % you may not use this file except in compliance with the License. tom@455: % You may obtain a copy of the License at tom@455: % tom@455: % http://www.apache.org/licenses/LICENSE-2.0 tom@455: % tom@455: % Unless required by applicable law or agreed to in writing, software tom@455: % distributed under the License is distributed on an "AS IS" BASIS, tom@455: % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. tom@455: % See the License for the specific language governing permissions and tom@455: % limitations under the License. tom@455: tom@455: function CF = CARFAC_Design(fs, CF_filter_params, ... tom@455: CF_AGC_params, ERB_break_freq, ERB_Q, CF_IHC_params) tom@455: % function CF = CARFAC_Design(fs, CF_filter_params, ... tom@455: % CF_AGC_params, ERB_break_freq, ERB_Q, CF_IHC_params) tom@455: % tom@455: % This function designs the CARFAC (Cascade of Asymmetric Resonators with tom@455: % Fast-Acting Compression); that is, it take bundles of parameters and tom@455: % computes all the filter coefficients needed to run it. tom@455: % tom@455: % fs is sample rate (per second) tom@455: % CF_filter_params bundles all the pole-zero filter cascade parameters tom@455: % CF_AGC_params bundles all the automatic gain control parameters tom@455: % CF_IHC_params bundles all the inner hair cell parameters tom@455: % tom@455: % See other functions for designing and characterizing the CARFAC: tom@455: % [naps, CF] = CARFAC_Run(CF, input_waves) tom@455: % transfns = CARFAC_Transfer_Functions(CF, to_channels, from_channels) tom@455: % tom@455: % Defaults to Glasberg & Moore's ERB curve: tom@455: % ERB_break_freq = 1000/4.37; % 228.833 tom@455: % ERB_Q = 1000/(24.7*4.37); % 9.2645 tom@455: % tom@455: % All args are defaultable; for sample/default args see the code; they tom@455: % make 96 channels at default fs = 22050, 114 channels at 44100. tom@455: tom@455: if nargin < 6 tom@455: % HACK: these constant control the defaults tom@455: one_cap = 0; % bool; 0 for new two-cap hack tom@455: just_hwr = 0; % book; 0 for normal/fancy IHC; 1 for HWR tom@455: if just_hwr tom@455: CF_IHC_params = struct('just_hwr', 1); % just a simple HWR tom@455: else tom@455: if one_cap tom@455: CF_IHC_params = struct( ... dicklyon@462: 'just_hwr', just_hwr, ... % not just a simple HWR tom@455: 'one_cap', one_cap, ... % bool; 0 for new two-cap hack tom@455: 'tau_lpf', 0.000080, ... % 80 microseconds smoothing twice tom@455: 'tau_out', 0.0005, ... % depletion tau is pretty fast tom@455: 'tau_in', 0.010 ); % recovery tau is slower tom@455: else tom@455: CF_IHC_params = struct( ... dicklyon@462: 'just_hwr', just_hwr, ... % not just a simple HWR tom@455: 'one_cap', one_cap, ... % bool; 0 for new two-cap hack tom@455: 'tau_lpf', 0.000080, ... % 80 microseconds smoothing twice tom@455: 'tau1_out', 0.020, ... % depletion tau is pretty fast tom@455: 'tau1_in', 0.020, ... % recovery tau is slower tom@455: 'tau2_out', 0.005, ... % depletion tau is pretty fast tom@455: 'tau2_in', 0.005 ); % recovery tau is slower tom@455: end tom@455: end tom@455: end tom@455: tom@455: if nargin < 5 tom@455: % Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138 tom@455: % ERB = 24.7 * (1 + 4.37 * CF_Hz / 1000); tom@455: ERB_Q = 1000/(24.7*4.37); % 9.2645 tom@455: if nargin < 4 tom@455: ERB_break_freq = 1000/4.37; % 228.833 tom@455: end tom@455: end tom@455: tom@455: if nargin < 3 tom@455: CF_AGC_params = struct( ... tom@455: 'n_stages', 4, ... tom@455: 'time_constants', [1, 4, 16, 64]*0.002, ... tom@455: 'AGC_stage_gain', 2, ... % gain from each stage to next slower stage dicklyon@462: 'decimation', [8, 2, 2, 2], ... % how often to update the AGC states dicklyon@463: 'AGC1_scales', [1, 2, 4, 6]*1, ... % in units of channels dicklyon@463: 'AGC2_scales', [1, 2, 4, 6]*1.5, ... % spread more toward base tom@455: 'detect_scale', 0.15, ... % the desired damping range dicklyon@462: 'AGC_mix_coeff', 0.5); tom@455: end tom@455: tom@455: if nargin < 2 tom@455: CF_filter_params = struct( ... dicklyon@462: 'velocity_scale', 0.2, ... % for the "cubic" velocity nonlinearity dicklyon@462: 'v_offset', 0.01, ... % offset gives a quadratic part dicklyon@462: 'v2_corner', 0.2, ... % corner for essential nonlin dicklyon@462: 'v_damp_max', 0.01, ... % damping delta damping from velocity nonlin tom@455: 'min_zeta', 0.12, ... dicklyon@467: 'first_pole_theta', 0.85*pi, ... dicklyon@467: 'zero_ratio', sqrt(2), ... % how far zero is above pole dicklyon@467: 'ERB_per_step', 0.5, ... % assume G&M's ERB formula dicklyon@467: 'min_pole_Hz', 30 ); tom@455: end tom@455: tom@455: if nargin < 1 tom@455: fs = 22050; tom@455: end tom@455: tom@455: % first figure out how many filter stages (PZFC/CARFAC channels): tom@455: pole_Hz = CF_filter_params.first_pole_theta * fs / (2*pi); tom@455: n_ch = 0; tom@455: while pole_Hz > CF_filter_params.min_pole_Hz tom@455: n_ch = n_ch + 1; tom@455: pole_Hz = pole_Hz - CF_filter_params.ERB_per_step * ... tom@455: ERB_Hz(pole_Hz, ERB_break_freq, ERB_Q); tom@455: end tom@455: % Now we have n_ch, the number of channels, so can make the array tom@455: % and compute all the frequencies again to put into it: tom@455: pole_freqs = zeros(n_ch, 1); tom@455: pole_Hz = CF_filter_params.first_pole_theta * fs / (2*pi); tom@455: for ch = 1:n_ch tom@455: pole_freqs(ch) = pole_Hz; tom@455: pole_Hz = pole_Hz - CF_filter_params.ERB_per_step * ... tom@455: ERB_Hz(pole_Hz, ERB_break_freq, ERB_Q); tom@455: end tom@455: % now we have n_ch, the number of channels, and pole_freqs array tom@455: dicklyon@467: max_channels_per_octave = log(2) / log(pole_freqs(1)/pole_freqs(2)); dicklyon@467: tom@455: CF = struct( ... tom@455: 'fs', fs, ... dicklyon@467: 'max_channels_per_octave', max_channels_per_octave, ... tom@455: 'filter_params', CF_filter_params, ... tom@455: 'AGC_params', CF_AGC_params, ... tom@455: 'IHC_params', CF_IHC_params, ... tom@455: 'n_ch', n_ch, ... tom@455: 'pole_freqs', pole_freqs, ... tom@455: 'filter_coeffs', CARFAC_DesignFilters(CF_filter_params, fs, pole_freqs), ... tom@455: 'AGC_coeffs', CARFAC_DesignAGC(CF_AGC_params, fs), ... tom@455: 'IHC_coeffs', CARFAC_DesignIHC(CF_IHC_params, fs), ... tom@455: 'n_mics', 0 ); tom@455: tom@455: % adjust the AGC_coeffs to account for IHC saturation level to get right tom@455: % damping change as specified in CF.AGC_params.detect_scale tom@455: CF.AGC_coeffs.detect_scale = CF.AGC_params.detect_scale / ... tom@455: (CF.IHC_coeffs.saturation_output * CF.AGC_coeffs.AGC_gain); tom@455: tom@455: %% Design the filter coeffs: tom@455: function filter_coeffs = CARFAC_DesignFilters(filter_params, fs, pole_freqs) tom@455: tom@455: n_ch = length(pole_freqs); tom@455: tom@455: % the filter design coeffs: tom@455: dicklyon@462: filter_coeffs = struct('velocity_scale', filter_params.velocity_scale, ... dicklyon@462: 'v_offset', filter_params.v_offset, ... dicklyon@462: 'v2_corner', filter_params.v2_corner, ... dicklyon@462: 'v_damp_max', filter_params.v_damp_max ... dicklyon@462: ); tom@455: tom@455: filter_coeffs.r_coeffs = zeros(n_ch, 1); tom@455: filter_coeffs.a_coeffs = zeros(n_ch, 1); tom@455: filter_coeffs.c_coeffs = zeros(n_ch, 1); tom@455: filter_coeffs.h_coeffs = zeros(n_ch, 1); tom@455: filter_coeffs.g_coeffs = zeros(n_ch, 1); tom@455: tom@455: % zero_ratio comes in via h. In book's circuit D, zero_ratio is 1/sqrt(a), tom@455: % and that a is here 1 / (1+f) where h = f*c. tom@455: % solve for f: 1/zero_ratio^2 = 1 / (1+f) tom@455: % zero_ratio^2 = 1+f => f = zero_ratio^2 - 1 tom@455: f = filter_params.zero_ratio^2 - 1; % nominally 1 for half-octave tom@455: tom@455: % Make pole positions, s and c coeffs, h and g coeffs, etc., tom@455: % which mostly depend on the pole angle theta: tom@455: theta = pole_freqs .* (2 * pi / fs); tom@455: tom@455: % different possible interpretations for min-damping r: tom@455: % r = exp(-theta * CF_filter_params.min_zeta). tom@455: % Using sin gives somewhat higher Q at highest thetas. dicklyon@467: ff = 5; % fudge factor for theta distortion; at least 1.0 dicklyon@467: r = (1 - ff*sin(theta/ff) * filter_params.min_zeta); tom@455: filter_coeffs.r_coeffs = r; tom@455: tom@455: % undamped coupled-form coefficients: tom@455: filter_coeffs.a_coeffs = cos(theta); tom@455: filter_coeffs.c_coeffs = sin(theta); tom@455: tom@455: % the zeros follow via the h_coeffs tom@455: h = sin(theta) .* f; tom@455: filter_coeffs.h_coeffs = h; tom@455: dicklyon@467: % % unity gain at min damping, radius r: dicklyon@467: g = (1 - 2*r.*cos(theta) + r.^2) ./ ... dicklyon@463: (1 - 2*r .* cos(theta) + h .* r .* sin(theta) + r.^2); dicklyon@467: % or assume r is 1, for the zero-damping gain g0: dicklyon@467: g0 = (2 - 2*cos(theta)) ./ ... dicklyon@467: (2 - 2 * cos(theta) + h .* sin(theta)); tom@455: dicklyon@467: filter_coeffs.g_coeffs = g0; dicklyon@467: % make coeffs that can correct g0 to make g based on (1 - r).^2: dicklyon@467: filter_coeffs.gr_coeffs = ((g ./ g0) - 1) ./ ((1 - r).^2); tom@455: tom@455: %% the AGC design coeffs: tom@455: function AGC_coeffs = CARFAC_DesignAGC(AGC_params, fs) tom@455: dicklyon@462: AGC_coeffs = struct('AGC_stage_gain', AGC_params.AGC_stage_gain); tom@455: tom@455: % AGC1 pass is smoothing from base toward apex; tom@455: % AGC2 pass is back, which is done first now tom@455: AGC1_scales = AGC_params.AGC1_scales; tom@455: AGC2_scales = AGC_params.AGC2_scales; tom@455: tom@455: n_AGC_stages = AGC_params.n_stages; tom@455: AGC_coeffs.AGC_epsilon = zeros(1, n_AGC_stages); % the 1/(tau*fs) roughly dicklyon@462: decim = 1; dicklyon@462: AGC_coeffs.decimation = AGC_params.decimation; dicklyon@462: dicklyon@462: total_DC_gain = 0; tom@455: for stage = 1:n_AGC_stages dicklyon@464: tau = AGC_params.time_constants(stage); % time constant in seconds dicklyon@464: decim = decim * AGC_params.decimation(stage); % net decim to this stage tom@455: % epsilon is how much new input to take at each update step: tom@455: AGC_coeffs.AGC_epsilon(stage) = 1 - exp(-decim / (tau * fs)); dicklyon@462: % effective number of smoothings in a time constant: dicklyon@464: ntimes = tau * (fs / decim); % typically 5 to 50 dicklyon@463: dicklyon@463: % decide on target spread (variance) and delay (mean) of impulse dicklyon@463: % response as a distribution to be convolved ntimes: dicklyon@464: % TODO (dicklyon): specify spread and delay instead of scales??? dicklyon@463: delay = (AGC2_scales(stage) - AGC1_scales(stage)) / ntimes; dicklyon@463: spread_sq = (AGC1_scales(stage)^2 + AGC2_scales(stage)^2) / ntimes; dicklyon@463: dicklyon@464: % get pole positions to better match intended spread and delay of dicklyon@464: % [[geometric distribution]] in each direction (see wikipedia) dicklyon@463: u = 1 + 1 / spread_sq; % these are based on off-line algebra hacking. dicklyon@463: p = u - sqrt(u^2 - 1); % pole that would give spread if used twice. dicklyon@463: dp = delay * (1 - 2*p +p^2)/2; dicklyon@463: polez1 = p - dp; dicklyon@463: polez2 = p + dp; dicklyon@462: AGC_coeffs.AGC_polez1(stage) = polez1; dicklyon@462: AGC_coeffs.AGC_polez2(stage) = polez2; dicklyon@462: dicklyon@464: % try a 3- or 5-tap FIR as an alternative to the double exponential: dicklyon@464: n_taps = 0; dicklyon@464: FIR_OK = 0; dicklyon@464: n_iterations = 1; dicklyon@464: while ~FIR_OK dicklyon@464: switch n_taps dicklyon@464: case 0 dicklyon@464: % first attempt a 3-point FIR to apply once: dicklyon@464: n_taps = 3; dicklyon@464: case 3 dicklyon@464: % second time through, go wider but stick to 1 iteration dicklyon@464: n_taps = 5; dicklyon@464: case 5 dicklyon@464: % apply FIR multiple times instead of going wider: dicklyon@464: n_iterations = n_iterations + 1; dicklyon@464: if n_iterations > 16 dicklyon@464: error('Too many n_iterations in CARFAC_DesignAGC'); dicklyon@464: end dicklyon@464: otherwise dicklyon@464: % to do other n_taps would need changes in CARFAC_Spatial_Smooth dicklyon@464: % and in Design_FIR_coeffs dicklyon@464: error('Bad n_taps in CARFAC_DesignAGC'); dicklyon@462: end dicklyon@464: [AGC_spatial_FIR, FIR_OK] = Design_FIR_coeffs( ... dicklyon@464: n_taps, spread_sq, delay, n_iterations); dicklyon@462: end dicklyon@464: % when FIR_OK, store the resulting FIR design in coeffs: dicklyon@462: AGC_coeffs.AGC_spatial_iterations(stage) = n_iterations; dicklyon@462: AGC_coeffs.AGC_spatial_FIR(:,stage) = AGC_spatial_FIR; dicklyon@462: AGC_coeffs.AGC_n_taps(stage) = n_taps; dicklyon@462: dicklyon@464: % accumulate DC gains from all the stages, accounting for stage_gain: dicklyon@462: total_DC_gain = total_DC_gain + AGC_params.AGC_stage_gain^(stage-1); dicklyon@462: dicklyon@464: % TODO (dicklyon) -- is this the best binaural mixing plan? dicklyon@462: if stage == 1 dicklyon@462: AGC_coeffs.AGC_mix_coeffs(stage) = 0; dicklyon@462: else dicklyon@462: AGC_coeffs.AGC_mix_coeffs(stage) = AGC_params.AGC_mix_coeff / ... dicklyon@462: (tau * (fs / decim)); dicklyon@462: end tom@455: end tom@455: dicklyon@463: AGC_coeffs.AGC_gain = total_DC_gain; dicklyon@462: dicklyon@464: % % print some results dicklyon@464: % AGC_coeffs dicklyon@464: % AGC_spatial_FIR = AGC_coeffs.AGC_spatial_FIR dicklyon@464: dicklyon@464: dicklyon@464: %% dicklyon@464: function [FIR, OK] = Design_FIR_coeffs(n_taps, var, mn, n_iter) dicklyon@464: % function [FIR, OK] = Design_FIR_coeffs(n_taps, spread_sq, delay, n_iter) dicklyon@464: dicklyon@464: % reduce mean and variance of smoothing distribution by n_iterations: dicklyon@464: mn = mn / n_iter; dicklyon@464: var = var / n_iter; dicklyon@464: switch n_taps dicklyon@464: case 3 dicklyon@464: % based on solving to match mean and variance of [a, 1-a-b, b]: dicklyon@464: a = (var + mn*mn - mn) / 2; dicklyon@464: b = (var + mn*mn + mn) / 2; dicklyon@464: FIR = [a, 1 - a - b, b]; dicklyon@464: OK = FIR(2) >= 0.2; dicklyon@464: case 5 dicklyon@464: % based on solving to match [a/2, a/2, 1-a-b, b/2, b/2]: dicklyon@464: a = ((var + mn*mn)*2/5 - mn*2/3) / 2; dicklyon@464: b = ((var + mn*mn)*2/5 + mn*2/3) / 2; dicklyon@464: % first and last coeffs are implicitly duplicated to make 5-point FIR: dicklyon@464: FIR = [a/2, 1 - a - b, b/2]; dicklyon@464: OK = FIR(2) >= 0.1; dicklyon@464: otherwise dicklyon@464: error('Bad n_taps in AGC_spatial_FIR'); dicklyon@464: end dicklyon@462: tom@455: tom@455: %% the IHC design coeffs: tom@455: function IHC_coeffs = CARFAC_DesignIHC(IHC_params, fs) tom@455: tom@455: if IHC_params.just_hwr tom@455: IHC_coeffs = struct('just_hwr', 1); tom@455: IHC_coeffs.saturation_output = 10; % HACK: assume some max out tom@455: else tom@455: if IHC_params.one_cap tom@455: IHC_coeffs = struct(... tom@455: 'just_hwr', 0, ... tom@455: 'lpf_coeff', 1 - exp(-1/(IHC_params.tau_lpf * fs)), ... tom@455: 'out_rate', 1 / (IHC_params.tau_out * fs), ... tom@455: 'in_rate', 1 / (IHC_params.tau_in * fs), ... tom@455: 'one_cap', IHC_params.one_cap); tom@455: else tom@455: IHC_coeffs = struct(... tom@455: 'just_hwr', 0, ... tom@455: 'lpf_coeff', 1 - exp(-1/(IHC_params.tau_lpf * fs)), ... tom@455: 'out1_rate', 1 / (IHC_params.tau1_out * fs), ... tom@455: 'in1_rate', 1 / (IHC_params.tau1_in * fs), ... tom@455: 'out2_rate', 1 / (IHC_params.tau2_out * fs), ... tom@455: 'in2_rate', 1 / (IHC_params.tau2_in * fs), ... tom@455: 'one_cap', IHC_params.one_cap); tom@455: end tom@455: tom@455: % run one channel to convergence to get rest state: tom@455: IHC_coeffs.rest_output = 0; tom@455: IHC_state = struct( ... tom@455: 'cap_voltage', 0, ... tom@455: 'cap1_voltage', 0, ... tom@455: 'cap2_voltage', 0, ... tom@455: 'lpf1_state', 0, ... tom@455: 'lpf2_state', 0, ... tom@455: 'ihc_accum', 0); tom@455: tom@455: IHC_in = 0; tom@455: for k = 1:30000 tom@455: [IHC_out, IHC_state] = CARFAC_IHCStep(IHC_in, IHC_coeffs, IHC_state); tom@455: end tom@455: tom@455: IHC_coeffs.rest_output = IHC_out; tom@455: IHC_coeffs.rest_cap = IHC_state.cap_voltage; tom@455: IHC_coeffs.rest_cap1 = IHC_state.cap1_voltage; tom@455: IHC_coeffs.rest_cap2 = IHC_state.cap2_voltage; tom@455: tom@455: LARGE = 2; tom@455: IHC_in = LARGE; % "Large" saturating input to IHC; make it alternate tom@455: for k = 1:30000 tom@455: [IHC_out, IHC_state] = CARFAC_IHCStep(IHC_in, IHC_coeffs, IHC_state); tom@455: prev_IHC_out = IHC_out; tom@455: IHC_in = -IHC_in; tom@455: end tom@455: tom@455: IHC_coeffs.saturation_output = (IHC_out + prev_IHC_out) / 2; tom@455: end tom@455: tom@455: %% tom@455: % default design result, running this function with no args, should look tom@455: % like this, before CARFAC_Init puts state storage into it: tom@455: % dicklyon@462: % tom@455: % CF = CARFAC_Design tom@455: % CF.filter_params tom@455: % CF.AGC_params tom@455: % CF.filter_coeffs tom@455: % CF.AGC_coeffs tom@455: % CF.IHC_coeffs tom@455: % tom@455: % CF = tom@455: % fs: 22050 tom@455: % filter_params: [1x1 struct] tom@455: % AGC_params: [1x1 struct] tom@455: % IHC_params: [1x1 struct] tom@455: % n_ch: 96 tom@455: % pole_freqs: [96x1 double] tom@455: % filter_coeffs: [1x1 struct] tom@455: % AGC_coeffs: [1x1 struct] tom@455: % IHC_coeffs: [1x1 struct] tom@455: % n_mics: 0 tom@455: % ans = tom@455: % velocity_scale: 0.2000 dicklyon@462: % v_offset: 0.0100 dicklyon@462: % v2_corner: 0.2000 dicklyon@462: % v_damp_max: 0.0100 tom@455: % min_zeta: 0.1200 tom@455: % first_pole_theta: 2.4504 tom@455: % zero_ratio: 1.4142 tom@455: % ERB_per_step: 0.3333 tom@455: % min_pole_Hz: 40 tom@455: % ans = tom@455: % n_stages: 4 tom@455: % time_constants: [0.0020 0.0080 0.0320 0.1280] tom@455: % AGC_stage_gain: 2 dicklyon@462: % decimation: [8 2 2 2] dicklyon@462: % AGC1_scales: [1 2 4 8] dicklyon@462: % AGC2_scales: [1.5000 3 6 12] tom@455: % detect_scale: 0.1500 dicklyon@462: % AGC_mix_coeff: 0.3500 tom@455: % ans = tom@455: % velocity_scale: 0.2000 dicklyon@462: % v_offset: 0.0100 dicklyon@462: % v2_corner: 0.2000 dicklyon@462: % v_damp_max: 0.0100 tom@455: % r_coeffs: [96x1 double] tom@455: % a_coeffs: [96x1 double] tom@455: % c_coeffs: [96x1 double] tom@455: % h_coeffs: [96x1 double] tom@455: % g_coeffs: [96x1 double] tom@455: % ans = dicklyon@462: % AGC_stage_gain: 2 dicklyon@462: % AGC_epsilon: [0.1659 0.0867 0.0443 0.0224] dicklyon@462: % decimation: [8 2 2 2] dicklyon@462: % AGC_spatial_iterations: [1 1 2 3] dicklyon@462: % AGC_spatial_FIR: [3x4 double] dicklyon@462: % AGC_n_taps: [3 5 5 5] dicklyon@462: % AGC_mix_coeffs: [0 0.0317 0.0159 0.0079] dicklyon@462: % AGC_gain: 15 dicklyon@462: % detect_scale: 0.0664 tom@455: % ans = dicklyon@462: % just_hwr: 0 tom@455: % lpf_coeff: 0.4327 tom@455: % out1_rate: 0.0023 tom@455: % in1_rate: 0.0023 tom@455: % out2_rate: 0.0091 tom@455: % in2_rate: 0.0091 tom@455: % one_cap: 0 tom@455: % rest_output: 0.0365 tom@455: % rest_cap: 0 tom@455: % rest_cap1: 0.9635 tom@455: % rest_cap2: 0.9269 dicklyon@462: % saturation_output: 0.1507 tom@455: tom@455: