Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/SimulMaskingFft.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2010 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 "../common_source.h"
00020 #include "SimulMaskingFft.h"
00021 #include <algorithm>
00022 //#define MTLB_DBG_LOG
00023 
00024 using std::min;
00025 using std::max;
00026 using namespace Marsyas;
00027 
00028 static const mrs_natural h2bIdx = 3;
00029 static const mrs_real   lowFreq = 80.;
00030 static const mrs_real   upFreq  = 18000.;
00031 
00032 static const mrs_real   truncTresh = 10;
00033 
00034 static inline mrs_real IntPow (mrs_real a, mrs_natural b)
00035 {
00036   if (b == 0)
00037     return 1.;
00038   mrs_real  dResult = a;
00039   while (--b > 0)
00040     dResult *= a;
00041   return (dResult < 1e-30)? 0 : dResult;
00042 }
00043 
00044 SimulMaskingFft::SimulMaskingFft(mrs_string name):MarSystem("SimulMaskingFft", name)
00045 {
00046   //Add any specific controls needed by SimulMaskingFft
00047   //(default controls all MarSystems should have
00048   //were already added by MarSystem::addControl(),
00049   //called by :MarSystem(name) constructor).
00050   //If no specific controls are needed by a MarSystem
00051   //there is no need to implement and call this addControl()
00052   //method (see for e.g. Rms.cpp)
00053   addControls();
00054 
00055   numBands_ = 0;
00056   freqBounds_   = 0;
00057   numBands_ = 0;
00058 
00059 }
00060 
00061 SimulMaskingFft::SimulMaskingFft(const SimulMaskingFft& a) : MarSystem(a)
00062 {
00063   // For any MarControlPtr in a MarSystem
00064   // it is necessary to perform this getctrl
00065   // in the copy constructor in order for cloning to work
00066   ctrl_SimulMaskingFft_ = getctrl("mrs_real/SimulMaskingFft");
00067 }
00068 
00069 SimulMaskingFft::~SimulMaskingFft()
00070 {
00071   if (freqBounds_)
00072     delete freqBounds_;
00073   freqBounds_   = 0;
00074 }
00075 
00076 MarSystem*
00077 SimulMaskingFft::clone() const
00078 {
00079   SimulMaskingFft *New  = new SimulMaskingFft(*this);
00080 
00081   if (this->numBands_ > 0)
00082   {
00083     New->freqBounds_    = new FrequencyBands_t [numBands_];
00084     New->ComputeTables ();
00085   }
00086   else
00087     New->freqBounds_    = 0;
00088 
00089 
00090   return New;
00091 }
00092 
00093 void
00094 SimulMaskingFft::addControls()
00095 {
00096   //Add specific controls needed by this MarSystem.
00097   addctrl("mrs_real/listeningLevelInDbSpl", 92.0);
00098   setctrlState("mrs_real/listeningLevelInDbSpl", true);
00099 }
00100 
00101 void
00102 SimulMaskingFft::myUpdate(MarControlPtr sender)
00103 {
00104   // no need to do anything SimulMaskingFft-specific in myUpdate
00105   MarSystem::myUpdate(sender);
00106 
00107 
00108   // compute audio samplerate, numBands, and normalization
00109   audiosrate_   = israte_*(mrs_real)(inObservations_-1)*2;
00110   barkRes_  = hertz2bark (lowFreq+israte_, h2bIdx)-hertz2bark (lowFreq,h2bIdx); // Delta f / N
00111   numBands_ = (mrs_natural)((hertz2bark (upFreq, h2bIdx) - hertz2bark (lowFreq, h2bIdx) + .5)/barkRes_);
00112   normFactor_   = (1<<15)*sqrt(8./3.)*2*20e-6*pow (10.,.05*getctrl("mrs_real/listeningLevelInDbSpl")->to<mrs_real>());
00113 
00114   // alloc memory
00115   if (numBands_ <= 0)
00116     return;
00117   processBuff_.stretch(inObservations_);
00118   processBuff_.setval (0);
00119   helpBuff_.stretch(inObservations_);
00120   helpBuff_.setval (0);
00121   outerEar_.stretch(inObservations_);
00122   outerEar_.setval (0);
00123   barkSpec_.stretch(numBands_);
00124   barkSpec_.setval (0);
00125   excPattern_.stretch(numBands_);
00126   excPattern_.setval (0);
00127   maskingThresh_.stretch(numBands_);
00128   maskingThresh_.setval (0);
00129   intNoise_.stretch(numBands_);
00130   intNoise_.setval (0);
00131   slopeSpread_.stretch(numBands_);
00132   slopeSpread_.setval (0);
00133   normSpread_.stretch(numBands_);
00134   normSpread_.setval (0);
00135 
00136   if (freqBounds_)
00137     delete freqBounds_;
00138   freqBounds_   = new FrequencyBands_t [numBands_];;
00139 
00140   ComputeTables ();
00141 }
00142 
00143 
00144 void
00145 SimulMaskingFft::myProcess(realvec& in, realvec& out)
00146 {
00147   for (mrs_natural t = 0; t < inSamples_; t++)
00148   {
00149     in.getCol(t, processBuff_);
00150 
00151     // normalize spectrum
00152     processBuff_    *= normFactor_;
00153     processBuff_    *= processBuff_;
00154 
00155     // weight with outer ear transfer function
00156     processBuff_    *= outerEar_;
00157 
00158     // compute bark spectrum
00159     GetBandLevels (freqBounds_, barkSpec_, false);
00160 
00161     // add internal noise
00162     barkSpec_   += intNoise_;
00163 
00164     // compute spreading function
00165     CalcSpreading (barkSpec_, excPattern_);
00166 
00167     // apply  masking threshold
00168     excPattern_ *= maskingThresh_;
00169 
00170     // normalize input spectrum
00171     in.getCol(t, processBuff_);
00172     processBuff_    *= normFactor_;
00173     processBuff_    *= processBuff_;
00174 
00175     // compute difference
00176     ComputeDifference (out, processBuff_,  t);
00177 
00178   }
00179 }
00180 
00181 
00182 void
00183 SimulMaskingFft::ComputeDifference (mrs_realvec &out, mrs_realvec &in, mrs_natural t)
00184 {
00185   mrs_natural i;
00186   t = 0;
00187 
00188   for (i = 0; i < inObservations_; ++i)
00189     out(i,t) = 0;
00190 
00191   for (mrs_natural k = 0; k < numBands_; k++)
00192   {
00193     mrs_real   fLowFrac    = freqBounds_[k].fLowFreqBound/audiosrate_ * (inObservations_<<1),
00194                fHighFrac   = freqBounds_[k].fUpFreqBound/audiosrate_ *(inObservations_<<1);
00195     mrs_natural     iLowBin     = (mrs_natural)ceil (fLowFrac),
00196                     iHighBin    = (mrs_natural)floor(fHighFrac);
00197     for (i = iLowBin; i <= iHighBin; ++i)
00198     {
00199       if (excPattern_(k) <= 1./truncTresh*in(i))
00200         out(i,t)    = truncTresh;
00201       else if (excPattern_(k) >= truncTresh*in(i))
00202         out(i,t)    = 1./truncTresh;
00203       else
00204         out(i,t)    = in(i) / excPattern_(k);
00205     }
00206   }
00207 #ifdef MARSYAS_MATLAB
00208 #ifdef MTLB_DBG_LOG
00209   MATLAB_PUT(in, "a");
00210   MATLAB_PUT(out, "out");
00211   MATLAB_PUT(audiosrate_, "fs");
00212   MATLAB_PUT(normFactor_, "normFactor");
00213   MATLAB_EVAL("figure(1);clf ;semilogy((1:length(out))*fs/(2*length(out)),sqrt(a)/normFactor, (1:length(out))*fs/(2*length(out)),sqrt(a./out)/normFactor,'g'),grid on, axis([50 4000 1e-4 1])");
00214   //MATLAB_PUT(excPattern_, "exc");
00215   //MATLAB_EVAL("figure(2);clf ;plot(exc,'r'),grid on");
00216 #endif
00217 #endif
00218 }
00219 
00220 void
00221 SimulMaskingFft::GetBandLevels (FrequencyBands_t *pFrequencyValues, mrs_realvec &bandLevels, mrs_bool bDezibel)
00222 {
00223 
00224   for (mrs_natural i = 0; i < numBands_; ++i)
00225   {
00226     mrs_real   fLowFrac    = pFrequencyValues[i].fLowFreqBound/audiosrate_ * (inObservations_<<1),
00227                fHighFrac   = pFrequencyValues[i].fUpFreqBound/audiosrate_ *(inObservations_<<1);
00228     mrs_natural     iLowBin     = (mrs_natural)ceil (fLowFrac),
00229                     iHighBin    = (mrs_natural)floor(fHighFrac);
00230 
00231     fLowFrac            = iLowBin - fLowFrac;
00232     fHighFrac           = fHighFrac - iHighBin;
00233     bandLevels(i)   = fLowFrac * processBuff_(max(0L,iLowBin-1));
00234     bandLevels(i)  += fHighFrac * processBuff_(min((mrs_natural)(inObservations_ - .5),iHighBin+1));
00235     for (mrs_natural j = iLowBin; j < iHighBin; j++)
00236       bandLevels(i)  += processBuff_(j);
00237     if (bDezibel)
00238     {
00239       bandLevels(i)   = max ((mrs_real)bandLevels(i), (mrs_real)1e-20);
00240       bandLevels(i)   = 10./log(10.) * log ((bandLevels(i)));
00241     }
00242   }
00243 
00244   return;
00245 }
00246 
00247 void
00248 SimulMaskingFft::CalcSpreading (mrs_realvec &bandLevels, mrs_realvec &result)
00249 {
00250   // this is level dependent adapted from ITU-R BS.1387
00251 
00252   mrs_natural    iBarkj,                 // Masker
00253                  iBarkk;                 // Maskee
00254 
00255   mrs_real  fTmp1,
00256             fTmp2,
00257             fSlope,
00258             fScale      = sqrt(8./3.),
00259             fBRes       = barkRes_,//hertz2bark (.5*audiosrate_, h2bIdx)/numBands_,
00260             *pfEnPowTmp = processBuff_.getData (),
00261              *pfSlopeUp  = helpBuff_.getData ();
00262   mrs_real  *pfSlope    = slopeSpread_.getData (),
00263              *pfNorm     = normSpread_.getData ();
00264 
00265   // initialize pfResult
00266   result.setval(0);
00267 
00268   fSlope                      = exp ( -fBRes * 2.7 * 2.302585092994045684017991454684364207601101488628772976033);
00269   fTmp2                     = 1.0 / (1.0 - fSlope);
00270   for (iBarkj = 0; iBarkj < numBands_; iBarkj++)
00271   {
00272     pfSlopeUp[iBarkj]       = pfSlope[iBarkj] * pow (bandLevels(iBarkj)*fScale,  .2 * fBRes);
00273     fTmp1                   = (1.0 - IntPow (fSlope, iBarkj+1)) * fTmp2;
00274     fTmp2                   = (1.0 - IntPow(pfSlopeUp[iBarkj], numBands_ - iBarkj)) / (1.0 - pfSlopeUp[iBarkj]);
00275     if (bandLevels(iBarkj) < 1e-20)
00276     {
00277       pfSlopeUp[iBarkj]   = 0;
00278       pfEnPowTmp[iBarkj]  = 0;
00279       continue;
00280     }
00281     pfSlopeUp[iBarkj]       = exp (0.4 * log (pfSlopeUp[iBarkj]));
00282     pfEnPowTmp[iBarkj]      = exp (0.4 * log (bandLevels(iBarkj)/(fTmp1 + fTmp2 -1)));
00283   }
00284   fSlope                      = exp ( 0.4 * log (fSlope));
00285 
00286   // lower slope
00287   result(numBands_-1)     = pfEnPowTmp[numBands_-1];
00288   for (iBarkk = numBands_-2; iBarkk >= 0; iBarkk--)
00289     result(iBarkk)        = pfEnPowTmp[iBarkk] + result(iBarkk + 1) * fSlope;
00290 
00291   // upper slope
00292   for (iBarkj = 0; iBarkj < numBands_-1; iBarkj++)
00293   {
00294     fSlope                  = pfSlopeUp[iBarkj];
00295     fTmp1                   = pfEnPowTmp[iBarkj];
00296     for (iBarkk = iBarkj+1; iBarkk < numBands_; iBarkk++)
00297     {
00298       mrs_real  dTmp1       = fTmp1 * fSlope;
00299       fTmp1               = (dTmp1 < 1e-30)? 0 : (mrs_real)dTmp1;
00300       result(iBarkk)   += fTmp1;
00301     }
00302   }
00303 
00304   // normalization
00305   for (iBarkk = 0; iBarkk < numBands_; iBarkk++)
00306   {
00307     mrs_real  dTmp      = result(iBarkk);
00308     result(iBarkk)      = sqrt(dTmp) * dTmp * dTmp *pfNorm[iBarkk];
00309     //result(iBarkk)        = sqrt (result(iBarkk));
00310     //result(iBarkk)       *= result(iBarkk)*result(iBarkk)*result(iBarkk)*result(iBarkk);
00311     //result(iBarkk)       *= pfNorm[iBarkk];
00312   }
00313   return;
00314 }
00315 
00316 void
00317 SimulMaskingFft::ComputeTables ()
00318 {
00319   mrs_natural i;
00320 
00321 
00322   // outer ear transfer function
00323   {
00324     for (i = 0; i < inObservations_; ++i)
00325     {
00326       mrs_real dTmp;
00327       mrs_real fkFreq    = i * .5e-3 * audiosrate_ / inObservations_;
00328       if (fkFreq < 1e-10)
00329       {
00330         outerEar_(i)   = 0;
00331         continue;
00332       }
00333       dTmp  = -.2184 * pow (fkFreq, -.8);
00334       dTmp   += .65 * exp (-.6 * (fkFreq-3.3)*(fkFreq-3.3));
00335       //outerEar_(i)   = (dTmp > 1e-37)? (mrs_real)dTmp : 0;
00336       outerEar_(i)  = dTmp - 1e-4 * pow (fkFreq, 3.6);
00337       if (outerEar_(i) < -12)
00338         outerEar_(i) = 0;
00339       else
00340         outerEar_(i)   = pow (10.,outerEar_(i));
00341     }
00342   }
00343 
00344   // frequency bands and spreading
00345   {
00346     mrs_real fLowBark = hertz2bark (lowFreq, h2bIdx);
00347     mrs_real fMaxBark = hertz2bark (.5*audiosrate_-1, h2bIdx);
00348     for (i = 0; i < numBands_; ++i)
00349     {
00350       freqBounds_[i].fLowBarkBound  = min(fMaxBark,fLowBark + i*barkRes_);
00351       freqBounds_[i].fMidBark       = min((mrs_real)fMaxBark,(mrs_real)freqBounds_[i].fLowBarkBound + (mrs_real).5*barkRes_);
00352       freqBounds_[i].fUpBarkBound   = min(fMaxBark,freqBounds_[i].fLowBarkBound + barkRes_);
00353 
00354       freqBounds_[i].fLowFreqBound  = bark2hertz (freqBounds_[i].fLowBarkBound, h2bIdx);
00355       freqBounds_[i].fMidFreq       = bark2hertz (freqBounds_[i].fMidBark, h2bIdx);
00356       freqBounds_[i].fUpFreqBound   = bark2hertz (freqBounds_[i].fUpBarkBound, h2bIdx);
00357     }
00358 
00359     for (i = 0; i < numBands_; ++i)
00360     {
00361       mrs_natural     cBarkk;
00362       mrs_real   fAtt    = 1.0,
00363                  fNorm   = 1.0,
00364                  fSlope  = pow (10.0, -2.7 * barkRes_);
00365 
00366       slopeSpread_(i) = 24.0 + 230./freqBounds_[i].fMidFreq;
00367       slopeSpread_(i) = pow (10.,-0.1 * barkRes_ * slopeSpread_(i));
00368 
00369       processBuff_(i)    = 1.0;
00370       // lower slope
00371       for (cBarkk = i; cBarkk > 0; cBarkk--)
00372       {
00373         mrs_real  dTmp                = fAtt * fSlope;
00374         fAtt                        = (dTmp < 1e-30) ? 0 : (mrs_real) dTmp;
00375 
00376         processBuff_(cBarkk-1) = fAtt;                                             // nonnormalized masking threshold
00377         fNorm                      += fAtt;                                             // normalization factor
00378       }
00379 
00380 
00381       // initialize new
00382       fAtt                            = 1.0F;
00383 
00384       // upper slope
00385       for (cBarkk = i; cBarkk < numBands_-1; cBarkk++)
00386       {
00387         mrs_real  dTmp                = fAtt * slopeSpread_(i);
00388         fAtt                        = (dTmp < 1e-30) ? 0 : (mrs_real) dTmp;
00389 
00390         processBuff_(cBarkk+1) = fAtt;
00391         fNorm                      += fAtt;                                             // normalization factor
00392       }
00393 
00394       fNorm                           = 1.0/fNorm;
00395 
00396       // nonlinear superposition
00397       for (cBarkk = 0; cBarkk < numBands_; cBarkk++)
00398       {
00399         processBuff_(cBarkk)                  *= fNorm;
00400         normSpread_(cBarkk)    += pow (processBuff_(cBarkk), 0.4);
00401       }
00402     }
00403     for (i = 0; i < numBands_; ++i)
00404       normSpread_(i)  = pow (normSpread_(i), -2.5);            // resulting normalization pattern (eq. 19)
00405   }
00406 
00407   {
00408     // masking (eq.25)
00409     for ( i = 0; i <numBands_; ++i )
00410     {
00411       intNoise_(i)     = .1456 * pow (.001 * freqBounds_[i].fMidFreq, -0.8);    // noise in dB
00412       intNoise_(i)     = pow (10., intNoise_(i));                          // -> in energy domain
00413     }
00414   }
00415 
00416   {
00417     // masking (eq.25)
00418     mrs_real    v           = pow (.1, .3);
00419     for ( i = 0; i < 12.0/barkRes_; ++i )
00420       maskingThresh_(i)   = v;
00421 
00422 
00423     while (i < numBands_)
00424     {
00425       maskingThresh_(i)   = pow (.1, .025 * barkRes_ * i);
00426       ++i;
00427     }
00428   }
00429 
00430   return;
00431 }