Marsyas
0.6.0-alpha
|
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 }