Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/LyonPassiveEar.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 
00020 #include "LyonPassiveEar.h"
00021 #include "../common_source.h"
00022 
00023 #include "Series.h"
00024 #include "Filter.h"
00025 #include "Cascade.h"
00026 #include "HalfWaveRectifier.h"
00027 #include <marsyas/basis.h>
00028 #include <sstream>
00029 
00030 using std::ostringstream;
00031 using std::abs;
00032 
00033 using namespace Marsyas;
00034 
00035 
00036 namespace Marsyas
00037 {
00049 class LyonAgc: public MarSystem
00050 {
00051 private:
00052   static mrs_real lyonEpsilonFromTauFS (mrs_real tau, mrs_real fs)
00053   {
00054     return 1.0 - exp (-1/tau/fs);
00055   }
00056   // copied and reformatted from auditory toolbox
00057   static void agc (const mrs_realvec input, mrs_realvec &output, mrs_realvec &state, mrs_real epsilon, mrs_real target, mrs_natural n)
00058   {
00059     mrs_natural i;
00060     mrs_real f,
00061              stateLimit             = 1. - 1e-1;
00062     mrs_real oneMinusEpsOverThree   = (1.0 - epsilon)/3.0;
00063     mrs_real epsOverTarget          = epsilon/target;
00064     mrs_real prevState              = state(0);
00065 
00066     for ( i = 0; i < n-1; ++i)
00067     {
00068       output(i) = fabs (input(i) * (1.0 - state(i)));
00069       f         = output(i) * epsOverTarget + oneMinusEpsOverThree * (prevState + state(i) + state(i+1));
00070 
00071       if (f > stateLimit)
00072         f = stateLimit;
00073 
00074       prevState = state(i);
00075       state(i)  = f;
00076     }
00077 
00078     // do the last iteration with a different state input
00079     output(i)   = fabs (input(i) * (1.0 - state(i)));
00080     f           = output(i) * epsOverTarget + oneMinusEpsOverThree * (prevState + state(i) + state(i));
00081 
00082     if (f > stateLimit)
00083       f = stateLimit;
00084 
00085     state(i)    = f;
00086   };
00087 
00088   void myUpdate(MarControlPtr sender)
00089   {
00091     MarSystem::myUpdate(sender);
00092 
00093     const mrs_natural   nStates = 4;
00094     mrs_real            fs      = getctrl ("mrs_real/israte")->to<mrs_real>();
00095     mrs_natural         nBands  = getctrl ("mrs_natural/numBands")->to<mrs_natural>();
00096 
00097     // allocate vectors
00098     agcParms_.create(2, nStates);
00099     state_.create(nBands, nStates);
00100     tmpOut_.create(nBands, 1);
00101 
00102     // set targets and epsilons
00103     agcParms_(0,0)  = .0032;
00104     agcParms_(0,1)  = .0016;
00105     agcParms_(0,2)  = .0008;
00106     agcParms_(0,3)  = .0004;
00107     agcParms_(1,0)  = lyonEpsilonFromTauFS (.64, fs);
00108     agcParms_(1,1)  = lyonEpsilonFromTauFS (.16, fs);
00109     agcParms_(1,2)  = lyonEpsilonFromTauFS (.04, fs);
00110     agcParms_(1,3)  = lyonEpsilonFromTauFS (.01, fs);
00111 
00112     // initialize
00113     state_.setval (0);
00114   };
00115 
00116   mrs_realvec   agcParms_,
00117               state_,
00118               tmpOut_;
00119 
00120   friend class LyonPassiveEar;
00121 
00122 public:
00123   LyonAgc(std::string name):MarSystem("LyonAgc", name)
00124   {
00125     addControls ();
00126   };
00127   ~LyonAgc() {};
00128 
00129   void addControls()
00130   {
00131     addctrl("mrs_natural/numBands", 1);
00132     setctrlState("mrs_natural/numBands", true);
00133   }
00134 
00135   MarSystem* clone() const
00136   {
00137     return new LyonAgc(*this);
00138   };
00139 
00140   void myProcess(realvec& in, realvec& out)
00141   {
00142     // outer agc loop (compare file agc.c from the auditory toolbox)
00143     for (mrs_natural t = 0; t < inSamples_; t++)
00144     {
00145       mrs_natural i,
00146                   nStages = agcParms_.getCols (),
00147                   nRows = in.getRows ();
00148       mrs_realvec state;
00149       in.getCol (t, tmpOut_);
00150 
00151       for (i = 0; i < nStages; ++i)
00152       {
00153         state_.getCol (i, state);
00154         agc (tmpOut_, tmpOut_, state, agcParms_(1,i), agcParms_(0,i), nRows);
00155         state_.setCol (i,state); // update agc state
00156       }
00157       out.setCol (t, tmpOut_);
00158     }
00159   };
00160 };
00161 
00173 class LyonChannelDiff: public MarSystem
00174 {
00175 private:
00176 
00177   void myUpdate(MarControlPtr sender)
00178   {
00179 
00181     MarSystem::myUpdate(sender);
00182 
00183     numBands_   = getctrl ("mrs_natural/numBands")->to<mrs_natural>();
00184 
00185     //   allocate memory
00186     procBuf1_.create(numBands_-1,1);
00187     procBuf2_.create(numBands_-1,1);
00188   };
00189 
00190   mrs_realvec   procBuf1_,
00191               procBuf2_;
00192   mrs_natural   numBands_;
00193 public:
00194   LyonChannelDiff(std::string name):MarSystem("LyonChannelDiff", name)
00195   {
00196     addControls ();
00197   };
00198   ~LyonChannelDiff() {};
00199 
00200   void addControls()
00201   {
00202     addctrl("mrs_natural/numBands", 1);
00203     setctrlState("mrs_natural/numBands", true);
00204   }
00205 
00206   MarSystem* clone() const
00207   {
00208     return new LyonChannelDiff(*this);
00209   };
00210 
00211   void myProcess(realvec& in, realvec& out)
00212   {
00213 
00214     mrs_natural t;
00215     for (t = 0; t < inSamples_; t++)
00216     {
00217       in.getSubMatrix (0,t, procBuf1_);
00218       in.getSubMatrix (1,t, procBuf2_);
00219       procBuf1_ -=procBuf2_;            // calc the difference over the filter channels
00220       out.setSubMatrix (1,t,procBuf1_);
00221       out(0,t)  = in(0,t);              // leave the first element the same
00222     }
00223   };
00224 };
00225 
00237 class LyonZeroOutPreEmph: public MarSystem
00238 {
00239 private:
00240 
00241   void myUpdate(MarControlPtr sender)
00242   {
00243 
00245     MarSystem::myUpdate(sender);
00246   };
00247 
00248 public:
00249   LyonZeroOutPreEmph(std::string name):MarSystem("LyonZeroOutPreEmph", name)
00250   {
00251   };
00252   ~LyonZeroOutPreEmph() {};
00253 
00254   MarSystem* clone() const
00255   {
00256     return new LyonZeroOutPreEmph(*this);
00257   };
00258 
00259   void myProcess(realvec& in, realvec& out)
00260   {
00261     mrs_natural t,o;
00262     for (t = 0; t < inSamples_; t++)
00263     {
00264       out(0,t)  = 0;
00265       out(1,t)  = 0;
00266       for (o = 2; o < onObservations_; o++)
00267         out(o,t)    = in(o,t);
00268     }
00269   };
00270 };
00271 }//namespace Marsyas
00272 
00273 
00274 
00277 
00278 LyonPassiveEar::LyonPassiveEar(mrs_string name)
00279   : MarSystem("LyonPassiveEar", name),
00280     fs_ (0),
00281     currDecimState_(0),
00282     numFilterChannels_ (0),
00283     passiveEar_ (0)
00284 {
00285   addControls();
00286 }
00287 
00288 
00289 LyonPassiveEar::~LyonPassiveEar()
00290 {
00291   if (passiveEar_)
00292   {
00293     delete passiveEar_;
00294     passiveEar_ = 0;
00295   }
00296 }
00297 
00298 MarSystem*
00299 LyonPassiveEar::clone() const
00300 {
00301   LyonPassiveEar    *newEar = new LyonPassiveEar(*this);
00302   if (passiveEar_)
00303     newEar->passiveEar_ = (Series*)passiveEar_->clone ();
00304 
00305   return newEar;
00306 }
00307 
00308 
00309 void
00310 LyonPassiveEar::addControls()
00311 {
00312   addctrl("mrs_natural/decimFactor", 1);
00313   addctrl("mrs_real/earQ", 8.0F);
00314   addctrl("mrs_real/stepFactor", 0.25F);
00315   addctrl("mrs_bool/channelDiffActive", true);
00316   addctrl("mrs_bool/agcActive", true);
00317   addctrl("mrs_real/decimTauFactor", 3.0F);
00318 
00319   addctrl("mrs_realvec/centerFreqs", centerFreqs_);
00320 
00321   setctrlState("mrs_natural/decimFactor", true);
00322   setctrlState("mrs_real/earQ", true);
00323   setctrlState("mrs_real/stepFactor", true);
00324   setctrlState("mrs_bool/channelDiffActive", true);
00325   setctrlState("mrs_bool/agcActive", true);
00326   setctrlState("mrs_real/decimTauFactor", true);
00327 
00328   setctrlState("mrs_realvec/centerFreqs", false);
00329 }
00330 
00331 void
00332 LyonPassiveEar::myUpdate(MarControlPtr sender)
00333 {
00334   (void) sender;  //suppress warning of unused parameter(s)
00335 
00336   if (!setParametersIntern ())
00337   {
00338     this->updateControlsIntern ();
00339     return;
00340   }
00341 
00342   MRSDIAG("LyonPassiveEar.cpp - LyonPassiveEar:myUpdate");
00343 
00344   //FilterBank creation
00345   // variable names more or less directly from matlab implementation...
00346   const mrs_real    Eb              = 1000.0;
00347   const mrs_real    EarZeroOffset   = 1.5;
00348   const mrs_real    EarSharpness    = 5.0;
00349   const mrs_real    EarPremphCorner = 300.0;
00350   mrs_natural       i,numChans;
00351   ostringstream name;
00352   mrs_realvec       a(3),
00353                 b(3);       
00354   mrs_real      lowFreq,
00355               topFreq           = .5 * getctrl ("mrs_real/israte")->to<mrs_real>();
00356   Cascade           *filterBank     = 0;
00357 
00358   //% Find top frequency, allowing space for first cascade filter.
00359   //    topf = fs/2.0;
00360   //  topf = topf - (sqrt(topf^2+Eb^2)/EarQ*StepFactor*EarZeroOffset)+ ...
00361   //    sqrt(topf^2+Eb^2)/EarQ*StepFactor;
00362   topFreq   = topFreq -
00363             (sqrt (sqr(topFreq) + sqr(Eb)) / earQ_ * stepFactor_ * EarZeroOffset) +
00364             sqrt (sqr(topFreq) + sqr(Eb)) / earQ_ * stepFactor_;
00365 
00366   //% Find place where CascadePoleQ < .5
00367   //    lowf = Eb/sqrt(4*EarQ^2-1);
00368   //  NumberOfChannels = floor((EarQ*(-log(lowf + sqrt(lowf^2 + Eb^2)) + ...
00369   //    log(topf + sqrt(Eb^2 + topf^2))))/StepFactor);
00370   lowFreq       = Eb / sqrt (4 * sqr(earQ_) - 1);
00371   numChans  = (mrs_natural)(floor ((earQ_ * (-log (lowFreq + sqrt (sqr(lowFreq) + sqr(Eb))) +
00372                                     log (topFreq + sqrt (sqr(Eb) + sqr(topFreq)))))/stepFactor_) + .1); // add .1 to ensure correct cast...
00373 
00374   //% Now make an array of CenterFreqs..... This expression was derived by
00375   //    % Mathematica by integrating 1/EarBandwidth(cf) and solving for f as a
00376   //    % function of filterctrl number.
00377   //    cn = 1:NumberOfChannels;
00378   //  CenterFreqs = (-((exp((cn*StepFactor)/EarQ)*Eb^2)/ ...
00379   //    (topf + sqrt(Eb^2 + topf^2))) + ...
00380   //    (topf + sqrt(Eb^2 + topf^2))./exp((cn*StepFactor)/EarQ))/2;
00381   centerFreqs_.create (numChans);
00382   for (i = 0; i < numChans; ++i)
00383     centerFreqs_(i) = (-((exp(((i+1)*stepFactor_)/earQ_) * sqr(Eb))/
00384                          (topFreq + sqrt(sqr(Eb) + sqr(topFreq)))) + (topFreq + sqrt(sqr(Eb) + sqr(topFreq)))/ exp(((i+1)*stepFactor_)/earQ_))*.5;
00385 
00386   // free memory if necessary...
00387   if (passiveEar_)
00388   {
00389     delete passiveEar_;
00390     passiveEar_ = 0;
00391   }
00392   // now create the processing chain
00393   passiveEar_ = new Series ("LyonPassiveEar");
00394   passiveEar_->addMarSystem (filterBank = new Cascade("LyonFilterBank"));
00395   passiveEar_->addMarSystem (new HalfWaveRectifier("hwr1"));
00396   passiveEar_->addMarSystem (new LyonZeroOutPreEmph("ZeroOut"));
00397   if (agcActive_)
00398     passiveEar_->addMarSystem (new LyonAgc("agc"));
00399   if (channelDiffActive_)
00400   {
00401     passiveEar_->addMarSystem (new LyonChannelDiff("differ"));
00402     passiveEar_->addMarSystem (new HalfWaveRectifier("hwr2"));
00403   }
00404 
00405   if (decimFactor_ > 1)
00406   {
00407     b(0)    = b(1) = 0;
00408     b(2)    = a(0) = 1;
00409     a(1)    = -2.0*(1-LyonAgc::lyonEpsilonFromTauFS (decimFactor_/fs_*decimTauFactor_, fs_));
00410     a(2)    = sqr(a(1)/(-2.0));
00411     b       *= lyonSetGain (b, a, 1.0, 0, fs_);
00412   }
00413   else
00414   {
00415     // otherwise set the filter to "bypass"
00416     b(0)    = a(0)  = 1;
00417     b(1)    = b(2)  = a(1)  = a(2)  = 0;
00418   }
00419 
00420   passiveEar_->addMarSystem (lyonCreateFilter (b,a,"decim"));
00421 
00422   //% Finally, let's design the front filters.
00423   //    front(1,:) = SetGain([0 1 -exp(-2*pi*EarPremphCorner/fs) 0 0], 1, fs/4, fs);
00424   //  topPoles = SecondOrderFilter(topf,CascadePoleQ(1), fs);
00425   //  front(2,:) = SetGain([1 0 -1 topPoles(2:3)], 1, fs/4, fs);
00426   b(0)  = a(1) = a(2) = 0;
00427   b(1)  = a(0) = 1;
00428   b(2)  = -exp(-TWOPI*EarPremphCorner/fs_);
00429   b     *= lyonSetGain (b, a, 1.0, fs_*.25, fs_);
00430 
00431   name.str("");
00432   name << "front_" << 0;
00433   filterBank->addMarSystem(lyonCreateFilter (b, a, name.str()));        // preemphasis 1
00434 
00435   b(0)  = 1;
00436   b(1)  = 0;
00437   b(2)  = -1;
00438   a     = lyonSecondOrderFilter (topFreq, centerFreqs_(0)/(sqrt(sqr(centerFreqs_(0)) + sqr(Eb))/earQ_), fs_);
00439   b     *= lyonSetGain (b, a, 1.0F, fs_*.25F, fs_);
00440 
00441   name.str("");
00442   name << "front_" << 1;
00443   filterBank->addMarSystem(lyonCreateFilter (b, a, name.str()));        // preemphasis 2
00444 
00445   for (i = 0; i < numChans; ++i)
00446   {
00447     mrs_real EarBandwidth, CascadeZeroCF, CascadeZeroQ, CascadePoleCF, CascadePoleQ;
00448 
00449     //% OK, now we can figure out the parameters of each stage filter.
00450     //  EarBandwidth = sqrt(CenterFreqs.^2+Eb^2)/EarQ;
00451     //  CascadeZeroCF = CenterFreqs +   EarBandwidth * StepFactor * EarZeroOffset;
00452     //  CascadeZeroQ = EarSharpness*CascadeZeroCF./EarBandwidth;
00453     //  CascadePoleCF = CenterFreqs;
00454     //  CascadePoleQ = CenterFreqs./EarBandwidth;
00455     EarBandwidth    = sqrt(sqr(centerFreqs_(i)) + sqr(Eb))/earQ_;
00456     CascadeZeroCF   = centerFreqs_(i) + EarBandwidth * stepFactor_ * EarZeroOffset;
00457     CascadeZeroQ    = EarSharpness * CascadeZeroCF / EarBandwidth;
00458     CascadePoleCF   = centerFreqs_(i);
00459     CascadePoleQ    = centerFreqs_(i)/EarBandwidth;
00460 
00461     // compute filter coeffs
00462     b   = lyonSecondOrderFilter (CascadeZeroCF, CascadeZeroQ, fs_);
00463     a   = lyonSecondOrderFilter (CascadePoleCF, CascadePoleQ, fs_);
00464 
00465     //% Now we can set the DC gain of each stage.
00466     //  dcgain(2:NumberOfChannels)=CenterFreqs(1:NumberOfChannels-1)./ ...
00467     //  CenterFreqs(2:NumberOfChannels);
00468     //  dcgain(1) = dcgain(2);
00469     //  for i=1:NumberOfChannels
00470     //  filters(i,:) = SetGain(filters(i,:), dcgain(i), 0, fs);
00471     //  end
00472     b   *= (i == 0) ? lyonSetGain (b, a, centerFreqs_(i)/centerFreqs_(i+1), 0, fs_) :
00473          lyonSetGain (b, a, centerFreqs_(i-1)/centerFreqs_(i), 0, fs_);
00474 
00475     // create the filter and add it
00476     name.str("");
00477     name << "filter_" << i;
00478     filterBank->addMarSystem(lyonCreateFilter (b, a, name.str()));
00479   }
00480 
00481   numFilterChannels_    = numChans + 2; // plus two front filters...
00482 
00483   this->updateControlsIntern ();
00484 
00485   // alloc some internal memory
00486   tmpOut_.create (numFilterChannels_, inSamples_);
00487   decimOut_.create (numFilterChannels_-2, inSamples_/decimFactor_); // integer division intentionally;
00488 }
00489 
00490 void
00491 LyonPassiveEar::myProcess(realvec& in, realvec& out)
00492 {
00493   if(getctrl("mrs_bool/mute")->to<mrs_bool>()) return;
00494 
00495   // check number of out observations
00496   mrs_natural   i,
00497               currCount     = -currDecimState_-1,
00498                numOutSamples    = (inSamples_+currDecimState_)/decimFactor_;// integer division intended
00499 
00500   MRSASSERT(currDecimState_ <= decimFactor_);
00501 
00502   // update the number of output samples if necessary
00503   if (onSamples_ != numOutSamples)
00504     updControl ("mrs_natural/onSamples", numOutSamples);
00505 
00506   decimOut_.stretch (numFilterChannels_-2, numOutSamples);
00507 
00508   // process the series
00509   passiveEar_->process(in, tmpOut_);
00510 
00511   // post-"processing": do the sample rate decimation and remove the pre-emphasis filters
00512   for (i = 0; i < numOutSamples; ++i)
00513   {
00514     mrs_realvec decimTmp(numFilterChannels_-2,1);
00515     currCount   = currCount + decimFactor_;
00516     MRSASSERT(currCount <= inSamples_);
00517     tmpOut_.getSubMatrix (2, currCount, decimTmp);
00518     decimOut_.setSubMatrix (0, i, decimTmp);
00519   }
00520   currDecimState_   = inSamples_ - currCount-1;
00521   out   = decimOut_;
00522 }
00523 
00524 
00526 mrs_bool LyonPassiveEar::setParametersIntern ()
00527 {
00528   mrs_bool  updateMe            = false;
00529 
00530   updateMe  |= (passiveEar_ == 0);
00531 
00532   // update controls
00533   if (decimFactor_      != getctrl ("mrs_natural/decimFactor")->to<mrs_natural>())
00534   {
00535     updateMe            = true;
00536     decimFactor_        = getctrl ("mrs_natural/decimFactor")->to<mrs_natural>();
00537   }
00538   if (earQ_             != getctrl ("mrs_real/earQ")->to<mrs_real>())
00539   {
00540     updateMe            = true;
00541     earQ_               = getctrl ("mrs_real/earQ")->to<mrs_real>();
00542   }
00543   if (stepFactor_           != getctrl ("mrs_real/stepFactor")->to<mrs_real>())
00544   {
00545     updateMe            = true;
00546     stepFactor_         = getctrl ("mrs_real/stepFactor")->to<mrs_real>();
00547   }
00548   if (channelDiffActive_    != getctrl ("mrs_bool/channelDiffActive")->to<mrs_bool>())
00549   {
00550     updateMe            = true;
00551     channelDiffActive_  = getctrl ("mrs_bool/channelDiffActive")->to<mrs_bool>();
00552   }
00553   if (agcActive_            != getctrl ("mrs_bool/agcActive")->to<mrs_bool>())
00554   {
00555     updateMe            = true;
00556     agcActive_          = getctrl ("mrs_bool/agcActive")->to<mrs_bool>();
00557   }
00558   if (decimTauFactor_       != getctrl ("mrs_real/decimTauFactor")->to<mrs_real>())
00559   {
00560     updateMe            = true;
00561     decimTauFactor_     = getctrl ("mrs_real/decimTauFactor")->to<mrs_real>();
00562   }
00563   if (fs_                   != getctrl ("mrs_real/israte")->to<mrs_real>())
00564   {
00565     updateMe            = true;
00566     fs_                 = getctrl ("mrs_real/israte")->to<mrs_real>();
00567   }
00568 
00569   return updateMe;
00570 }
00571 
00572 void LyonPassiveEar::updateControlsIntern ()
00573 {
00574   passiveEar_->updControl("mrs_natural/inObservations", getctrl("mrs_natural/inObservations")->to<mrs_natural>());
00575   passiveEar_->updControl("mrs_natural/inSamples", getctrl("mrs_natural/inSamples")->to<mrs_natural>());
00576   passiveEar_->updControl("mrs_real/israte", getctrl("mrs_real/israte")->to<mrs_real>());
00577 
00578   ctrl_onSamples_->setValue(inSamples_/decimFactor_); // integer division intended
00579   ctrl_osrate_->setValue(israte_*1.0/decimFactor_);
00580 
00581   if (numFilterChannels_)
00582   {
00583     updControl ("mrs_realvec/centerFreqs", centerFreqs_);
00584     ctrl_onObservations_->setValue((numFilterChannels_-2)*getctrl("mrs_natural/inObservations")->to<mrs_natural>());
00585 
00586     passiveEar_->setctrl("mrs_natural/onObservations", getctrl("mrs_natural/onObservations")->to<mrs_natural>());
00587     if (agcActive_)
00588       passiveEar_->updControl("LyonAgc/agc/mrs_natural/numBands", numFilterChannels_);
00589     if (channelDiffActive_)
00590       passiveEar_->updControl("LyonChannelDiff/differ/mrs_natural/numBands", numFilterChannels_);
00591   }
00592 }
00593 
00595 mrs_realvec LyonPassiveEar::lyonSecondOrderFilter (mrs_real midFreq, mrs_real q, mrs_real sRate)
00596 {
00597   //function filts = SecondOrderFilter(f,q,fs)
00598   //cft = f'/fs;
00599   //rho = exp(- pi * cft ./ q');
00600   //          theta = 2 * pi * cft .* sqrt(1-1 ./(4*q'.^2));
00601   //          filts = [ ones(size(rho)) -2*rho.*cos(theta) rho.^2 ];
00602   mrs_realvec   result (3);
00603   mrs_real cft  = midFreq / sRate,
00604             rho     = exp (-PI * cft / q);
00605 
00606   result(0) = 1;
00607   result(1) = -2 * rho * cos (TWOPI * cft * sqrt (1 - 1.0/(4*sqr(q))));
00608   result(2) = sqr(rho);
00609   return result;
00610 }
00611 
00613 mrs_real LyonPassiveEar::lyonFreqResp (mrs_realvec firCoeffs, mrs_realvec iirCoeffs, mrs_real freq, mrs_real sRate, mrs_bool inDb)
00614 {
00615   //function mag=FreqResp(filter,f,fs)
00616   //    cf = exp(i*2*pi*f/fs);
00617   //mag = (filter(3) + filter(2)*cf + filter(1)*cf.^2) ./ ...
00618   //    (filter(5) + filter(4)*cf + cf.^2);
00619   //mag = 20*log10(abs(mag));
00620   mrs_complex cf    = mrs_complex (cos (TWOPI * freq / sRate), sin (TWOPI * freq / sRate));
00621   mrs_complex mag   = (firCoeffs(2) + firCoeffs(1) * cf + firCoeffs(0) * (cf*cf)) / (iirCoeffs(2) + iirCoeffs(1) * cf + (cf*cf));
00622   mrs_real  res = sqrt (sqr(mag.real ()) + sqr (mag.imag ()));
00623 
00624   return inDb? 20.0/log(10.0) * log (res) : res;
00625 }
00626 
00628 mrs_real LyonPassiveEar::lyonSetGain (mrs_realvec firCoeffs, mrs_realvec iirCoeffs, mrs_real newGain, mrs_real freq, mrs_real sRate)
00629 {
00630   //function filter = SetGain(filter, desired, f, fs)
00631   //    oldGain = 10^(FreqResp(filter, f, fs)/20);
00632   //filter(1:3) = filter(1:3)*desired/oldGain;
00633 
00634   return newGain / lyonFreqResp (firCoeffs, iirCoeffs, freq, sRate, false);
00635 }
00636 
00638 Filter*     LyonPassiveEar::lyonCreateFilter (mrs_realvec b, mrs_realvec a, mrs_string name)
00639 {
00640   Filter *filter = new Filter(name);
00641   filter->setctrl("mrs_realvec/ncoeffs", b);
00642   filter->setctrl("mrs_realvec/dcoeffs", a);
00643 
00644   return filter;
00645 }