Mercurial > hg > aimc

changeset 601:d838de2ce1b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added AGC::designAGC() This new method is not debugged ... that is the next step.
author flatmax
date Tue, 02 Apr 2013 08:38:23 +0000
parents 2d18671876c8
children 256caa099e32
files C++/AGC.C C++/AGC.H C++/AGCCoeff.H C++/CAR.H C++/CARFACCommon.H C++/Ear.C C++/codeblocks/CARFAC.layout
diffstat 7 files changed, 164 insertions(+), 34 deletions(-) [+]
line wrap: on
line diff
--- a/C++/AGC.C	Tue Feb 26 10:43:26 2013 +0000
+++ b/C++/AGC.C	Tue Apr 02 08:38:23 2013 +0000
@@ -23,12 +23,129 @@
 
 #include "AGC.H"
 
-AGC::AGC()
-{
-    //ctor
+AGC::AGC() {
 }
 
-AGC::~AGC()
-{
-    //dtor
+AGC::~AGC() {
 }
+
+void AGC::designAGC(FP_TYPE fs, int n_ch) {
+    int n_AGC_stages = 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 (in double exp. version)
+//AGC1_scales = AGC_params.AGC1_scales;
+//AGC2_scales = AGC_params.AGC2_scales;
+
+    coeffs.AGC_epsilon = Array<FP_TYPE, 1, Dynamic>::Zero(1, n_AGC_stages);  // the 1/(tau*fs) roughly
+    FP_TYPE decim = 1.;
+//AGC_coeffs.decimation = AGC_params.decimation;
+
+    FP_TYPE total_DC_gain = 0.;
+    for (int stage = 1; stage<=n_AGC_stages; stage++) {
+        FP_TYPE tau = params.time_constants(stage-1); // time constant in seconds
+        decim = decim * params.decimation(stage-1); // net decim to this stage
+        // epsilon is how much new input to take at each update step:
+        coeffs.AGC_epsilon(stage-1) = 1. - exp(-decim / (tau * fs));
+        // effective number of smoothings in a time constant:
+        FP_TYPE ntimes = tau * (fs / decim);  // typically 5 to 50
+
+        // decide on target spread (variance) and delay (mean) of impulse
+        // response as a distribution to be convolved ntimes:
+        // TODO (dicklyon): specify spread and delay instead of scales???
+        FP_TYPE delay = (param.AGC2_scales(stage-1) - param.AGC1_scales(stage-1)) / ntimes;
+        FP_TYPE spread_sq = (param.AGC1_scales(stage-1).pow(2.) + param.AGC2_scales(stage-1).pow(2)) / ntimes;
+
+        // get pole positions to better match intended spread and delay of
+        // [[geometric distribution]] in each direction (see wikipedia)
+        FP_TYPE u = 1. + 1. / spread_sq; // these are based on off-line algebra hacking.
+        FP_TYPE p = u - sqrt(pow(u,2.) - 1.); // pole that would give spread if used twice.
+        FP_TYPE dp = delay * (1. - 2.*p +pow(p,2.))/2.;
+                              FP_TYPE polez1 = p - dp;
+                              FP_TYPE polez2 = p + dp;
+                              coeffs.AGC_polez1(stage) = polez1;
+                              coeffs.AGC_polez2(stage) = polez2;
+
+                              // try a 3- or 5-tap FIR as an alternative to the double exponential:
+                              Array<FP_TYPE,1, Dynamic> AGC_spatial_FIR;
+                              int n_taps = 0;
+                              int FIR_OK = 0;
+                              int n_iterations = 1;
+        while (~FIR_OK) {
+        switch (n_taps) {
+            case 0:
+                // first attempt a 3-point FIR to apply once:
+                n_taps = 3;
+                break;
+            case 3:
+                // second time through, go wider but stick to 1 iteration
+                n_taps = 5;
+                break;
+            case 5:
+                // apply FIR multiple times instead of going wider:
+                n_iterations = n_iterations + 1;
+                if (n_iterations > 16) {
+                    cerr<<"Too many n_iterations in CARFAC_DesignAGC"<<endl;
+                    exit(AGC_DESIGN_ITERATION_ERROR);
+                }
+                break;
+            default:
+                // to do other n_taps would need changes in CARFAC_Spatial_Smooth
+                // and in Design_FIR_coeffs
+                cerr<<"Bad n_taps in CARFAC_DesignAGC"<<endl;
+                exit(AGC_DESIGN_TAPS_OOB_ERROR);
+                break;
+            }
+            FIR_OK = Design_FIR_coeffs(n_taps, spread_sq, delay, n_iterations,AGC_spatial_FIR);
+        }
+        // when FIR_OK, store the resulting FIR design in coeffs:
+        coeff.AGC_spatial_iterations(stage-1) = n_iterations;
+        coeff.AGC_spatial_FIR.col(stage-1).block(0,AGC_spatial_FIR.size()) = AGC_spatial_FIR;
+        coeff.AGC_spatial_n_taps(stage-1) = n_taps;
+
+        // accumulate DC gains from all the stages, accounting for stage_gain:
+        total_DC_gain = total_DC_gain + params.AGC_stage_gain.pow(stage-1);
+
+        // TODO (dicklyon) -- is this the best binaural mixing plan?
+        if (stage == 1)
+        coeff.AGC_mix_coeffs(stage-1) = 0.;
+        else
+            coeff.AGC_mix_coeffs(stage-1) = param.AGC_mix_coeff / (tau * (fs / decim));
+        }
+
+coeff.AGC_gain = total_DC_gain;
+
+// adjust the detect_scale to be the reciprocal DC gain of the AGC filters:
+AGC_coeffs.detect_scale = 1. / total_DC_gain;
+
+}
+
+int OK AGC::Design_FIR_coeffs(int n_taps, FP_TYPE var, FP_TYPE mn, int n_iter, Array<FP_TYPE,Dynamic,1> &FIR) {
+// reduce mean and variance of smoothing distribution by n_iterations:
+    mn = mn / (FP_TYPE)n_iter;
+    var = var / (FP_TYPE)n_iter;
+    switch (n_taps) {
+    case 3:
+        // based on solving to match mean and variance of [a, 1-a-b, b]:
+        a = (var + mn*mn - mn) / 2.;
+        b = (var + mn*mn + mn) / 2.;
+        FIR.resize(3,1);
+        FIR<<a, 1. - a - b, b;
+        OK = FIR(2) >= 0.2;
+    case 5
+            // based on solving to match [a/2, a/2, 1-a-b, b/2, b/2]:
+            a = ((var + mn*mn)*2./5. - mn*2./3.) / 2.;
+        b = ((var + mn*mn)*2./5. + mn*2./3.) / 2.;
+        // first and last coeffs are implicitly duplicated to make 5-point FIR:
+        FIR.resize(5,1);
+        FIR<<a/2., 1. - a - b, b/2.;
+        OK = FIR(2) >= 0.1;
+    default:
+        cerr<<"Bad n_taps in AGC_spatial_FIR"<<endl;
+        exit(AGC_FIR_TAP_COUNT_ERROR);
+    }
+}
--- a/C++/AGC.H	Tue Feb 26 10:43:26 2013 +0000
+++ b/C++/AGC.H	Tue Apr 02 08:38:23 2013 +0000
@@ -31,6 +31,12 @@
     public:
         AGC();
         virtual ~AGC();
+
+    /** Method to design the automatic gain control coefficients
+    \param fs The sample rate in Hz
+    \param n_ch The number of channels. The CAR.coeff are one row per channel
+    */
+    void designAGC(FP_TYPE fs, int n_ch);
 };
 
 #endif // AGC_H_
--- a/C++/AGCCoeff.H	Tue Feb 26 10:43:26 2013 +0000
+++ b/C++/AGCCoeff.H	Tue Apr 02 08:38:23 2013 +0000
@@ -1,5 +1,5 @@
 // Copyright 2013 Matt R. Flax <flatmax\@> All Rights Reserved.
-// Author Matt Flax <flatmax@>
+// Author Matt Flax <flatmax@>s
 //
 // This C++ file is part of an implementation of Lyon's cochlear model:
 // "Cascade of Asymmetric Resonators with Fast-Acting Compression"
@@ -29,6 +29,8 @@
 */
 class AGCCoeff : public Coefficients {
 public:
+    Array<FP_TYPE, 1, Dynamic> AGC_epsilon; ///< the 1/(tau*fs) roughly
+
     AGCCoeff();
     virtual ~AGCCoeff();
 };
--- a/C++/CAR.H	Tue Feb 26 10:43:26 2013 +0000
+++ b/C++/CAR.H	Tue Apr 02 08:38:23 2013 +0000
@@ -36,6 +36,7 @@
 
     /** Method to design the auditory filter coefficients
     \param fs The sample rate in Hz
+    \param n_ch The number of channels. The CAR.coeff are one row per channel
     */
     void designFilters(FP_TYPE fs, int n_ch);
 
--- a/C++/CARFACCommon.H	Tue Feb 26 10:43:26 2013 +0000
+++ b/C++/CARFACCommon.H	Tue Apr 02 08:38:23 2013 +0000
@@ -30,6 +30,10 @@
 
 #include "PsychoAcoustics.H"
 
+#define AGC_DESIGN_ITERATION_ERROR -100 ///< Error when designing the AGC filter taps
+#define AGC_DESIGN_TAPS_OOB_ERROR -101 ///< The number of taps requested have not been accounted for in the code yet.
+#define AGC_FIR_TAP_COUNT_ERROR -102 ///< The number of taps passed to the AGC::Design_FIR_coeffs method are not handled
+
 /**
     \mainpage CARFAC C++
 
--- a/C++/Ear.C	Tue Feb 26 10:43:26 2013 +0000
+++ b/C++/Ear.C	Tue Apr 02 08:38:23 2013 +0000
@@ -60,7 +60,7 @@
     max_channels_per_octave = (FP_TYPE)(log(2.) / log(car.pole_freqs[0]/car.pole_freqs[1]));
 
     // convert to include an ear_array, each w coeffs and state...
-    car.designFilters(fs, n_ch); //(car.param. fs, car.pole_freqs);
-    //AGC.designAGC(CF_AGC_params, fs, n_ch);
+    car.designFilters(fs, n_ch);
+    AGC.designAGC(fs, n_ch);
     //IHC.designIHC(CF_IHC_params, fs, n_ch);
 }
--- a/C++/codeblocks/CARFAC.layout	Tue Feb 26 10:43:26 2013 +0000
+++ b/C++/codeblocks/CARFAC.layout	Tue Apr 02 08:38:23 2013 +0000
@@ -1,23 +1,23 @@
 <?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
 <CodeBlocks_layout_file>
 	<ActiveTarget name="Debug" />
-	<File name="../AGC.C" open="0" top="0" tabpos="2">
-		<Cursor position="94" topLine="0" />
+	<File name="../AGC.C" open="1" top="1" tabpos="7">
+		<Cursor position="6353" topLine="125" />
 	</File>
-	<File name="../AGC.H" open="0" top="0" tabpos="3">
-		<Cursor position="94" topLine="0" />
+	<File name="../AGC.H" open="1" top="0" tabpos="6">
+		<Cursor position="1419" topLine="14" />
 	</File>
 	<File name="../AGCCoeff.C" open="0" top="0" tabpos="4">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../AGCCoeff.H" open="0" top="0" tabpos="5">
-		<Cursor position="94" topLine="0" />
+	<File name="../AGCCoeff.H" open="1" top="0" tabpos="9">
+		<Cursor position="1122" topLine="10" />
 	</File>
 	<File name="../AGCParam.C" open="0" top="0" tabpos="5">
 		<Cursor position="94" topLine="1" />
 	</File>
-	<File name="../AGCParam.H" open="0" top="0" tabpos="5">
-		<Cursor position="1149" topLine="13" />
+	<File name="../AGCParam.H" open="1" top="0" tabpos="8">
+		<Cursor position="2004" topLine="25" />
 	</File>
 	<File name="../AGCState.C" open="0" top="0" tabpos="8">
 		<Cursor position="94" topLine="0" />
@@ -25,23 +25,23 @@
 	<File name="../AGCState.H" open="0" top="0" tabpos="9">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../CAR.C" open="1" top="1" tabpos="1">
-		<Cursor position="3388" topLine="65" />
+	<File name="../CAR.C" open="0" top="0" tabpos="1">
+		<Cursor position="2218" topLine="46" />
 	</File>
-	<File name="../CAR.H" open="1" top="0" tabpos="2">
-		<Cursor position="1679" topLine="26" />
+	<File name="../CAR.H" open="0" top="0" tabpos="2">
+		<Cursor position="1537" topLine="21" />
 	</File>
 	<File name="../CARCoeff.C" open="0" top="0" tabpos="12">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../CARCoeff.H" open="1" top="0" tabpos="8">
+	<File name="../CARCoeff.H" open="0" top="0" tabpos="8">
 		<Cursor position="1171" topLine="20" />
 	</File>
-	<File name="../CARFACCommon.H" open="1" top="0" tabpos="7">
-		<Cursor position="2156" topLine="46" />
+	<File name="../CARFACCommon.H" open="1" top="0" tabpos="5">
+		<Cursor position="1559" topLine="15" />
 	</File>
 	<File name="../CARParam.C" open="0" top="0" tabpos="7">
-		<Cursor position="94" topLine="16" />
+		<Cursor position="1321" topLine="16" />
 	</File>
 	<File name="../CARParam.H" open="0" top="0" tabpos="3">
 		<Cursor position="2340" topLine="31" />
@@ -55,19 +55,19 @@
 	<File name="../Coefficients.C" open="0" top="0" tabpos="18">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../Coefficients.H" open="1" top="0" tabpos="10">
+	<File name="../Coefficients.H" open="0" top="0" tabpos="10">
 		<Cursor position="1139" topLine="12" />
 	</File>
-	<File name="../Ear.C" open="1" top="0" tabpos="3">
-		<Cursor position="1709" topLine="40" />
+	<File name="../Ear.C" open="1" top="0" tabpos="1">
+		<Cursor position="2303" topLine="40" />
 	</File>
-	<File name="../Ear.H" open="1" top="0" tabpos="4">
+	<File name="../Ear.H" open="1" top="0" tabpos="2">
 		<Cursor position="1086" topLine="18" />
 	</File>
 	<File name="../EarComponent.C" open="0" top="0" tabpos="6">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../EarComponent.H" open="1" top="0" tabpos="9">
+	<File name="../EarComponent.H" open="0" top="0" tabpos="9">
 		<Cursor position="1541" topLine="26" />
 	</File>
 	<File name="../Ears.C" open="0" top="0" tabpos="24">
@@ -100,13 +100,13 @@
 	<File name="../IHCState.H" open="0" top="0" tabpos="33">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../Makefile" open="1" top="0" tabpos="6">
+	<File name="../Makefile" open="1" top="0" tabpos="4">
 		<Cursor position="1159" topLine="9" />
 	</File>
 	<File name="../Parameters.C" open="0" top="0" tabpos="35">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../Parameters.H" open="1" top="0" tabpos="11">
+	<File name="../Parameters.H" open="0" top="0" tabpos="11">
 		<Cursor position="1125" topLine="12" />
 	</File>
 	<File name="../PsychoAcoustics.C" open="0" top="0" tabpos="8">
@@ -118,10 +118,10 @@
 	<File name="../State.C" open="0" top="0" tabpos="37">
 		<Cursor position="94" topLine="0" />
 	</File>
-	<File name="../State.H" open="1" top="0" tabpos="12">
+	<File name="../State.H" open="0" top="0" tabpos="12">
 		<Cursor position="1095" topLine="12" />
 	</File>
-	<File name="../test/EarTest.C" open="1" top="0" tabpos="5">
+	<File name="../test/EarTest.C" open="1" top="0" tabpos="3">
 		<Cursor position="989" topLine="4" />
 	</File>
 </CodeBlocks_layout_file>