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 #include <algorithm> 00019 00020 #include "../common_source.h" 00021 #include <marsyas/common_header.h> 00022 #include "PeakConvert2.h" 00023 #include "Peaker.h" 00024 #include "MaxArgMax.h" 00025 #include "SimulMaskingFft.h" 00026 #include <marsyas/peakView.h> 00027 00028 00029 //#define MTLB_DBG_LOG 00030 //#define ORIGINAL_VERSION 00031 //#define LOG2FILE 00032 00033 00034 #ifdef LOG2FILE 00035 #include <iomanip> 00036 static std::ofstream pFDbgFile; 00037 static const std::string kDbgFilePath = "d:/temp/peaks.new.txt"; 00038 #endif 00039 00040 00041 00042 using std::ostringstream; 00043 using std::min; 00044 using std::max; 00045 using std::abs; 00046 00047 00048 using namespace Marsyas; 00049 00050 static const mrs_real gaussianStd = 0.42466090014401; // results in output of .5 for input of .5 00051 00052 static void FreqSmear (mrs_realvec &spectrum) 00053 { 00054 mrs_natural length = spectrum.getSize (); 00055 mrs_real buf[3] = {0,0,0}, 00056 coeff[3] = {.25, .5, .25}; 00057 00058 spectrum(0) = spectrum(length-1) = 0; 00059 00060 for (mrs_natural i = 0; i < length-1; i++) 00061 { 00062 buf[(i+1)%3] = spectrum(i+1); 00063 spectrum(i) = coeff[0]*buf[(i-1+3)%3] + coeff[1]*buf[(i)%3] + coeff[2]*buf[(i+1)%3]; 00064 } 00065 00066 return; 00067 } 00068 mrs_real 00069 PeakConvert2::GaussianPdf (mrs_real x, mrs_real std) 00070 { 00071 return exp (-(x*x)/(2*(std*std)));// / sqrt (TWOPI*std); 00072 } 00073 00074 mrs_real princArg (mrs_real phase) 00075 { 00076 mrs_real fx = phase + PI; 00077 return PI + (fx + TWOPI*floor (fx*(-1.0/TWOPI))); 00078 } 00079 00080 PeakConvert2::PeakConvert2(mrs_string name):MarSystem("PeakConvert2",name), 00081 peaker_(0), 00082 max_(0), 00083 masking_(0) 00084 { 00085 psize_ = 0; 00086 size_ = 0; 00087 nbParameters_ = peakView::nbPkParameters; 00088 skip_=0; 00089 frame_ = 0; 00090 N_ = 0; 00091 00092 fundamental_ = 0.0; 00093 factor_ = 0.0; 00094 nbPeaks_ = 0; 00095 frameMaxNumPeaks_ = 0; 00096 instFreqHopSize_ = 1; 00097 00098 useStereoSpectrum_ = false; 00099 00100 peaker_ = new Peaker("Peaker"); 00101 max_ = new MaxArgMax("MaxArgMax"); 00102 masking_ = new SimulMaskingFft("masking"); 00103 00104 addControls(); 00105 } 00106 00107 PeakConvert2::PeakConvert2(const PeakConvert2& a) : MarSystem(a) 00108 { 00109 psize_ = a.psize_; 00110 size_ = a.size_; 00111 nbParameters_ = a.nbParameters_; 00112 skip_ = a.skip_; 00113 frame_ = a.frame_; 00114 N_ = a.N_; 00115 00116 fundamental_ = a.fundamental_; 00117 factor_ = a.factor_; 00118 nbPeaks_ = a.nbPeaks_; 00119 frameMaxNumPeaks_ = a.frameMaxNumPeaks_; 00120 hopSize_ = a.hopSize_; 00121 instFreqHopSize_ = 1; 00122 00123 useStereoSpectrum_ = a.useStereoSpectrum_; 00124 00125 peaker_ = (Peaker*)a.peaker_->clone (); 00126 max_ = (MaxArgMax*)a.max_->clone (); 00127 masking_ = (SimulMaskingFft*)a.masking_->clone (); 00128 00129 ctrl_totalNumPeaks_ = getctrl("mrs_natural/totalNumPeaks"); 00130 ctrl_frameMaxNumPeaks_ = getctrl("mrs_natural/frameMaxNumPeaks"); 00131 00132 #ifdef LOG2FILE 00133 pFDbgFile.open (kDbgFilePath.c_str (), std::ios::out); 00134 #endif 00135 } 00136 00137 PeakConvert2::~PeakConvert2() 00138 { 00139 delete peaker_; 00140 delete max_; 00141 if (masking_) 00142 delete masking_; 00143 #ifdef LOG2FILE 00144 pFDbgFile.close (); 00145 #endif 00146 } 00147 00148 MarSystem* 00149 PeakConvert2::clone() const 00150 { 00151 return new PeakConvert2(*this); 00152 } 00153 00154 void 00155 PeakConvert2::addControls() 00156 { 00157 realvec tmp(3); 00158 addctrl("mrs_natural/frameMaxNumPeaks", 0); 00159 setctrlState("mrs_natural/frameMaxNumPeaks", true); 00160 00161 addctrl("mrs_string/frequencyInterval", "MARSYAS_EMPTY"); 00162 setctrlState("mrs_string/frequencyInterval", true); 00163 00164 addctrl("mrs_natural/nbFramesSkipped", 0); 00165 setctrlState("mrs_natural/nbFramesSkipped", true); 00166 00167 addctrl("mrs_bool/improvedPrecision", true); 00168 setctrlState("mrs_bool/improvedPrecision", true); 00169 00170 addctrl("mrs_bool/picking", true); 00171 setctrlState("mrs_bool/picking", true); 00172 00173 addctrl("mrs_natural/hopSize", 1); 00174 setctrlState("mrs_natural/hopSize", true); 00175 00176 addctrl("mrs_real/probabilityTresh" , .5); 00177 setctrlState("mrs_real/probabilityTresh", true); 00178 00179 addctrl("mrs_natural/totalNumPeaks", 0, ctrl_totalNumPeaks_); 00180 00181 #ifdef ORIGINAL_VERSION 00182 addctrl("mrs_bool/useMasking", false); 00183 setctrlState("mrs_bool/useMasking", true); 00184 00185 tmp(0) = 0; // probability weight for peak being a masker 00186 tmp(1) = 0; // probability weight for peak being stationary 00187 tmp(2) = 1; // probability weight for peak being tonal 00188 addctrl("mrs_realvec/peakProbabilityWeight", tmp); 00189 setctrlState("mrs_realvec/peakProbabilityWeight", true); 00190 00191 addctrl( "mrs_real/peakSmearingTimeInS" , .0); // check with other hopsizes 00192 setctrlState( "mrs_real/peakSmearingTimeInS", true); 00193 #else 00194 addctrl("mrs_bool/useMasking", true); 00195 setctrlState("mrs_bool/useMasking", true); 00196 00197 00198 tmp(0) = 1; // probability weight for peak being a masker 00199 tmp(1) = 1; // probability weight for peak being stationary 00200 tmp(2) = 1; // probability weight for peak being tonal 00201 00202 addctrl("mrs_realvec/peakProbabilityWeight", tmp); 00203 setctrlState("mrs_realvec/peakProbabilityWeight", true); 00204 00205 addctrl( "mrs_real/peakSmearingTimeInS" , 0.03); // check with other hopsizes 00206 setctrlState( "mrs_real/peakSmearingTimeInS", true); 00207 #endif 00208 } 00209 00210 void 00211 PeakConvert2::myUpdate(MarControlPtr sender) 00212 { 00213 //(void) sender; //suppress warning of unused parameter(s) 00214 MarSystem::myUpdate (sender); 00215 00216 hopSize_ = getctrl ("mrs_natural/hopSize")->to<mrs_natural>(); 00217 mrs_real timeSrate = israte_*(mrs_real)N_;//israte_*(mrs_real)inObservations_/2.0; 00218 //check the input to see if we are also getting stereo information 00219 //(N_ is the FFT size) 00220 if (fmod(inObservations_, 2.0) == 0.0) 00221 { 00222 //we just have the two shifted spectra stacked vertically 00223 //(input has N + N observations) 00224 N_ = inObservations_/2; 00225 useStereoSpectrum_ = false; 00226 } 00227 else if(fmod(inObservations_-1, 2.5) == 0.0) 00228 { 00229 //we also have stereo spectrum info at the bottom 00230 //(input has N + N + N/2+1 observations) 00231 N_ = (mrs_natural)((inObservations_-1) / 2.5); 00232 useStereoSpectrum_ = true; 00233 } 00234 size_ = N_/2+1;//inObservations_ /4 +1; 00235 00236 00237 skip_ = getctrl("mrs_natural/nbFramesSkipped")->to<mrs_natural>(); 00238 prec_ = getctrl("mrs_bool/improvedPrecision")->to<mrs_bool>(); 00239 pick_ = getctrl("mrs_bool/picking")->to<mrs_bool>(); 00240 if(getctrl("mrs_string/frequencyInterval")->to<mrs_string>() != "MARSYAS_EMPTY") 00241 { 00242 realvec conv(2); 00243 string2parameters(getctrl("mrs_string/frequencyInterval")->to<mrs_string>(), conv, '_'); //[!] 00244 downFrequency_ = (mrs_natural) floor(conv(0)/timeSrate*size_*2) ; 00245 upFrequency_ = min(size_,(mrs_natural) floor(conv(1)/timeSrate*size_*2)); 00246 } 00247 else 00248 { 00249 downFrequency_ = 0; 00250 upFrequency_ = size_; 00251 } 00252 00253 00254 //if picking is disabled (==false), the number of sinusoids should be set 00255 //to the number of unique bins of the spectrums at the input (i.e. N/2+1) 00256 if(!pick_ /*&& ctrl_frameMaxNumPeaks_->to<mrs_natural>() == 0*/) 00257 { 00258 //frameMaxNumPeaks_ = N_/2+1; //inObservations_/4+1; 00259 //downFrequency_ = 0; 00260 //upFrequency_ = size_; 00261 frameMaxNumPeaks_ = upFrequency_-downFrequency_; 00262 } 00263 else 00264 frameMaxNumPeaks_ = ctrl_frameMaxNumPeaks_->to<mrs_natural>(); 00265 00266 setctrl(ctrl_onSamples_, ctrl_inSamples_); 00267 setctrl(ctrl_onObservations_, frameMaxNumPeaks_*nbParameters_); 00268 setctrl(ctrl_osrate_, ctrl_israte_); 00269 00270 ostringstream oss; 00271 for(mrs_natural j=0; j< nbParameters_; ++j) //j = param index 00272 { 00273 for (mrs_natural i=0; i < frameMaxNumPeaks_; i++) //i = peak index 00274 oss << peakView::getParamName(j) << "_" << i+j*frameMaxNumPeaks_ << ","; 00275 } 00276 ctrl_onObsNames_->setValue(oss.str(), NOUPDATE); 00277 00278 if (getctrl("mrs_real/peakSmearingTimeInS")->to<mrs_real>() == 0 || !pick_) 00279 lpCoeff_ = 0; 00280 else 00281 lpCoeff_ = exp(-2.2/(timeSrate/hopSize_*getctrl("mrs_real/peakSmearingTimeInS")->to<mrs_real>())); 00282 00283 if (size_ != psize_) 00284 { 00285 tmpBuff_.stretch(inObservations_); 00286 lastphase_.stretch(size_); 00287 phase_.stretch(size_); 00288 mag_.stretch(size_); 00289 masked_.stretch(size_,1); 00290 lpPeakerRes_.stretch(size_,1); 00291 magCorr_.stretch(size_); 00292 frequency_.stretch(size_); 00293 lastmag_.stretch(size_); 00294 lastfrequency_.stretch(size_); 00295 deltamag_.stretch(size_); 00296 deltafrequency_.stretch(size_); 00297 psize_ = size_; 00298 00299 lpPeakerRes_.setval (1.); 00300 } 00301 00302 factor_ = timeSrate / TWOPI / instFreqHopSize_; 00303 fundamental_ = israte_; 00304 00305 peakProb_.stretch (3,1); 00306 peakProbWeight_ = getctrl("mrs_realvec/peakProbabilityWeight")->to<mrs_realvec>(); 00307 if (peakProbWeight_.getRows () > peakProbWeight_.getCols ()) 00308 peakProbWeight_.transpose (); 00309 peakProbWeight_ /= peakProbWeight_.sum (); 00310 00311 // init lastfreq to avoid distance inconsistencies 00312 for (mrs_natural k = 0; k < size_; k++) 00313 lastfrequency_(k) = k*fundamental_; 00314 } 00315 00316 mrs_real 00317 PeakConvert2::lobe_value_compute(mrs_real f, mrs_natural type, mrs_natural size) 00318 { 00319 mrs_real re ; 00320 00321 // size par size-2 !!! 00322 switch (type) 00323 { 00324 case 1: 00325 { 00326 re= fabs (0.5*lobe_value_compute(f, 0, size)+ 00327 0.25*lobe_value_compute(f-2.*PI/size, 0, size)+ 00328 0.25*lobe_value_compute(f+2.*PI/size, 0, size))/size ; 00329 return fabs(re); 00330 } 00331 case 0: 00332 return (mrs_real) (f == 0) ? size : (sin(f*0.5*(size))/sin(f*0.5)); 00333 default: 00334 { 00335 return 0.0 ; 00336 } 00337 } 00338 } 00339 00340 void 00341 PeakConvert2::getShortBinInterval(realvec& interval, realvec& index, realvec& mag) 00342 { 00343 const unsigned int maxLobeWidth = 6; 00344 unsigned int nbP=index.getSize(); 00345 unsigned int minIndex = 0, 00346 endLoop, 00347 length = mag.getSize (); 00348 00349 // we could also use the instantaneous frequency here...? 00350 00351 for(unsigned int i=0 ; i<nbP ; i++) 00352 { 00353 unsigned int idx = (unsigned int)(index(i)+.1); 00354 00355 if (idx <= 0) 00356 continue; 00357 00358 endLoop = min(length,idx + maxLobeWidth); 00359 minIndex = endLoop; 00360 // look for the next valley location upward 00361 for (unsigned int j = idx ; j < endLoop ; j++) 00362 { 00363 if(mag(j) < mag(j+1)) 00364 { 00365 minIndex = j; 00366 break; 00367 } 00368 } 00369 00370 interval(2*i+1) = minIndex; 00371 00372 endLoop = max((unsigned int)0,idx - maxLobeWidth); 00373 minIndex = endLoop; 00374 00375 // look for the next valley location downward 00376 for (unsigned int j= idx ; j > endLoop ; j--) 00377 { 00378 if(mag(j) < mag(j-1)) 00379 { 00380 minIndex = j; 00381 break; 00382 } 00383 } 00384 00385 interval(2*i) = minIndex; 00386 } 00387 } 00388 00389 00390 void PeakConvert2::ComputeMasking (realvec& in) 00391 { 00392 (void) in; // [!] what was this supposed to do? 00393 masking_->updControl ("mrs_natural/inObservations", size_); 00394 masking_->updControl ("mrs_natural/inSamples", 1); 00395 masking_->updControl ("mrs_real/israte", israte_); 00396 00397 mag_.transpose(); 00398 masking_->myProcess (mag_,masked_); 00399 mag_.transpose(); 00400 } 00401 00402 void PeakConvert2::ComputeMagnitudeAndPhase (mrs_realvec in) 00403 { 00404 mrs_real a, c; 00405 mrs_real b, d; 00406 mrs_real phasediff; 00407 for (mrs_natural o=0; o < size_; o++) 00408 { 00409 if (o==0) //DC bins 00410 { 00411 a = in(0); 00412 b = 0.0; 00413 c = in(N_); 00414 d = 0.0; 00415 } 00416 else if (o == size_-1) //Nyquist bins 00417 { 00418 a = in(1); 00419 b = 0.0; 00420 c = in(N_+1); 00421 d = 0.0; 00422 } 00423 else //all other bins 00424 { 00425 a = in(2*o); 00426 b = in(2*o+1); 00427 c = in(N_+2*o); 00428 d = in(N_+2*o+1); 00429 } 00430 00431 if (o < downFrequency_ || o > upFrequency_) 00432 { 00433 frequency_(o) = 0; 00434 mag_(o) = sqrt((a*a + b*b))*2; 00435 continue; 00436 } 00437 if ( a == .0 || c == .0) 00438 { 00439 frequency_(o) = o*fundamental_; 00440 } 00441 else 00442 { 00443 if(prec_ && pick_) 00444 { 00445 mrs_real Omega = TWOPI*o*instFreqHopSize_/N_; // now works with hopsizes != 1 as well (AL) 00446 00447 // compute phase 00448 phase_(o) = atan2(b,a); 00449 00450 // compute precise frequency using the phase difference 00451 lastphase_(o) = atan2(d,c); 00452 phasediff = princArg(phase_(o)-lastphase_(o) - Omega) + Omega; 00453 frequency_(o) = abs(phasediff * factor_ ); 00454 } 00455 else 00456 frequency_(o) = o*fundamental_; 00457 } 00458 00459 00460 // compute precise amplitude 00461 mag_(o) = sqrt((a*a + b*b))*2; 00462 if (pick_) 00463 { 00464 mrs_real mag = lobe_value_compute((o * fundamental_-frequency_(o))/factor_, 1, N_); 00465 magCorr_(o) = mag_(o)/mag; 00466 } 00467 else 00468 { 00469 magCorr_(o) = mag_(o); 00470 } 00471 //mrs_real freq = frequency_(o); 00472 00473 if(frequency_(o) != 0.0) 00474 { 00475 const mrs_natural searchRange = 8; 00476 mrs_natural k, 00477 kd = o, 00478 ku = o, 00479 kEndSearch; 00480 mrs_real diff[2]; 00481 00482 // find closest preceding frequency 00483 kEndSearch = max ((mrs_natural)0, o-searchRange); 00484 for (k = o-1; k > kEndSearch; k--) 00485 { 00486 diff[0] = abs(frequency_(o) - lastfrequency_(k)); 00487 diff[1] = abs(frequency_(o) - lastfrequency_(kd)); 00488 if (diff[0] < diff[1]) 00489 kd = k; 00490 } 00491 kEndSearch = min (size_, o+searchRange); 00492 for (k = o+1; k < kEndSearch; k++) 00493 { 00494 diff[0] = abs(frequency_(o) - lastfrequency_(k)); 00495 diff[1] = abs(frequency_(o) - lastfrequency_(ku)); 00496 if (diff[0] < diff[1]) 00497 ku = k; 00498 } 00499 diff[0] = abs(frequency_(o) - lastfrequency_(kd)); 00500 diff[1] = abs(frequency_(o) - lastfrequency_(ku)); 00501 k = (diff[0] < diff[1])? kd : ku; 00502 00503 deltafrequency_(o) = (lastfrequency_(k) == 0)? .0 : log10(lastfrequency_(k)/frequency_(o)); 00504 //deltafrequency_(o) = (frequency_(o)-lastfrequency_(o))/(frequency_(o)+lastfrequency_(o)); 00505 } 00506 else 00507 deltafrequency_(o) = .0; 00508 00509 mrs_real lastmag = lastmag_(o); 00510 if (o > 0) 00511 lastmag = max(lastmag_(o-1), lastmag); 00512 if (o < size_-1) 00513 lastmag = max(lastmag_(o+1), lastmag); 00514 00515 if (mag_(o) > 0) 00516 deltamag_(o) = (mag_(o)-lastmag)/mag_(o); 00517 else if (lastmag > 0) 00518 deltamag_(o) = (mag_(o)-lastmag)/lastmag; 00519 else 00520 deltamag_(o) = 0; 00521 } 00522 lastfrequency_ = frequency_; 00523 lastmag_ = mag_; 00524 } 00525 00526 void PeakConvert2::ComputePeaker (mrs_realvec in, realvec& out) 00527 { 00528 #ifdef ORIGINAL_VERSION 00529 peaker_->updControl("mrs_real/peakStrength", 0.2);// to be set as a control [!] 00530 #else 00531 peaker_->updControl("mrs_real/peakStrength",1e-1); 00532 peaker_->updControl("mrs_real/peakStrengthRelMax" ,1e-2); 00533 peaker_->updControl("mrs_real/peakStrengthAbs",1e-10 ); 00534 peaker_->updControl("mrs_real/peakStrengthTreshLpParam" ,0.95); 00535 peaker_->updControl("mrs_real/peakStrengthRelThresh" , 1.); 00536 #endif 00537 00538 peaker_->updControl("mrs_real/peakSpacing", 2e-3); // 0 00539 peaker_->updControl("mrs_natural/peakStart", downFrequency_); // 0 00540 peaker_->updControl("mrs_natural/peakEnd", upFrequency_); // size_ 00541 peaker_->updControl("mrs_natural/inSamples", in.getCols()); 00542 peaker_->updControl("mrs_natural/inObservations", in.getRows()); 00543 peaker_->updControl("mrs_natural/onSamples", out.getCols()); 00544 peaker_->updControl("mrs_natural/onObservations", out.getRows()); 00545 00546 peaker_->process(in, out); 00547 } 00548 00549 void 00550 PeakConvert2::myProcess(realvec& in, realvec& out) 00551 { 00552 mrs_natural o,i; 00553 out.setval(0); 00554 peakView pkViewOut(out); 00555 00556 const mrs_bool useMasking = getctrl("mrs_bool/useMasking")->to<mrs_bool>(); 00557 const mrs_real probThresh = getctrl("mrs_real/probabilityTresh")->to<mrs_real>(); 00558 00559 max_->updControl("mrs_natural/nMaximums", frameMaxNumPeaks_); 00560 00561 max_->setctrl("mrs_natural/inSamples", size_); 00562 max_->setctrl("mrs_natural/inObservations", 1); 00563 max_->update(); 00564 tmp_.stretch(frameMaxNumPeaks_*2); 00565 00566 for(mrs_natural f=0 ; f < inSamples_; ++f) 00567 { 00568 //we should avoid the first empty frames, 00569 //that will contain silence and consequently create 00570 //discontinuities in the signal, ruining the peak calculation! 00571 //only process if we have a full data vector (i.e. no zeros) 00572 if(frame_ >= skip_) 00573 { 00574 // get pair of ffts 00575 in.getCol (f, tmpBuff_); 00576 00577 // compute magnitude, phase, and instantaneous frequency 00578 this->ComputeMagnitudeAndPhase (tmpBuff_); 00579 00580 // compute masking threshold 00581 if (useMasking && pick_) 00582 ComputeMasking (tmpBuff_); 00583 else 00584 masked_.setval(10.); 00585 00586 // select bins with local maxima in magnitude (--> peaks) 00587 peaks_ = mag_; 00588 if(pick_) 00589 this->ComputePeaker (mag_, peaks_); 00590 else 00591 { 00592 for (o = 0 ; o < downFrequency_ ; o++) 00593 peaks_(o)=0.0; 00594 for (o = upFrequency_ ; o < (mrs_natural)peaks_.getSize() ; o++) 00595 peaks_(o)=0.0; 00596 } 00597 00598 if (lpCoeff_ > 0) 00599 FreqSmear (lpPeakerRes_); 00600 00601 //compute the probability of a peak being a peak 00602 for(o=0 ; o < size_ ; o++) 00603 { 00604 if (peaks_(o) <= 0) 00605 { 00606 frequency_(o) = .0; 00607 //lastmag_(o) = .0; 00608 lastfrequency_(o) = .0; 00609 // time smearing if no new peak 00610 lpPeakerRes_(o) *=lpCoeff_; 00611 continue; 00612 } 00613 #ifdef ORIGINAL_VERSION 00614 // probability of peak being a masker 00615 peakProb_(0) = 0; 00616 // probability of peak being stationary 00617 peakProb_(1) = 0; 00618 // probability of peak being tonal 00619 peakProb_(2) = (abs(frequency_(o)/fundamental_-o) > .5)? 0 : 1; 00620 #else 00621 // probability of peak being a masker 00622 peakProb_(0) = max((mrs_real).1, (mrs_real).5 * (mrs_real)(log10(masked_(o)) +1.)); 00623 // probability of peak being stationary 00624 peakProb_(1) = max((mrs_real).1, (mrs_real)lpPeakerRes_(o)); 00625 // probability or peak being tonal 00626 peakProb_(2) = GaussianPdf (frequency_(o)/fundamental_-o, gaussianStd); 00627 #endif 00628 00629 // reset lpPeakerRes with peaker results 00630 lpPeakerRes_(o) = 1; 00631 00632 peakProb_ *= peakProbWeight_; 00633 if ((peakProb_.sum() < probThresh) && pick_) 00634 { 00635 peaks_(o) = .0; 00636 frequency_(o) = .0; 00637 //lastmag_(o) = .0; 00638 lastfrequency_(o) = .0; 00639 } 00640 } 00641 00642 // keep only the frameMaxNumPeaks_ highest amplitude local maxima 00643 tmp_.setval(0.); 00644 max_->process(peaks_, tmp_); 00645 00646 nbPeaks_=tmp_.getSize()/2; 00647 realvec index_(nbPeaks_); //[!] make member to avoid reallocation at each tick! 00648 for (i=0 ; i<nbPeaks_ ; i++) 00649 index_(i) = tmp_(2*i+1); 00650 00651 // search for bins interval 00652 realvec interval_(nbPeaks_*2); //[!] make member to avoid reallocation at each tick! 00653 interval_.setval(0); 00654 if(pick_) 00655 getShortBinInterval(interval_, index_, mag_); 00656 else 00657 { 00658 for (i=0 ; i<nbPeaks_ ; i++) 00659 interval_(2*i+1) = index_(i); 00660 } 00661 00662 #ifdef LOG2FILE 00663 for (i=0 ; i<nbPeaks_ ; i++) 00664 { 00665 mrs_real value = frequency_((mrs_natural) (index_(i)+.1)); 00666 pFDbgFile << std::scientific << std::setprecision(4) << value << "\t"; 00667 } 00668 pFDbgFile << std::endl; 00669 #endif 00670 #ifdef MARSYAS_MATLAB 00671 #ifdef MTLB_DBG_LOG 00672 MATLAB_PUT(mag_, "peaks"); 00673 MATLAB_PUT(peaks_, "k"); 00674 MATLAB_PUT(tmp_, "tmp"); 00675 MATLAB_PUT(interval_, "int"); 00676 MATLAB_PUT(frequency_, "freq"); 00677 // MATLAB_EVAL("figure(1);clf;hold on ;plot(peaks);stem(k);stem(tmp(2:2:end)+1, peaks(tmp(2:2:end)+1), 'r')"); 00678 // MATLAB_EVAL("stem(int+1, peaks(int+1), 'k')"); 00679 MATLAB_EVAL("figure(1);hold on ;stem(freq(tmp(2:2:end)+1), peaks(tmp(2:2:end)+1), 'r');hold off"); 00680 #endif 00681 #endif 00682 00683 00684 // fill output with peaks data 00685 interval_ /= N_; 00686 00687 for (i = 0; i < nbPeaks_; i++) 00688 { 00689 mrs_natural index = (mrs_natural) (index_(i)+.1); 00690 pkViewOut(i, peakView::pkFrequency, f) = frequency_(index); 00691 pkViewOut(i, peakView::pkAmplitude, f) = magCorr_(index); 00692 pkViewOut(i, peakView::pkPhase, f) = -phase_(index); 00693 pkViewOut(i, peakView::pkDeltaFrequency, f) = deltafrequency_(index); 00694 pkViewOut(i, peakView::pkDeltaAmplitude, f) = /*abs*/(deltamag_(index)); 00695 pkViewOut(i, peakView::pkFrame, f) = frame_; 00696 pkViewOut(i, peakView::pkGroup, f) = 0.;//(pick_)?-1.:0.; //This should be -1!!!! [TODO] 00697 pkViewOut(i, peakView::pkVolume, f) = 1.0; 00698 pkViewOut(i, peakView::pkBinLow, f) = interval_(2*i); 00699 pkViewOut(i, peakView::pkBin, f) = index_(i); 00700 pkViewOut(i, peakView::pkBinHigh, f) = interval_(2*i+1); 00701 pkViewOut(i, peakView::pkTrack, f) = -1.0; //null-track ID 00702 00703 MRSASSERT((index_(i) <= interval_(2*i)) || (interval_(2*i+1) <= index_(i))); 00704 00705 if(useStereoSpectrum_) 00706 pkViewOut(i, peakView::pkPan, f) = in((mrs_natural)index_(i)+2*N_, f); 00707 else 00708 pkViewOut(i, peakView::pkPan, f) = 0.0; 00709 } 00710 } 00711 else //if not yet reached "skip" number of frames... 00712 { 00713 for(mrs_natural i=0; i< frameMaxNumPeaks_; ++i) 00714 { 00715 //pkViewOut(i, peakView::pkFrequency, f) = 0; 00716 //pkViewOut(i, peakView::pkAmplitude, f) = 0; 00717 //pkViewOut(i, peakView::pkPhase, f) = 0; 00718 //pkViewOut(i, peakView::pkDeltaFrequency, f) = 0; 00719 //pkViewOut(i, peakView::pkDeltaAmplitude, f) = 0; 00720 pkViewOut(i, peakView::pkFrame, f) = frame_; 00721 //pkViewOut(i, peakView::pkGroup, f) = -1; 00722 //pkViewOut(i, peakView::pkVolume, f) = 0; 00723 //pkViewOut(i, peakView::pkPan, f) = 0; 00724 //pkViewOut(i, peakView::pkBinLow, f) = 0; 00725 //pkViewOut(i, peakView::pkBin, f) = 0; 00726 //pkViewOut(i, peakView::pkBinHigh, f) = 0; 00727 } 00728 } 00729 frame_++; 00730 } 00731 00732 //count the total number of existing peaks (i.e. peak freq != 0) 00733 ctrl_totalNumPeaks_->setValue(pkViewOut.getTotalNumPeaks()); 00734 }