Marsyas  0.5.0-beta1
/Users/jleben/code/marsyas/src/marsyas/marsystems/PeakConvert2.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 #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 }