Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/PeakConvert.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2005 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 "PeakConvert.h"
00021 #include "Peaker.h"
00022 #include "MaxArgMax.h"
00023 #include <marsyas/peakView.h>
00024 
00025 #include <algorithm>
00026 
00027 using std::ostringstream;
00028 using std::abs;
00029 
00030 using namespace Marsyas;
00031 
00032 PeakConvert::PeakConvert(mrs_string name):MarSystem("PeakConvert",name)
00033 {
00034   psize_ = 0;
00035   size_ = 0;
00036   nbParameters_ = peakView::nbPkParameters;
00037   skip_=0;
00038   frame_ = 0;
00039   N_ = 0;
00040 
00041   fundamental_ = 0.0;
00042   factor_ = 0.0;
00043   nbPeaks_ = 0;
00044   frameMaxNumPeaks_ = 0;
00045 
00046   useStereoSpectrum_ = false;
00047 
00048   peaker_ = new Peaker("Peaker");
00049   max_ = new MaxArgMax("MaxArgMax");
00050 
00051   addControls();
00052 }
00053 
00054 PeakConvert::PeakConvert(const PeakConvert& a) : MarSystem(a)
00055 {
00056   psize_ = a.psize_;
00057   size_ = a.size_;
00058   nbParameters_ = a.nbParameters_;
00059   skip_ = a.skip_;
00060   frame_ = a.frame_;
00061   N_ = a.N_;
00062 
00063   fundamental_ = a.fundamental_;
00064   factor_ = a.factor_;
00065   nbPeaks_ = a.nbPeaks_;
00066   frameMaxNumPeaks_  = a.frameMaxNumPeaks_;
00067 
00068   useStereoSpectrum_ = a.useStereoSpectrum_;
00069 
00070   peaker_ = new Peaker("Peaker");
00071   max_ = new MaxArgMax("MaxArgMax");
00072 
00073   ctrl_totalNumPeaks_ = getctrl("mrs_natural/totalNumPeaks");
00074   ctrl_frameMaxNumPeaks_ = getctrl("mrs_natural/frameMaxNumPeaks");
00075 }
00076 
00077 PeakConvert::~PeakConvert()
00078 {
00079   delete peaker_;
00080   delete max_;
00081 }
00082 
00083 MarSystem*
00084 PeakConvert::clone() const
00085 {
00086   return new PeakConvert(*this);
00087 }
00088 
00089 void
00090 PeakConvert::addControls()
00091 {
00092   addctrl("mrs_natural/frameMaxNumPeaks", 0);
00093   setctrlState("mrs_natural/frameMaxNumPeaks", true);
00094 
00095   addctrl("mrs_string/frequencyInterval", "MARSYAS_EMPTY");
00096   setctrlState("mrs_string/frequencyInterval", true);
00097 
00098   addctrl("mrs_natural/nbFramesSkipped", 0);
00099   setctrlState("mrs_natural/nbFramesSkipped", true);
00100 
00101   addctrl("mrs_bool/improvedPrecision", true);
00102   setctrlState("mrs_bool/improvedPrecision", true);
00103 
00104   addctrl("mrs_bool/picking", true);
00105   setctrlState("mrs_bool/picking", true);
00106 
00107   addctrl("mrs_natural/totalNumPeaks", 0, ctrl_totalNumPeaks_);
00108 }
00109 
00110 void
00111 PeakConvert::myUpdate(MarControlPtr sender)
00112 {
00113   (void) sender;  //suppress warning of unused parameter(s)
00114 
00115   //check the input to see if we are also getting stereo information
00116   //(N_ is the FFT size)
00117   if (fmod(inObservations_, 2.0) == 0.0)
00118   {
00119     //we just have the two shifted spectrums stacked vertically
00120     //(input has N + N observations)
00121     N_ = inObservations_/2;
00122     useStereoSpectrum_ = false;
00123   }
00124   else if(fmod(inObservations_-1, 2.5) == 0.0)
00125   {
00126     //we also have stereo spectrum info at the bottom
00127     //(input has N + N + N/2+1 observations)
00128     N_ = (mrs_natural)((inObservations_-1) / 2.5);
00129     useStereoSpectrum_ = true;
00130   }
00131 
00132   //if picking is disabled (==false), the number of sinusoids should be set
00133   //to the number of unique bins of the spectrums at the input (i.e. N/2+1)
00134   pick_ = getctrl("mrs_bool/picking")->to<mrs_bool>();
00135   if(!pick_ && ctrl_frameMaxNumPeaks_->to<mrs_natural>() == 0)
00136     frameMaxNumPeaks_ = N_/2+1; //inObservations_/4+1;
00137   else
00138     frameMaxNumPeaks_ = ctrl_frameMaxNumPeaks_->to<mrs_natural>();
00139 
00140   setctrl(ctrl_onSamples_, ctrl_inSamples_);
00141   setctrl(ctrl_onObservations_, frameMaxNumPeaks_*nbParameters_);
00142   setctrl(ctrl_osrate_, ctrl_israte_);
00143 
00144   ostringstream oss;
00145   for(mrs_natural j=0; j< nbParameters_; ++j) //j = param index
00146   {
00147     for (mrs_natural i=0; i < frameMaxNumPeaks_; ++i) //i = peak index
00148       oss << peakView::getParamName(j) << "_" << i+j*frameMaxNumPeaks_ << ",";
00149   }
00150   ctrl_onObsNames_->setValue(oss.str(), NOUPDATE);
00151 
00152   mrs_real timeSrate = israte_*(mrs_real)N_;//israte_*(mrs_real)inObservations_/2.0;
00153 
00154   size_ = N_/2+1;//inObservations_ /4 +1;
00155   if (size_ != psize_)
00156   {
00157     lastphase_.stretch(size_);
00158     phase_.stretch(size_);
00159     mag_.stretch(size_);
00160     magCorr_.stretch(size_);
00161     frequency_.stretch(size_);
00162     lastmag_.stretch(size_);
00163     lastfrequency_.stretch(size_);
00164     deltamag_.stretch(size_);
00165     deltafrequency_.stretch(size_);
00166     psize_ = size_;
00167   }
00168 
00169   factor_ = timeSrate / TWOPI;
00170   fundamental_ = israte_;
00171 
00172   skip_ = getctrl("mrs_natural/nbFramesSkipped")->to<mrs_natural>();
00173   prec_ = getctrl("mrs_bool/improvedPrecision")->to<mrs_bool>();
00174 
00175   if(getctrl("mrs_string/frequencyInterval")->to<mrs_string>() != "MARSYAS_EMPTY")
00176   {
00177     realvec conv(2);
00178     string2parameters(getctrl("mrs_string/frequencyInterval")->to<mrs_string>(), conv, '_'); //[!]
00179     downFrequency_ = (mrs_natural) floor(conv(0)/timeSrate*size_*2) ;
00180     upFrequency_ = (mrs_natural) floor(conv(1)/timeSrate*size_*2);
00181   }
00182   else
00183   {
00184     downFrequency_ = 0;
00185     upFrequency_ = size_;
00186   }
00187 }
00188 
00189 mrs_real
00190 PeakConvert::lobe_value_compute(mrs_real f, mrs_natural type, mrs_natural size)
00191 {
00192   mrs_real re ;
00193 
00194   // size par size-2 !!!
00195   switch (type)
00196   {
00197   case 1:
00198   {
00199     re= fabs (0.5*lobe_value_compute(f, 0, size)+
00200               0.25*lobe_value_compute(f-2.*PI/size, 0, size)+
00201               0.25*lobe_value_compute(f+2.*PI/size, 0, size))/size ;
00202     return fabs(re);
00203   }
00204   case 0:
00205     return (mrs_real) (f == 0) ? size : (sin(f*0.5*(size))/sin(f*0.5));
00206   default:
00207   {
00208     return 0.0 ;
00209   }
00210   }
00211 }
00212 
00213 /*
00214 void
00215 PeakConvert::getBinInterval(realvec& interval, realvec& index, realvec& mag) //[WTF] is this being used for anything?!?
00216 {
00217     mrs_natural k=0, start=0, nbP=index.getSize();
00218     mrs_natural minIndex = 0;
00219 
00220     // getting rid of padding zeros
00221     while(start<index.getSize() && !index(start))
00222         start++;
00223 
00224     for(mrs_natural i=start ; i<nbP ; ++i, k++)
00225     {
00226         interval(2*k) = index(i)-1;
00227         interval(2*k+1) = index(i);
00228     }
00229 }*/
00230 
00231 void
00232 PeakConvert::getShortBinInterval(realvec& interval, realvec& index, realvec& mag)
00233 {
00234   mrs_natural k=0, start=0, nbP=index.getSize();
00235   mrs_natural minIndex = 0;
00236 
00237   // getting rid of padding zeros
00238   while(start<index.getSize() && !index(start))
00239     start++;
00240 
00241   for(mrs_natural i=start ; i<nbP ; ++i, ++k)
00242   {
00243     minIndex = 0;
00244     // look for the next valley location upward
00245     for (mrs_natural j= (mrs_natural)index(i) ; j<mag.getSize()-1 ; ++j)
00246     {
00247       if(mag(j) < mag(j+1))
00248       {
00249         minIndex = j;
00250         break;
00251       }
00252     }
00253 //      if(!minIndex) //arght!!! I hate using logic with integers!!! Makes code so difficult to read!!! [!]
00254 //      {
00255 //          cout << "pb while looking for bin intervals" << endl; //[?]
00256 //      }
00257 
00258     interval(2*k+1) = minIndex;
00259 
00260     // look for the next valley location downward
00261     for (mrs_natural j= (mrs_natural)index(i) ; j>1 ; --j)
00262     {
00263       if(mag(j) < mag(j-1))
00264       {
00265         minIndex = j;
00266         break;
00267       }
00268     }
00269 //      if(!minIndex) //arght!!! I hate using logic with integers!!! Makes code so difficult to read!!! [!]
00270 //      {
00271 //          cout << "pb while looking for bin intervals" << endl; //[?]
00272 //      }
00273     interval(2*k) = minIndex;
00274   }
00275 }
00276 
00277 
00278 void
00279 PeakConvert::getLargeBinInterval(realvec& interval, realvec& index, realvec& mag)
00280 {
00281   mrs_natural k=0, start=0, nbP=index.getSize();
00282 
00283   // handling the first case
00284   mrs_real minVal = HUGE_VAL;
00285   mrs_natural minIndex = 0;
00286 
00287   // getting rid of padding zeros
00288   while(!index(start))
00289     start++;
00290 
00291   for (mrs_natural j=0 ; j<index(start) ; j++) //is this foor loop like this?!?!?!?!?!?!?!?! [!]
00292   {
00293     if(minVal > mag(j))
00294     {
00295       minVal = mag(j);
00296       minIndex = j;
00297     }
00298   }
00299 //  if(!minIndex)
00300 //  {
00301 //      cout << "pb while looking for minimal bin intervals" << endl; //[WTF]
00302 //  }
00303   interval(0) = minIndex;
00304 
00305   for(mrs_natural i=start ; i<nbP-1 ; ++i, k++)
00306   {
00307     minVal = HUGE_VAL;
00308     minIndex = 0;
00309     // look for the minimal value among successive peaks
00310     for (mrs_natural j= (mrs_natural) index(i) ; j<index(i+1) ; j++) // is this for loop like this?!?!?! [?]
00311     {
00312       if(minVal > mag(j))
00313       {
00314         minVal = mag(j);
00315         minIndex = j;
00316       }
00317     }
00318 
00319 //      if(!minIndex)
00320 //      {
00321 //          cout << "pb while looking for bin intervals" << endl; //[WTF]
00322 //      }
00323     interval(2*k+1) = minIndex-1;
00324     interval(2*(k+1)) = minIndex;
00325   }
00326 
00327   // handling the last case
00328   minVal = HUGE_VAL;
00329   minIndex = 0;
00330   for (mrs_natural j= (mrs_natural)index(nbP-1) ; j<mag.getSize()-1 ; ++j)
00331   {
00332     if(minVal > mag(j))
00333     {
00334       minVal = mag(j);
00335       minIndex = j;
00336     }
00337     // consider stopping the search at the first valley
00338     if(minVal<mag(j+1))
00339       break;
00340   }
00341 //  if(!minIndex)
00342 //  {
00343 //      cout << "pb while looking for maximal bin intervals" << endl; //[WTF]
00344 //  }
00345   interval(2*(k)+1) = minIndex;
00346 }
00347 
00348 void
00349 PeakConvert::myProcess(realvec& in, realvec& out)
00350 {
00351   mrs_natural o;
00352   mrs_real a, c;
00353   mrs_real b, d;
00354   mrs_real phasediff;
00355 
00356   out.setval(0);
00357   peakView pkViewOut(out);
00358 
00359   for(mrs_natural f=0 ; f < inSamples_; ++f)
00360   {
00361     //we should avoid the first empty frames,
00362     //that will contain silence and consequently create
00363     //discontinuities in the signal, ruining the peak calculation!
00364     //only process if we have a full data vector (i.e. no zeros)
00365     if(frame_ >= skip_)
00366     {
00367       // handle amplitudes from shifted spectrums at input
00368       for (o=0; o < size_; o++)
00369       {
00370         if (o==0) //DC bins
00371         {
00372           a = in(0,f);
00373           b = 0.0;
00374           c = in(N_, f);
00375           d = 0.0;
00376         }
00377         else if (o == size_-1) //Nyquist bins
00378         {
00379           a = in(1, f);
00380           b = 0.0;
00381           c = in(N_+1, f);
00382           d = 0.0;
00383         }
00384         else //all other bins
00385         {
00386           a = in(2*o, f);
00387           b = in(2*o+1, f);
00388           c = in(N_+2*o, f);
00389           d = in(N_+2*o+1, f);
00390         }
00391 
00392         // computer magnitude value
00393         //mrs_real par = lobe_value_compute (0, 1, 2048); //[?]
00394 
00395         // compute phase
00396         phase_(o) = atan2(b,a);
00397 
00398         // compute precise frequency using the phase difference
00399         lastphase_(o)= atan2(d,c);
00400         if(phase_(o) >= lastphase_(o))
00401           phasediff = phase_(o)-lastphase_(o);
00402         else
00403           phasediff = phase_(o)-lastphase_(o)+TWOPI;
00404         if(prec_)
00405           frequency_(o) = phasediff * factor_ ;
00406         else
00407           frequency_(o) = o*fundamental_;
00408 
00409         // compute precise amplitude
00410         mag_(o) = sqrt((a*a + b*b))*2; //*4/0.884624;//*50/3); // [!]
00411         mrs_real mag = lobe_value_compute((o * fundamental_-frequency_(o))/factor_, 1, N_);
00412         magCorr_(o) = mag_(o)/mag;
00413 
00414         // computing precise frequency using the derivative method // use at your own risk  [?]
00415         /*  mrs_real lastmag = sqrt(c*c + d*d);
00416         mrs_real rap = (mag_(o)-lastmag)/(lastmag*2);
00417         f=asin(rap);
00418         f *= (getctrl("mrs_real/israte")->to<mrs_real>()*inObservations/2.0)/PI;
00419         */
00420         // rough frequency and amplitude
00421         //frequency_(o) = o * fundamental_;
00422         //magCorr_(o) = mag_(o);
00423 
00424         if(lastfrequency_(o) != 0.0)
00425           deltafrequency_(o) = (frequency_(o)-lastfrequency_(o))/(frequency_(o)+lastfrequency_(o));
00426 
00427         deltamag_(o) = (mag_(o)-lastmag_(o))/(mag_(o)+lastmag_(o));
00428 
00429         // remove potential peak if frequency too irrelevant
00430         if(abs(frequency_(o)-o*fundamental_) > 0.5*fundamental_)
00431           frequency_(o)=0.0;
00432 
00433         lastfrequency_(o) = frequency_(o);
00434         lastmag_(o) = mag_(o);
00435       }
00436 
00437       // select bins with local maxima in magnitude (--> peaks)
00438       realvec peaks_ = mag_;
00439       realvec tmp_;
00440       peaker_->updControl("mrs_real/peakStrength", 0.2);// to be set as a control [!]
00441       peaker_->updControl("mrs_natural/peakStart", downFrequency_);   // 0
00442       peaker_->updControl("mrs_natural/peakEnd", upFrequency_);  // size_
00443       peaker_->updControl("mrs_natural/inSamples", mag_.getCols());
00444       peaker_->updControl("mrs_natural/inObservations", mag_.getRows());
00445       peaker_->updControl("mrs_natural/onSamples", peaks_.getCols());
00446       peaker_->updControl("mrs_natural/onObservations", peaks_.getRows());
00447       if(pick_)
00448         peaker_->process(mag_, peaks_);
00449       else
00450       {
00451         //peaks_ = mag_;
00452         for (o = 0 ; o < downFrequency_ ; o++)
00453           peaks_(o)=0.0;
00454         for (o = upFrequency_ ; o < (mrs_natural)peaks_.getSize() ; ++o)
00455           peaks_(o)=0.0;
00456       }
00457 
00458       //discard bins whose frequency was set as irrelevant...
00459       for(o=0 ; o < size_ ; o++)
00460       {
00461         if(frequency_(o) == 0)
00462           peaks_(o) = 0.0;
00463       }
00464 
00465       // keep only the frameMaxNumPeaks_ highest amplitude local maxima
00466       if(ctrl_frameMaxNumPeaks_->to<mrs_natural>() != 0) //?????????????????? is this check needed?!? See myUpdate
00467       {
00468         tmp_.stretch(frameMaxNumPeaks_*2);
00469         max_->updControl("mrs_natural/nMaximums", frameMaxNumPeaks_);
00470       }
00471       else //?????????????????? is this check needed?!? See myUpdate
00472       {
00473         tmp_.stretch((upFrequency_-downFrequency_)*2);
00474         max_->updControl("mrs_natural/nMaximums", upFrequency_-downFrequency_);
00475       }
00476       max_->setctrl("mrs_natural/inSamples", size_);
00477       max_->setctrl("mrs_natural/inObservations", 1);
00478       max_->update();
00479       max_->process(peaks_, tmp_);
00480 
00481       nbPeaks_=tmp_.getSize()/2;
00482       realvec index_(nbPeaks_); //[!] make member to avoid reallocation at each tick!
00483       for (mrs_natural i=0 ; i<nbPeaks_ ; ++i)
00484         index_(i) = tmp_(2*i+1);
00485       realvec index2_ = index_;
00486       index2_.sort();
00487 
00488       // search for bins interval
00489       realvec interval_(nbPeaks_*2); //[!] make member to avoid reallocation at each tick!
00490       interval_.setval(0);
00491       if(pick_)
00492         getShortBinInterval(interval_, index2_, mag_);
00493 
00494       // fill output with peaks data
00495       /*
00496       MATLAB_PUT(mag_, "peaks");
00497       MATLAB_PUT(peaks_, "k");
00498       MATLAB_PUT(tmp_, "tmp");
00499       MATLAB_PUT(interval_, "int");
00500       MATLAB_EVAL("figure(1);clf;hold on ;plot(peaks);stem(k);stem(tmp(2:2:end)+1, peaks(tmp(2:2:end)+1), 'r')");
00501       MATLAB_EVAL("stem(int+1, peaks(int+1), 'k')");
00502       */
00503 
00504       interval_ /= N_*2;
00505 
00506       for (mrs_natural i = 0; i < nbPeaks_; ++i)
00507       {
00508         pkViewOut(i, peakView::pkFrequency, f) = frequency_((mrs_natural) index_(i));
00509         pkViewOut(i, peakView::pkAmplitude, f) = magCorr_((mrs_natural) index_(i));
00510         pkViewOut(i, peakView::pkPhase, f) = -phase_((mrs_natural) index_(i));
00511         pkViewOut(i, peakView::pkDeltaFrequency, f) = deltafrequency_((mrs_natural) index_(i));
00512         pkViewOut(i, peakView::pkDeltaAmplitude, f) = deltamag_((mrs_natural) index_(i));
00513         pkViewOut(i, peakView::pkFrame, f) = frame_;
00514         pkViewOut(i, peakView::pkGroup, f) = 0.0; //This should be -1!!!! [TODO]
00515         pkViewOut(i, peakView::pkVolume, f) = 1.0;
00516         pkViewOut(i, peakView::pkBinLow, f) = interval_(2*i);
00517         pkViewOut(i, peakView::pkBin, f) = index_(i);
00518         pkViewOut(i, peakView::pkBinHigh, f) = interval_(2*i+1);
00519         pkViewOut(i, peakView::pkTrack, f) = -1.0; //null-track ID
00520 
00521         if(useStereoSpectrum_)
00522           pkViewOut(i, peakView::pkPan, f) = in((mrs_natural)index_(i)+2*N_, f);
00523         else
00524           pkViewOut(i, peakView::pkPan, f) = 0.0;
00525       }
00526     }
00527     else //if not yet reached "skip" number of frames...
00528     {
00529       for(mrs_natural i=0; i< frameMaxNumPeaks_; ++i)
00530       {
00531         //pkViewOut(i, peakView::pkFrequency, f) = 0;
00532         //pkViewOut(i, peakView::pkAmplitude, f) = 0;
00533         //pkViewOut(i, peakView::pkPhase, f) = 0;
00534         //pkViewOut(i, peakView::pkDeltaFrequency, f) = 0;
00535         //pkViewOut(i, peakView::pkDeltaAmplitude, f) = 0;
00536         pkViewOut(i, peakView::pkFrame, f) = frame_;
00537         //pkViewOut(i, peakView::pkGroup, f) = -1;
00538         //pkViewOut(i, peakView::pkVolume, f) = 0;
00539         //pkViewOut(i, peakView::pkPan, f) = 0;
00540         //pkViewOut(i, peakView::pkBinLow, f) = 0;
00541         //pkViewOut(i, peakView::pkBin, f) = 0;
00542         //pkViewOut(i, peakView::pkBinHigh, f) = 0;
00543       }
00544     }
00545     frame_++;
00546   }
00547 
00548   //count the total number of existing peaks (i.e. peak freq != 0)
00549   ctrl_totalNumPeaks_->setValue(pkViewOut.getTotalNumPeaks());
00550 
00551   // MATLAB_PUT(out, "peaks");
00552   // MATLAB_PUT(frameMaxNumPeaks_, "k");
00553   // MATLAB_EVAL("figure(1);clf;plot(peaks(6*k+1:7*k));");
00554 
00555 }