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( ... tom@455: 'just_hwr', 0, ... % 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( ... tom@455: 'just_hwr', 0, ... % 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 tom@455: 'decimation', 16, ... % how often to update the AGC states tom@455: 'AGC1_scales', [1, 2, 3, 4]*1, ... % in units of channels tom@455: 'AGC2_scales', [1, 2, 3, 4]*1.25, ... % spread more toward base tom@455: 'detect_scale', 0.15, ... % the desired damping range tom@455: 'AGC_mix_coeff', 0.25); tom@455: end tom@455: tom@455: if nargin < 2 tom@455: CF_filter_params = struct( ... tom@455: 'velocity_scale', 0.2, ... % for the cubic nonlinearity tom@455: 'min_zeta', 0.12, ... tom@455: 'first_pole_theta', 0.78*pi, ... tom@455: 'zero_ratio', sqrt(2), ... tom@455: 'ERB_per_step', 0.3333, ... % assume G&M's ERB formula tom@455: 'min_pole_Hz', 40 ); 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: tom@455: CF = struct( ... tom@455: 'fs', fs, ... 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: tom@455: filter_coeffs = struct('velocity_scale', filter_params.velocity_scale); 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. tom@455: r = (1 - sin(theta) * 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: tom@455: r2 = r; % aim for unity DC gain at min damping, here; or could try r^2 tom@455: filter_coeffs.g_coeffs = 1 ./ (1 + h .* r2 .* sin(theta) ./ ... tom@455: (1 - 2 * r2 .* cos(theta) + r2 .^ 2)); tom@455: tom@455: tom@455: %% the AGC design coeffs: tom@455: function AGC_coeffs = CARFAC_DesignAGC(AGC_params, fs) tom@455: tom@455: AGC_coeffs = struct('AGC_stage_gain', AGC_params.AGC_stage_gain, ... tom@455: 'AGC_mix_coeff', AGC_params.AGC_mix_coeff); tom@455: 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 tom@455: decim = AGC_params.decimation; tom@455: gain = 0; tom@455: for stage = 1:n_AGC_stages tom@455: tau = AGC_params.time_constants(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)); tom@455: % and these are the smoothing scales and poles for decimated rate: tom@455: ntimes = tau * (fs / decim); % effective number of smoothings tom@455: % divide the spatial variance by effective number of smoothings: tom@455: t = (AGC1_scales(stage)^2) / ntimes; % adjust scale for diffusion tom@455: AGC_coeffs.AGC1_polez(stage) = 1 + 1/t - sqrt((1+1/t)^2 - 1); tom@455: t = (AGC2_scales(stage)^2) / ntimes; % adjust scale for diffusion tom@455: AGC_coeffs.AGC2_polez(stage) = 1 + 1/t - sqrt((1+1/t)^2 - 1); tom@455: gain = gain + AGC_params.AGC_stage_gain^(stage-1); tom@455: end tom@455: tom@455: AGC_coeffs.AGC_gain = gain; 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: % 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 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 tom@455: % decimation: 16 tom@455: % AGC1_scales: [1 2 3 4] tom@455: % AGC2_scales: [1.2500 2.5000 3.7500 5] tom@455: % detect_scale: 0.1500 tom@455: % AGC_mix_coeff: 0.2500 tom@455: % ans = tom@455: % velocity_scale: 0.2000 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 = tom@455: % AGC_stage_gain: 2 tom@455: % AGC_mix_coeff: 0.2500 tom@455: % AGC_epsilon: [0.3043 0.0867 0.0224 0.0057] tom@455: % AGC1_polez: [0.1356 0.1356 0.0854 0.0417] tom@455: % AGC2_polez: [0.1872 0.1872 0.1227 0.0623] tom@455: % AGC_gain: 15 tom@455: % detect_scale: 0.0630 tom@455: % ans = 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 tom@455: % saturation_output: 0.1587 tom@455: tom@455: tom@455: