Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/SeneffEar.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2006 George Tzanetakis <gtzan@cs.uvic.ca>
00003 **
00004 ** This program is free software; you can redistribute it and/or modify
00005 ** it under the terms of the GNU General Public License as published by
00006 ** the Free Software Foundation; either version 2 of the License, or
00007 ** (at your option) any later version.
00008 **
00009 ** This program is distributed in the hope that it will be useful,
00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 ** GNU General Public License for more details.
00013 **
00014 ** You should have received a copy of the GNU General Public License
00015 ** along with this program; if not, write to the Free Software
00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 */
00018 
00019 #include "SeneffEar.h"
00020 #include "../common_source.h"
00021 
00022 #include <marsyas/fft.h>
00023 
00024 #include <sstream>
00025 #include <algorithm>
00026 
00027 using namespace std;
00028 using namespace Marsyas;
00029 
00030 SeneffEar::SeneffEar(mrs_string name):MarSystem("SeneffEar",name)
00031 {
00032   //type_ = "SeneffEar";
00033   //name_ = name;
00034 
00035   firstUpdate = true;
00036   fs = 0.0f;
00037   stage = 4;
00038 
00039   addControls();
00040 }
00041 
00042 SeneffEar::~SeneffEar()
00043 {
00044 }
00045 
00046 MarSystem*
00047 SeneffEar::clone() const
00048 {
00049   return new SeneffEar(*this);
00050 }
00051 
00052 void
00053 SeneffEar::addControls()
00054 {
00055   addctrl("mrs_natural/stage", stage);
00056 }
00057 
00058 void
00059 SeneffEar::myUpdate(MarControlPtr sender)
00060 {
00061   (void) sender;  //suppress warning of unused parameter(s)
00062   MRSDIAG("SeneffEar.cpp - SeneffEar:myUpdate");
00063 
00064   ostringstream name;
00065 
00066   //first time update
00067   if (firstUpdate) {
00068     //PreemphasisRTheta coefficients initialization
00069     stringstream matrix;
00070     matrix << "# MARSYAS realvec" << endl;
00071     matrix << "# Size = 8" << endl;
00072     matrix << endl;
00073     matrix << endl;
00074     matrix << "# type: matrix" << endl;
00075     matrix << "# rows: 4" << endl;
00076     matrix << "# columns: 2" << endl;
00077     matrix << "0.86 3.1148863";
00078     matrix << "0.99 0.0";
00079     matrix << "0.5  0.0";
00080     matrix << "0.95 3.14159";
00081     matrix << endl;
00082     matrix << "# Size = 8" << endl;
00083     matrix << "# MARSYAS realvec" << endl;
00084     matrix >> PreemphasisRThetaCoeffs;
00085 
00086     //FilterBankRTheta coefficients initialization
00087     matrix.clear();
00088     matrix << "# MARSYAS realvec" << endl;
00089     matrix << "# Size = 200" << endl;
00090     matrix << endl;
00091     matrix << endl;
00092     matrix << "# type: matrix" << endl;
00093     matrix << "# rows: 40" << endl;
00094     matrix << "# columns: 5" << endl;
00095     //         R-z       Theta-z       Radius       Theta       R-z2
00096     matrix << "0.0       3.14159       0.740055     2.633909    0.8";
00097     matrix << "0.86      2.997077      0.753637     2.178169    0.8";
00098     matrix << "0.86      2.879267      0.775569     1.856744    0.8";
00099     matrix << "0.86      2.761458      0.798336     1.617919    0.8";
00100     matrix << "0.86      2.643648      0.819169     1.433496    0.8";
00101     matrix << "0.86      2.525839      0.837158     1.286795    0.8";
00102     matrix << "0.8       2.964876      0.852598     1.167321    0.8";
00103     matrix << "0.86      2.408029      0.865429     1.068141    0.8";
00104     matrix << "0.86      2.29022       0.876208     0.984489    0.8";
00105     matrix << "0.86      2.17241       0.885329     0.912985    0.8";
00106     matrix << "0.86      2.054601      0.893116     0.851162    0.8";
00107     matrix << "0.86      1.936791      0.899823     0.797179    0.8";
00108     matrix << "0.8       2.788161      0.906118     0.749633    0.8";
00109     matrix << "0.86      1.818981      0.911236     0.70744     0.8";
00110     matrix << "0.86      1.701172      0.915747     0.669742    0.8";
00111     matrix << "0.86      1.583362      0.919753     0.635858    0.8";
00112     matrix << "0.86      1.465552      0.923335     0.605237    0.8";
00113     matrix << "0.86      1.347743      0.926565     0.57743     0.8";
00114     matrix << "0.8       2.611447      0.929914     0.552065    0.8";
00115     matrix << "0.86      1.229933      0.932576     0.528834    0.8";
00116     matrix << "0.86      1.112123      0.944589     0.487783    0.75";
00117     matrix << "0.86      0.994314      0.957206     0.452645    0.660714";
00118     matrix << "0.86      0.876504      0.956548     0.42223     0.672143";
00119     matrix << "0.86      0.758694      0.956653     0.395644    0.682143";
00120     matrix << "0.8       2.434732      0.956518     0.372208    0.690966";
00121     matrix << "0.86      0.640885      0.956676     0.351393    0.69881";
00122     matrix << "0.86      0.523075      0.956741     0.316044    0.712143";
00123     matrix << "0.8       2.258018      0.956481     0.287157    0.723052";
00124     matrix << "0.8       2.081304      0.956445     0.263108    0.732143";
00125     matrix << "0.8       1.904589      0.956481     0.242776    0.739835";
00126     matrix << "0.86      0.405265      0.958259     0.217558    0.749384";
00127     matrix << "0.8       1.727875      0.963083     0.197086    0.757143";
00128     matrix << "0.8       1.55116       0.969757     0.175115    0.769048";
00129     matrix << "0.8       1.374446      0.97003      0.153697    0.780662";
00130     matrix << "0.8       1.197732      0.970382     0.134026    0.791337";
00131     matrix << "0.8       1.021017      0.970721     0.118819    0.799596";
00132     matrix << "0.8       1.5           0.970985     0.106711    0.8";
00133     matrix << "0.8       1.2           0.971222     0.096843    0.8";
00134     matrix << "0.8       1.0           0.97144      0.088645    0.8";
00135     matrix << "0.8       0.9           0.971645     0.081727    0.8";
00136     matrix << endl;
00137     matrix << "# Size = 200" << endl;
00138     matrix << "# MARSYAS realvec" << endl;
00139     matrix >> FilterBankRThetaCoeffs;
00140 
00141     realvec a,b,c;
00142 
00143     //SeneffPreemphasis coefficients computation
00144     a.create(3);
00145     b.create(3);
00146     a(0) = PreemphasisRThetaCoeffs(0,0)*PreemphasisRThetaCoeffs(0,0);
00147     a(1) = -2.0f*PreemphasisRThetaCoeffs(0,0)*cos(PreemphasisRThetaCoeffs(0,1));
00148     a(2) = 1.0f;
00149     b(0) = PreemphasisRThetaCoeffs(1,0)*PreemphasisRThetaCoeffs(1,0);
00150     b(1) = -2.0f*PreemphasisRThetaCoeffs(1,0)*cos(PreemphasisRThetaCoeffs(1,1));
00151     b(2) = 1.0f;
00152     polyConv(a, b, a);
00153     b(0) = PreemphasisRThetaCoeffs(2,0)*PreemphasisRThetaCoeffs(2,0);
00154     b(1) = -2.0f*PreemphasisRThetaCoeffs(2,0)*cos(PreemphasisRThetaCoeffs(2,1));
00155     b(2) = 1.0f;
00156     polyConv(a, b, a);
00157     b(0) = PreemphasisRThetaCoeffs(3,0)*PreemphasisRThetaCoeffs(3,0);
00158     b(1) = -2.0f*PreemphasisRThetaCoeffs(3,0)*cos(PreemphasisRThetaCoeffs(3,1));
00159     b(2) = 1.0f;
00160     polyConv(a, b, c);
00161     SeneffPreemphasisCoeffs.create(c.getSize());
00162 
00163     for (mrs_natural i = 0; i < c.getSize(); ++i)
00164       SeneffPreemphasisCoeffs(i) = c(c.getSize() - i - 1);
00165 
00166     //SeneffFilterBank coefficients computation
00167     channels = FilterBankRThetaCoeffs.getRows();
00168     SeneffFilterBankCoeffs.create(channels, 5);
00169     mrs_real sum;
00170     for (mrs_natural i = 0; i < channels; ++i) {
00171       SeneffFilterBankCoeffs(i,0) = 1.0f;
00172       SeneffFilterBankCoeffs(i,1) = -2.0f*FilterBankRThetaCoeffs(i,0)*cos(FilterBankRThetaCoeffs(i,1));
00173       SeneffFilterBankCoeffs(i,2) = FilterBankRThetaCoeffs(i,0)*FilterBankRThetaCoeffs(i,0);
00174       sum = SeneffFilterBankCoeffs(i,0) + SeneffFilterBankCoeffs(i,1) + SeneffFilterBankCoeffs(i,2);
00175       SeneffFilterBankCoeffs(i,0) /= sum;
00176       SeneffFilterBankCoeffs(i,1) /= sum;
00177       SeneffFilterBankCoeffs(i,2) /= sum;
00178     }
00179 
00180     //SeneffForward & SeneffBackward coefficients computation
00181     SeneffForwardCoeffs.create(mrs_natural(5), static_cast<mrs_natural>(channels));
00182     SeneffBackwardCoeffs.create(mrs_natural(5), static_cast<mrs_natural>(channels));
00183     a.create(3),b.create(3);
00184     for (mrs_natural j = 0; j < channels; j++) {
00185       a(0) = FilterBankRThetaCoeffs(j,4)*FilterBankRThetaCoeffs(j,4);
00186       a(1) = -2.0f*FilterBankRThetaCoeffs(j,4)*cos(FilterBankRThetaCoeffs(j,3)/2);
00187       a(2) = 1.0f;
00188       polyConv(a, a, c);
00189       for (mrs_natural i = 0; i < 5; ++i) SeneffForwardCoeffs(i,j) = c(5 - i - 1);
00190       b(0) = FilterBankRThetaCoeffs(j,2)*FilterBankRThetaCoeffs(j,2);
00191       b(1) = -2.0f*FilterBankRThetaCoeffs(j,2)*cos(FilterBankRThetaCoeffs(j,3));
00192       b(2) = 1.0f;
00193       polyConv(b, b, c);
00194       for (mrs_natural i = 0; i < 5; ++i) SeneffBackwardCoeffs(i,j) = c(5 - i - 1);
00195     }
00196 
00197     //SeneffPreemphasis filter creation and coefficients initialization
00198     SeneffPreemphasisFilter = new Filter("SeneffPreemphasisFilter");
00199     a.create(1);
00200     a(0) = 1.0f;
00201     SeneffPreemphasisFilter->setctrl("mrs_realvec/ncoeffs", SeneffPreemphasisCoeffs);
00202     SeneffPreemphasisFilter->setctrl("mrs_realvec/dcoeffs", a);
00203 
00204     //SeneffFilterBank filters cascade creation and coefficients initialization
00205     SeneffFilterBank = new Cascade("SeneffFilterBank");
00206     Filter* filter;
00207     a.create(3);
00208     b.create(3);
00209     for (mrs_natural i = 0; i < channels; ++i) {
00210       name.clear();
00211       name.str("");
00212       name << "filter_" << i;
00213       filter = new Filter(name.str());
00214       b(0) = SeneffFilterBankCoeffs(i,0);
00215       b(1) = SeneffFilterBankCoeffs(i,1);
00216       b(2) = SeneffFilterBankCoeffs(i,2);
00217       a(0) = 1.0f;
00218       a(1) = SeneffFilterBankCoeffs(i,3);
00219       a(2) = SeneffFilterBankCoeffs(i,4);
00220       filter->setctrl("mrs_realvec/ncoeffs", b);
00221       filter->setctrl("mrs_realvec/dcoeffs", a);
00222       SeneffFilterBank->addMarSystem(filter);
00223     }
00224 
00225     //resonatorFilter creation
00226     resonatorFilter = new Parallel("resonatorFilter");
00227     a.create(5);
00228     b.create(5);
00229     for (mrs_natural i = 0; i < channels; ++i) {
00230       name.clear();
00231       name.str("");
00232       name << "filter_" << i;
00233       filter = new Filter(name.str());
00234       for (mrs_natural j = 0; j < 5; j++) b(j) = SeneffForwardCoeffs(j,i);
00235       for (mrs_natural j = 0; j < 5; j++) a(j) = SeneffBackwardCoeffs(j,i);
00236       filter->setctrl("mrs_realvec/ncoeffs", b);
00237       filter->setctrl("mrs_realvec/dcoeffs", a);
00238       resonatorFilter->addMarSystem(filter);
00239     }
00240 
00241     //run an impulse through the preemphasis filter
00242     realvec impulse;
00243     impulse.create(mrs_natural(1), mrs_natural(256));
00244     impulse(0) = 1.0f;
00245     realvec y0;
00246     y0.create(impulse.getRows(), impulse.getCols());
00247     SeneffPreemphasisFilter->setctrl("mrs_natural/inSamples", impulse.getCols());
00248     SeneffPreemphasisFilter->setctrl("mrs_natural/inObservations", impulse.getRows());
00249     SeneffPreemphasisFilter->update();
00250     SeneffPreemphasisFilter->process(impulse, y0);
00251 
00252     //then through the cascade of zeros
00253     realvec y1;
00254     y1.create(y0.getRows()*channels, y0.getCols());
00255     SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inSamples", y0.getCols());
00256     SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inObservations", y0.getRows());
00257     SeneffFilterBank->update();
00258     SeneffFilterBank->process(y0, y1);
00259 
00260     //Finally through each of the resonator filters
00261     y.create(y1.getRows(), y1.getCols());
00262     resonatorFilter->setctrl("Filter/filter_0/mrs_natural/inSamples", y1.getCols());
00263     for (mrs_natural i = 0; i < channels; ++i) {
00264       name.clear();
00265       name.str("");
00266       name << "Filter/" << "filter_" << i << "/mrs_natural/inObservations";
00267       resonatorFilter->setctrl(name.str(), y0.getRows());
00268     }
00269     resonatorFilter->update();
00270     resonatorFilter->process(y1, y);
00271 
00272     //reset the state vector of each filter
00273     realvec state;
00274     state.create(SeneffPreemphasisFilter->getctrl("mrs_realvec/state")->to<mrs_realvec>().getRows(),SeneffPreemphasisFilter->getctrl("mrs_realvec/state")->to<mrs_realvec>().getCols());
00275     SeneffPreemphasisFilter->setctrl("mrs_natural/stateUpdate", mrs_natural(1));
00276     SeneffPreemphasisFilter->updControl("mrs_realvec/state", state);
00277     SeneffPreemphasisFilter->setctrl("mrs_natural/stateUpdate", mrs_natural(0));
00278 
00279     state.create(SeneffFilterBank->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getRows(),SeneffFilterBank->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getCols());
00280     for (mrs_natural i = 0; i < channels; ++i) {
00281       name.clear();
00282       name.str("");
00283       name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
00284       SeneffFilterBank->setctrl(name.str(), mrs_natural(1));
00285       name.clear();
00286       name.str("");
00287       name << "Filter/" << "filter_" << i << "/mrs_realvec/state";
00288       SeneffFilterBank->updControl(name.str(), state);
00289       name.clear();
00290       name.str("");
00291       name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
00292       SeneffFilterBank->setctrl(name.str(), mrs_natural(0));
00293     }
00294 
00295     state.create(resonatorFilter->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getRows(),resonatorFilter->getctrl("Filter/filter_0/mrs_realvec/state")->to<mrs_realvec>().getCols());
00296     for (mrs_natural i = 0; i < channels; ++i) {
00297       name.clear();
00298       name.str("");
00299       name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
00300       resonatorFilter->setctrl(name.str(), mrs_natural(1));
00301       name.clear();
00302       name.str("");
00303       name << "Filter/" << "filter_" << i << "/mrs_realvec/state";
00304       resonatorFilter->updControl(name.str(), state);
00305       name.clear();
00306       name.str("");
00307       name << "Filter/" << "filter_" << i << "/mrs_natural/stateUpdate";
00308       resonatorFilter->setctrl(name.str(), mrs_natural(0));
00309     }
00310 
00311     //constants initialization
00312     hwrA = 10.0f;
00313     hwrB = 65.0f;
00314     hwrG = 2.35f;
00315 
00316     Cn.create(channels);
00317 
00318     //lowPassFilter creation and coefficients initialization
00319     lpAlpha = 0.209611f;
00320     lowPassFilter = new Filter("lowPassFilter");
00321     a.create(2);
00322     a(0) = -lpAlpha;
00323     a(1) = 1.0f;
00324     polyConv(a, a, a);
00325     polyConv(a, a, a);
00326     polyFlip(a);
00327     b.create(1);
00328     b(0) = a.sum();
00329     lowPassFilter->setctrl("mrs_realvec/ncoeffs", b);
00330     lowPassFilter->setctrl("mrs_realvec/dcoeffs", a);
00331 
00332     //Seneff's Adaptive Gain Control filter creation and coefficients initialization
00333     initial_yn = (mrs_real)0.23071276;
00334     alpha_agc = (mrs_real)0.979382181;
00335     kagc = (mrs_real)0.002;
00336     AGCfilter = new Filter("AGCfilter");
00337     a.create(1);
00338     a(0) = alpha_agc;
00339     b.create(2);
00340     b(1) = 1.0f - alpha_agc;
00341     AGCfilter->setctrl("mrs_realvec/ncoeffs", b);
00342     AGCfilter->setctrl("mrs_realvec/dcoeffs", a);
00343     state.create(channels, mrs_natural(1));
00344     state.setval(initial_yn);
00345     AGCfilter->setctrl("mrs_natural/inObservations", channels);
00346     AGCfilter->setctrl("mrs_natural/stateUpdate", mrs_natural(1));
00347     AGCfilter->updControl("mrs_realvec/state", state);
00348     AGCfilter->setctrl("mrs_natural/stateUpdate", mrs_natural(0));
00349 
00350     firstUpdate = false;
00351   }//end first update
00352 
00353   //--------------------------------------------------------------------------------
00354   //compute the FFT then find the gain peak then divide each forward polynomial by the maximum gain (to normalize) and then multiply by the desired low frequency roll-off.
00355   if (fs != getctrl("mrs_real/israte")->to<mrs_real>()) {
00356     fs = getctrl("mrs_real/israte")->to<mrs_real>();
00357     fft fft;
00358     realvec Y(y.getCols());
00359     mrs_real *Yd = Y.getData();
00360     mrs_real abs2(1.0f),maxAbs2(0.0f);
00361     mrs_real gain;
00362     mrs_real rolloff;
00363     SeneffForwardCoeffsNormalized.create(SeneffForwardCoeffs.getRows(), SeneffForwardCoeffs.getCols());
00364     for (mrs_natural i = 0; i < channels; ++i) {
00365       for (mrs_natural j = 0; j < Y.getCols(); j++) Y(j) = y(i,j);
00366       fft.rfft(Yd, Y.getCols()/2, FFT_FORWARD);
00367       maxAbs2 = Y(0)*Y(0);
00368       for (mrs_natural j = 1; j < y.getCols()/2; j++) if ((abs2 = Y(2*j)*Y(2*j) + Y(2*j+1)*Y(2*j+1)) > maxAbs2) maxAbs2 = abs2;
00369       gain = 1/(y.getCols()*sqrt(maxAbs2));
00370       rolloff = min((mrs_real)((FilterBankRThetaCoeffs(i,3)/PI*fs/2)/1600),(mrs_real)1.0);
00371       for (mrs_natural j = 0; j < SeneffForwardCoeffsNormalized.getRows(); j++) SeneffForwardCoeffsNormalized(j,i) = SeneffForwardCoeffs(j,i)*gain*rolloff;
00372     }
00373     //then update the resonatorFilter
00374     realvec b(5);
00375     for (mrs_natural i = 0; i < channels; ++i) {
00376       name.clear();
00377       name.str("");
00378       name << "Filter/" << "filter_" << i << "/mrs_realvec/ncoeffs";
00379       for (mrs_natural j = 0; j < 5; j++) b(j) = SeneffForwardCoeffsNormalized(j,i);
00380       resonatorFilter->setctrl(name.str(), b);
00381     }
00382     resonatorFilter->update();
00383 
00384     Tua = (mrs_real)(58.3333/fs);
00385     Tub = (mrs_real)(8.3333/fs);
00386   }
00387 
00388   //--------------------------------------------------------------------------------
00389   //SeneffPreemphasis update
00390   SeneffPreemphasisFilter->setctrl("mrs_natural/inSamples", getctrl("mrs_natural/inSamples"));
00391   SeneffPreemphasisFilter->setctrl("mrs_natural/inObservations", getctrl("mrs_natural/inObservations"));
00392   SeneffPreemphasisFilter->setctrl("mrs_real/israte", getctrl("mrs_real/israte"));
00393   SeneffPreemphasisFilter->update();
00394   if ((int)slice_0.getSize() != SeneffPreemphasisFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * SeneffPreemphasisFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
00395     slice_0.create(SeneffPreemphasisFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), SeneffPreemphasisFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
00396   }
00397 
00398   //SeneffFilterBank filters cascade update
00399   SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inSamples", SeneffPreemphasisFilter->getctrl("mrs_natural/onSamples"));
00400   SeneffFilterBank->setctrl("Filter/filter_0/mrs_natural/inObservations", SeneffPreemphasisFilter->getctrl("mrs_natural/onObservations"));
00401   SeneffFilterBank->setctrl("Filter/filter_0/mrs_real/israte", SeneffPreemphasisFilter->getctrl("mrs_real/osrate"));
00402   SeneffFilterBank->update();
00403   if ((int)slice_1.getSize() != SeneffFilterBank->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * SeneffFilterBank->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
00404     slice_1.create(SeneffFilterBank->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), SeneffFilterBank->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
00405   }
00406 
00407   //resonatorFilter update
00408   resonatorFilter->setctrl("Filter/filter_0/mrs_natural/inSamples", SeneffFilterBank->getctrl("mrs_natural/onSamples"));
00409   for (mrs_natural i = 0; i < channels; ++i) {
00410     name.clear();
00411     name.str("");
00412     name << "Filter/" << "filter_" << i << "/mrs_natural/inObservations";
00413     resonatorFilter->setctrl(name.str(), SeneffFilterBank->getctrl("Filter/filter_0/mrs_natural/onObservations"));
00414   }
00415   resonatorFilter->setctrl("Filter/filter_0/mrs_real/israte", SeneffFilterBank->getctrl("mrs_real/osrate"));
00416   resonatorFilter->update();
00417   if ((int)slice_2.getSize() != resonatorFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * resonatorFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
00418     slice_2.create(resonatorFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), resonatorFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
00419   }
00420 
00421   //lowPassFilter update
00422   lowPassFilter->setctrl("mrs_natural/inSamples", resonatorFilter->getctrl("mrs_natural/onSamples"));
00423   lowPassFilter->setctrl("mrs_natural/inObservations", resonatorFilter->getctrl("mrs_natural/onObservations"));
00424   lowPassFilter->setctrl("mrs_real/israte", resonatorFilter->getctrl("mrs_real/osrate"));
00425   lowPassFilter->update();
00426   if ((int)slice_3.getSize() != lowPassFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>() * lowPassFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>()) {
00427     slice_3.create(lowPassFilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>(), lowPassFilter->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
00428   }
00429 
00430   //AGCfilter update
00431   AGCfilter->setctrl("mrs_natural/inSamples", lowPassFilter->getctrl("mrs_natural/onSamples"));
00432   AGCfilter->setctrl("mrs_natural/inObservations", lowPassFilter->getctrl("mrs_natural/onObservations"));
00433   AGCfilter->setctrl("mrs_real/israte", lowPassFilter->getctrl("mrs_real/osrate"));
00434   AGCfilter->update();
00435 
00436   //this update
00437   setctrl("mrs_natural/onSamples", AGCfilter->getctrl("mrs_natural/onSamples"));
00438   setctrl("mrs_natural/onObservations", AGCfilter->getctrl("mrs_natural/onObservations")->to<mrs_natural>());
00439   setctrl("mrs_real/osrate", AGCfilter->getctrl("mrs_real/osrate"));
00440 }
00441 
00442 void
00443 SeneffEar::polyConv(realvec& a, realvec& b, realvec& c)
00444 {
00445   mrs_natural la(a.getSize());
00446   mrs_natural lb(b.getSize());
00447   mrs_natural n = la + lb - 1;
00448 
00449   realvec ta(a); ta.stretch(n);
00450   realvec tb(b); tb.stretch(n);
00451   realvec tc;    tc.create(n);
00452 
00453   for (mrs_natural k = 0; k < n; k++) {
00454     for (mrs_natural i = 0; i <= k; ++i) {
00455       tc(k) += ta(i)*tb(k-i);
00456     }
00457   }
00458   if ((mrs_natural)c.getSize() != n) c.create(n);
00459   c = tc;
00460 }
00461 
00462 void
00463 SeneffEar::polyFlip(realvec& a)
00464 {
00465   mrs_natural la(a.getSize());
00466   realvec ta(a);
00467 
00468   for (mrs_natural i = 0; i < la; ++i) {
00469     a(i) = ta(la - i - 1);
00470   }
00471 }
00472 
00473 
00474 void
00475 SeneffEar::myProcess(realvec& in, realvec& out)
00476 {
00477   checkFlow(in, out);
00478   //lmartins: if (mute_) return;
00479   if(getctrl("mrs_bool/mute")->to<mrs_bool>()) return;
00480 
00481   mrs_natural s = 0;
00482   stage = getctrl("mrs_natural/stage")->to<mrs_natural>();
00483 
00484   SeneffPreemphasisFilter->process(in, slice_0);
00485   SeneffFilterBank->process(slice_0, slice_1);
00486   if (s++ == stage) {out = slice_1; return;}
00487   resonatorFilter->process(slice_1, slice_2);
00488   if (s++ == stage) {out = slice_2; return;}
00489   //Seneff's detector non-linearity
00490   for (mrs_natural i = 0; i < slice_2.getRows(); ++i) {
00491     for (mrs_natural j = 0; j < slice_2.getCols(); j++) {
00492       slice_2(i,j) = hwrA*atan(hwrB*max((mrs_real)0.0,slice_2(i,j))) + exp(hwrA*hwrB*min((mrs_real)0.0,slice_2(i,j)));
00493     }
00494   }
00495   //Seneff's short-term adaptation (a reservoir hair cell model)
00496   mrs_real flow;
00497   for (mrs_natural j = 0; j < slice_2.getCols(); j++) {
00498     for (mrs_natural i = 0; i < slice_2.getRows(); ++i) {
00499       flow = max((mrs_real)0.0,Tua*(slice_2(i,j)-Cn(i)));
00500       Cn(i) = Cn(i) + flow - Tub*Cn(i);
00501       slice_2(i,j) = flow;
00502     }
00503   }
00504   if (s++ == stage) {out = slice_2; return;}
00505   //Seneff's Low Pass filter (to model the loss of Synchrony)
00506   lowPassFilter->process(slice_2, slice_3);
00507   if (s++ == stage) {out = slice_3; return;}
00508   //Seneff's Adaptive Gain Control
00509   AGCfilter->process(slice_3, out);
00510   for (mrs_natural i = 0; i < out.getRows(); ++i) {
00511     for (mrs_natural j = 0; j < out.getCols(); j++) {
00512       out(i,j) = slice_3(i,j)/(1.0f + kagc*out(i,j));
00513     }
00514   }
00515 }