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 "AGC.H"
flatmax@597: 
flatmax@612: AGC::AGC() {
flatmax@597: }
flatmax@597: 
flatmax@612: AGC::~AGC() {
flatmax@597: }
flatmax@612: 
flatmax@612: void AGC::designAGC(FP_TYPE fs, int n_ch) {
flatmax@612:     int n_AGC_stages = params.n_stages;
flatmax@612: //AGC_coeffs = struct( ...
flatmax@612: //  'n_ch', n_ch, ...
flatmax@612: //  'n_AGC_stages', n_AGC_stages, ...
flatmax@612: //  'AGC_stage_gain', AGC_params.AGC_stage_gain);
flatmax@612: 
flatmax@612: // AGC1 pass is smoothing from base toward apex;
flatmax@612: // AGC2 pass is back, which is done first now (in double exp. version)
flatmax@612: //AGC1_scales = AGC_params.AGC1_scales;
flatmax@612: //AGC2_scales = AGC_params.AGC2_scales;
flatmax@612: 
flatmax@612:     coeffs.AGC_epsilon = Array<FP_TYPE, 1, Dynamic>::Zero(1, n_AGC_stages);  // the 1/(tau*fs) roughly
flatmax@612:     FP_TYPE decim = 1.;
flatmax@612: //AGC_coeffs.decimation = AGC_params.decimation;
flatmax@612: 
flatmax@612:     FP_TYPE total_DC_gain = 0.;
flatmax@612:     for (int stage = 1; stage<=n_AGC_stages; stage++) {
flatmax@612:         FP_TYPE tau = params.time_constants(stage-1); // time constant in seconds
flatmax@612:         decim = decim * params.decimation(stage-1); // net decim to this stage
flatmax@612:         // epsilon is how much new input to take at each update step:
flatmax@612:         coeffs.AGC_epsilon(stage-1) = 1. - exp(-decim / (tau * fs));
flatmax@612:         // effective number of smoothings in a time constant:
flatmax@612:         FP_TYPE ntimes = tau * (fs / decim);  // typically 5 to 50
flatmax@612: 
flatmax@612:         // decide on target spread (variance) and delay (mean) of impulse
flatmax@612:         // response as a distribution to be convolved ntimes:
flatmax@612:         // TODO (dicklyon): specify spread and delay instead of scales???
flatmax@612:         FP_TYPE delay = (param.AGC2_scales(stage-1) - param.AGC1_scales(stage-1)) / ntimes;
flatmax@612:         FP_TYPE spread_sq = (param.AGC1_scales(stage-1).pow(2.) + param.AGC2_scales(stage-1).pow(2)) / ntimes;
flatmax@612: 
flatmax@612:         // get pole positions to better match intended spread and delay of
flatmax@612:         // [[geometric distribution]] in each direction (see wikipedia)
flatmax@612:         FP_TYPE u = 1. + 1. / spread_sq; // these are based on off-line algebra hacking.
flatmax@612:         FP_TYPE p = u - sqrt(pow(u,2.) - 1.); // pole that would give spread if used twice.
flatmax@612:         FP_TYPE dp = delay * (1. - 2.*p +pow(p,2.))/2.;
flatmax@612:                               FP_TYPE polez1 = p - dp;
flatmax@612:                               FP_TYPE polez2 = p + dp;
flatmax@612:                               coeffs.AGC_polez1(stage) = polez1;
flatmax@612:                               coeffs.AGC_polez2(stage) = polez2;
flatmax@612: 
flatmax@612:                               // try a 3- or 5-tap FIR as an alternative to the double exponential:
flatmax@612:                               Array<FP_TYPE,1, Dynamic> AGC_spatial_FIR;
flatmax@612:                               int n_taps = 0;
flatmax@612:                               int FIR_OK = 0;
flatmax@612:                               int n_iterations = 1;
flatmax@612:         while (~FIR_OK) {
flatmax@612:         switch (n_taps) {
flatmax@612:             case 0:
flatmax@612:                 // first attempt a 3-point FIR to apply once:
flatmax@612:                 n_taps = 3;
flatmax@612:                 break;
flatmax@612:             case 3:
flatmax@612:                 // second time through, go wider but stick to 1 iteration
flatmax@612:                 n_taps = 5;
flatmax@612:                 break;
flatmax@612:             case 5:
flatmax@612:                 // apply FIR multiple times instead of going wider:
flatmax@612:                 n_iterations = n_iterations + 1;
flatmax@612:                 if (n_iterations > 16) {
flatmax@612:                     cerr<<"Too many n_iterations in CARFAC_DesignAGC"<<endl;
flatmax@612:                     exit(AGC_DESIGN_ITERATION_ERROR);
flatmax@612:                 }
flatmax@612:                 break;
flatmax@612:             default:
flatmax@612:                 // to do other n_taps would need changes in CARFAC_Spatial_Smooth
flatmax@612:                 // and in Design_FIR_coeffs
flatmax@612:                 cerr<<"Bad n_taps in CARFAC_DesignAGC"<<endl;
flatmax@612:                 exit(AGC_DESIGN_TAPS_OOB_ERROR);
flatmax@612:                 break;
flatmax@612:             }
flatmax@612:             FIR_OK = Design_FIR_coeffs(n_taps, spread_sq, delay, n_iterations,AGC_spatial_FIR);
flatmax@612:         }
flatmax@612:         // when FIR_OK, store the resulting FIR design in coeffs:
flatmax@612:         coeff.AGC_spatial_iterations(stage-1) = n_iterations;
flatmax@612:         coeff.AGC_spatial_FIR.col(stage-1).block(0,AGC_spatial_FIR.size()) = AGC_spatial_FIR;
flatmax@612:         coeff.AGC_spatial_n_taps(stage-1) = n_taps;
flatmax@612: 
flatmax@612:         // accumulate DC gains from all the stages, accounting for stage_gain:
flatmax@612:         total_DC_gain = total_DC_gain + params.AGC_stage_gain.pow(stage-1);
flatmax@612: 
flatmax@612:         // TODO (dicklyon) -- is this the best binaural mixing plan?
flatmax@612:         if (stage == 1)
flatmax@612:         coeff.AGC_mix_coeffs(stage-1) = 0.;
flatmax@612:         else
flatmax@612:             coeff.AGC_mix_coeffs(stage-1) = param.AGC_mix_coeff / (tau * (fs / decim));
flatmax@612:         }
flatmax@612: 
flatmax@612: coeff.AGC_gain = total_DC_gain;
flatmax@612: 
flatmax@612: // adjust the detect_scale to be the reciprocal DC gain of the AGC filters:
flatmax@612: AGC_coeffs.detect_scale = 1. / total_DC_gain;
flatmax@612: 
flatmax@612: }
flatmax@612: 
flatmax@612: int OK AGC::Design_FIR_coeffs(int n_taps, FP_TYPE var, FP_TYPE mn, int n_iter, Array<FP_TYPE,Dynamic,1> &FIR) {
flatmax@612: // reduce mean and variance of smoothing distribution by n_iterations:
flatmax@612:     mn = mn / (FP_TYPE)n_iter;
flatmax@612:     var = var / (FP_TYPE)n_iter;
flatmax@612:     switch (n_taps) {
flatmax@612:     case 3:
flatmax@612:         // based on solving to match mean and variance of [a, 1-a-b, b]:
flatmax@612:         a = (var + mn*mn - mn) / 2.;
flatmax@612:         b = (var + mn*mn + mn) / 2.;
flatmax@612:         FIR.resize(3,1);
flatmax@612:         FIR<<a, 1. - a - b, b;
flatmax@612:         OK = FIR(2) >= 0.2;
flatmax@612:     case 5
flatmax@612:             // based on solving to match [a/2, a/2, 1-a-b, b/2, b/2]:
flatmax@612:             a = ((var + mn*mn)*2./5. - mn*2./3.) / 2.;
flatmax@612:         b = ((var + mn*mn)*2./5. + mn*2./3.) / 2.;
flatmax@612:         // first and last coeffs are implicitly duplicated to make 5-point FIR:
flatmax@612:         FIR.resize(5,1);
flatmax@612:         FIR<<a/2., 1. - a - b, b/2.;
flatmax@612:         OK = FIR(2) >= 0.1;
flatmax@612:     default:
flatmax@612:         cerr<<"Bad n_taps in AGC_spatial_FIR"<<endl;
flatmax@612:         exit(AGC_FIR_TAP_COUNT_ERROR);
flatmax@612:     }
flatmax@612: }