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