flatmax@674: 
flatmax@598: // Author Matt Flax <flatmax@>
flatmax@597: //
flatmax@597: // This C++ file is part of an implementation of Lyon's cochlear model:
flatmax@597: // "Cascade of Asymmetric Resonators with Fast-Acting Compression"
flatmax@597: // to supplement Lyon's upcoming book "Human and Machine Hearing"
flatmax@597: //
flatmax@597: // Licensed under the Apache License, Version 2.0 (the "License");
flatmax@597: // you may not use this file except in compliance with the License.
flatmax@597: // You may obtain a copy of the License at
flatmax@597: //
flatmax@597: //     http://www.apache.org/licenses/LICENSE-2.0
flatmax@597: //
flatmax@597: // Unless required by applicable law or agreed to in writing, software
flatmax@597: // distributed under the License is distributed on an "AS IS" BASIS,
flatmax@597: // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
flatmax@597: // See the License for the specific language governing permissions and
flatmax@597: // limitations under the License.
flatmax@597: /**
flatmax@597:     \author {Matt Flax <flatmax\@>}
flatmax@597:     \date 2013.02.08
flatmax@597: */
flatmax@597: 
flatmax@597: #include "CAR.H"
flatmax@597: 
flatmax@597: CAR::CAR() {
flatmax@597:     //ctor
flatmax@597: }
flatmax@597: 
flatmax@597: CAR::~CAR() {
flatmax@597:     //dtor
flatmax@597: }
flatmax@598: 
flatmax@598: void CAR::designFilters(FP_TYPE fs, int n_ch) {
flatmax@598:     // don't really need these zero arrays, but it's a clue to what fields
flatmax@598:     // and types are need in ohter language implementations:
flatmax@600:     coeff.r1_coeffs=Matrix<FP_TYPE, Dynamic, 1>::Zero(n_ch,1); // resize and zero
flatmax@600:     coeff.a0_coeffs=Matrix<FP_TYPE, Dynamic, 1>::Zero(n_ch,1);
flatmax@600:     coeff.c0_coeffs=Matrix<FP_TYPE, Dynamic, 1>::Zero(n_ch,1);
flatmax@600:     coeff.h_coeffs=Matrix<FP_TYPE, Dynamic, 1>::Zero(n_ch,1);
flatmax@600:     coeff.g0_coeffs=Matrix<FP_TYPE, Dynamic, 1>::Zero(n_ch,1);
flatmax@598:     // zr_coeffs is not zeroed ... perhaps it should be ?
flatmax@598: 
flatmax@598: // zero_ratio comes in via h.  In book's circuit D, zero_ratio is 1/sqrt(a),
flatmax@598: // and that a is here 1 / (1+f) where h = f*c.
flatmax@598: // solve for f:  1/zero_ratio^2 = 1 / (1+f)
flatmax@598: // zero_ratio^2 = 1+f => f = zero_ratio^2 - 1
flatmax@598:     FP_TYPE f = pow(param.zero_ratio,2.) - 1.;  // nominally 1 for half-octave
flatmax@598: 
flatmax@598: // Make pole positions, s and c coeffs, h and g coeffs, etc.,
flatmax@598: // which mostly depend on the pole angle theta:
flatmax@598:     Array<FP_TYPE, Dynamic, 1> theta = pole_freqs * (2. * M_PI / fs);
flatmax@598: 
flatmax@598: 
flatmax@598: // undamped coupled-form coefficients:
flatmax@598:     coeff.c0_coeffs = theta.sin();
flatmax@598:     coeff.a0_coeffs = theta.cos();
flatmax@598: 
flatmax@598: // different possible interpretations for min-damping r:
flatmax@598: // r = exp(-theta * CF_CAR_params.min_zeta).
flatmax@598: // Compress theta to give somewhat higher Q at highest thetas:
flatmax@598:     FP_TYPE ff = param.high_f_damping_compression; // 0 to 1 typ. 0.5
flatmax@598:     Array<FP_TYPE, Dynamic,1> x = theta/M_PI;
flatmax@598: 
flatmax@598:     coeff.zr_coeffs = M_PI * (x - ff * x.pow(3.));  // when ff is 0, this is just theta,
flatmax@598: //                       and when ff is 1 it goes to zero at theta = pi.
flatmax@598:     coeff.r1_coeffs = (1. - coeff.zr_coeffs.array() * param.max_zeta);  // "r1" for the max-damping condition
flatmax@598: 
flatmax@598: // Increase the min damping where channels are spaced out more, by pulling
flatmax@598: // 25% of the way toward ERB_Hz/pole_freqs (close to 0.1 at high f)
flatmax@598:     Array<FP_TYPE, Dynamic, 1> min_zetas = param.min_zeta + 0.25*(PsychoAcoustics::Hz2ERB(pole_freqs, param.ERB_break_freq, param.ERB_Q).array() / pole_freqs - param.min_zeta);
flatmax@598:     coeff.zr_coeffs = coeff.zr_coeffs.array() * (param.max_zeta - min_zetas); // how r relates to undamping
flatmax@598: 
flatmax@598: // the zeros follow via the h_coeffs
flatmax@598:     coeff.h_coeffs = coeff.c0_coeffs * f;
flatmax@598: 
flatmax@598: // for unity gain at min damping, radius r; only used in CARFAC_Init:
flatmax@598:     Array<FP_TYPE, Dynamic,1> relative_undamping(n_ch, 1);
flatmax@600:     relative_undamping=Array<FP_TYPE, Dynamic, 1>::Zero(n_ch, 1).cos();
flatmax@598: 
flatmax@598: // this function needs to take CAR_coeffs even if we haven't finished
flatmax@598: // constucting it by putting in the g0_coeffs:
flatmax@598:     coeff.g0_coeffs = stageG(relative_undamping);
flatmax@598: }
flatmax@598: 
flatmax@598: Array<FP_TYPE, Dynamic, 1> CAR::stageG(Array<FP_TYPE, Dynamic, 1> &relative_undamping) {
flatmax@598:     // at max damping
flatmax@598:     Array<FP_TYPE, Dynamic, 1> r  = coeff.r1_coeffs.array() + coeff.zr_coeffs.array() * relative_undamping;
flatmax@598:     return (1. - 2.*r*coeff.a0_coeffs.array() + r.pow(2.)) / (1. - 2.*r*coeff.a0_coeffs.array() + coeff.h_coeffs.array()*r*coeff.c0_coeffs.array() + pow(r,2.));
flatmax@598: }