Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/Peaker.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 "Peaker.h"
00020 #include "../common_source.h"
00021 #include <algorithm>
00022 
00023 using std::ostringstream;
00024 using std::min;
00025 using std::max;
00026 using std::cout;
00027 using std::endl;
00028 
00029 
00030 using namespace Marsyas;
00031 
00032 Peaker::Peaker(mrs_string name):MarSystem("Peaker",name)
00033 {
00034   addControls();
00035 }
00036 
00037 Peaker::~Peaker()
00038 {
00039 }
00040 
00041 MarSystem*
00042 Peaker::clone() const
00043 {
00044   return new Peaker(*this);
00045 }
00046 
00047 
00048 
00049 void
00050 Peaker::addControls()
00051 {
00052   addctrl("mrs_real/peakSpacing", 0.0);
00053   addctrl("mrs_real/peakStrength", 0.0);
00054   addctrl("mrs_real/peakStrengthRelMax", 0.0);
00055   addctrl("mrs_real/peakStrengthRelThresh", 0.0);
00056   addctrl("mrs_real/peakStrengthAbs", 0.0);
00057   addctrl("mrs_real/peakStrengthThreshLpParam", 0.95);
00058   addctrl("mrs_natural/peakStart", (mrs_natural)0);
00059   addctrl("mrs_natural/peakEnd", (mrs_natural)0);
00060   addctrl("mrs_natural/interpolation", (mrs_natural)0);
00061   //addctrl("mrs_real/peakGain", 1.0);
00062   addctrl("mrs_bool/peakHarmonics", false);
00063   addctrl("mrs_bool/rmsNormalize", false);
00064   addctrl("mrs_natural/peakNeighbors", 2);
00065 }
00066 
00067 
00068 void
00069 Peaker::myUpdate(MarControlPtr sender)
00070 {
00071   (void) sender;  //suppress warning of unused parameter(s)
00072   MRSDIAG("Peaker.cpp - Peaker:myUpdate");
00073 
00074   MarSystem::myUpdate (sender);
00075 
00076   lpThresh_.stretch (getctrl ("mrs_natural/inSamples")->to<mrs_natural>());
00077 }
00078 
00079 void
00080 Peaker::myProcess(realvec& in, realvec& out)
00081 {
00082   mrs_natural t,o;
00083 
00084   const mrs_natural peakSpacing  = (mrs_natural)(inSamples_ *getctrl("mrs_real/peakSpacing")->to<mrs_real>() + .5);
00085   mrs_real peakStrengthRelRms,
00086            peakStrengthRelMax,
00087            peakStrengthRelThresh,
00088            peakStrengthAbs;
00089   //mrs_real peakGain;
00090   mrs_bool peakHarmonics;
00091   mrs_bool rmsNormalize;
00092 
00093   mrs_natural peakStart;
00094   mrs_natural peakEnd;
00095   mrs_natural interpolationMode;
00096   mrs_natural peakNeighbors;
00097 
00098 
00099   peakStrengthRelRms = getctrl("mrs_real/peakStrength")->to<mrs_real>();
00100   peakStrengthRelMax = getctrl("mrs_real/peakStrengthRelMax")->to<mrs_real>();
00101   peakStrengthRelThresh = getctrl("mrs_real/peakStrengthRelThresh")->to<mrs_real>();
00102   lpCoeff_ = getctrl("mrs_real/peakStrengthThreshLpParam")->to<mrs_real>();
00103   peakStrengthAbs = getctrl("mrs_real/peakStrengthAbs")->to<mrs_real>();
00104   peakStart = getctrl("mrs_natural/peakStart")->to<mrs_natural>();
00105   peakEnd = getctrl("mrs_natural/peakEnd")->to<mrs_natural>();
00106   interpolationMode = getctrl("mrs_natural/interpolation")->to<mrs_natural>();
00107   //peakGain = getctrl("mrs_real/peakGain")->to<mrs_real>();
00108   peakHarmonics = getctrl("mrs_bool/peakHarmonics")->to<mrs_bool>();
00109   rmsNormalize = getctrl("mrs_bool/rmsNormalize")->to<mrs_bool>();
00110   peakNeighbors = getctrl("mrs_natural/peakNeighbors")->to<mrs_natural>();
00111 
00112 
00113 
00114   if (peakEnd == 0)
00115     peakEnd = inSamples_;
00116   // FIXME This line defines an unused variable
00117   // mrs_real srate = getctrl("mrs_real/israte")->to<mrs_real>();
00118 
00119   out.setval(0.0);
00120 
00121 
00122   //peakStrengthRelRms = 0.0;
00123 
00124 
00125 
00126   for (o = 0; o < inObservations_; o++)
00127   {
00128     rms_    = 0.0;
00129     max_    = -1e37;
00130 
00131     for (t=peakStart; t < peakEnd; t++)
00132     {
00133       rms_ += in(o,t) * in(o,t);
00134       if (max_ < in(o,t))
00135         max_    = in(o,t);
00136     }
00137     if (rms_ != 0.0)
00138       rms_ /= (peakEnd - peakStart);
00139     rms_ = sqrt(rms_);
00140 
00141     mrs_real max;
00142     mrs_natural maxIndex;
00143 
00144     bool peakFound = false;
00145 
00146     if (peakStrengthRelThresh > .0)
00147     {
00148       in.getRow (o,lpThresh_);
00149       compLpThresh (lpThresh_, lpThresh_);  // do it inplace to avoid another copy...
00150     }
00151 
00152     for (t=peakStart; t < peakEnd; t++)
00153     {
00154       peakFound = true;
00155 
00156       // peak has to be larger than neighbors
00157       for (int j = 1; j < peakNeighbors; j++)
00158       {
00159         mrs_natural index=t-j;
00160         if (index<0) index=0;
00161         if (in(o,index) >= in(o,t))
00162         {
00163           peakFound = false;
00164           break;
00165         }
00166         index=t+j;
00167         if (index>=inSamples_) index=inSamples_-1;
00168         if (in(o,index) >= in(o,t))
00169         {
00170           peakFound = false;
00171           break;
00172         }
00173       }
00174 
00175       if (peakFound)
00176       {
00177         currThresh_ = lpThresh_(t);
00178         peakFound   = doThresholding (in(o,t), peakStrengthRelRms, peakStrengthRelMax, peakStrengthRelThresh, peakStrengthAbs);
00179       }
00180 
00181       if (peakFound)
00182       {
00183         // check for another peak in the peakSpacing area
00184         max = in(o,t);
00185         maxIndex = t;
00186 
00187 
00188         for (int j=0; j < peakSpacing; j++)
00189         {
00190           if (t+j < peakEnd-1)
00191             if (in(o,t+j) > max)
00192             {
00193               max = in(o,t+j);
00194               maxIndex = t+j;
00195             }
00196         }
00197 
00198         t += peakSpacing;
00199 
00200         if (rmsNormalize)
00201         {
00202           out(o,maxIndex) = in(o,maxIndex) / rms_;
00203           if(interpolationMode && maxIndex > 0 && maxIndex < inSamples_)
00204           {
00205             out(o,maxIndex-1) = in(o,maxIndex-1) /rms_;
00206             out(o,maxIndex+1) = in(o,maxIndex+1) / rms_;
00207           }
00208         }
00209         else
00210         {
00211           out(o,maxIndex) = in(o,maxIndex);
00212           if(interpolationMode && maxIndex > 0 && maxIndex < inSamples_)
00213           {
00214             out(o,maxIndex-1) = in(o,maxIndex-1);
00215             out(o,maxIndex+1) = in(o,maxIndex+1);
00216           }
00217         }
00218 
00219 
00220         // peakNeighbors = 0;
00221         // rms_ = 1.0;
00222 
00223 
00224         if (peakHarmonics)
00225         {
00226           twice_ = 2 * maxIndex;
00227           half_ = (mrs_natural) (0.5 * maxIndex + 0.5);
00228           triple_ = 3 * maxIndex;
00229           third_ = (mrs_natural) (0.33 * maxIndex + 0.5);
00230 
00231           if (twice_ < (peakEnd - peakStart))
00232           {
00233             peakFound = true;
00234 
00235             // peak has to be larger than neighbors
00236             for (int j = 1; j < peakNeighbors; j++)
00237             {
00238               if (in(o,twice_-j) >= in(o,twice_))
00239               {
00240                 peakFound = false;
00241                 break;
00242               }
00243               if (in(o,twice_+j) >= in(o,twice_))
00244               {
00245                 peakFound = false;
00246                 break;
00247               }
00248             }
00249 
00250 
00251             if (peakFound)
00252             {
00253               out(o, maxIndex) += (in(o, twice_)/rms_);
00254               out(o, twice_) = in(o,twice_)/rms_ + 0.5 * out(o, maxIndex);
00255             }
00256 
00257           }
00258 
00259           if (half_ < (peakEnd - peakStart))
00260           {
00261             peakFound = true;
00262 
00263             // peak has to be larger than neighbors
00264             for (int j = 1; j < peakNeighbors; j++)
00265             {
00266               if (in(o,half_-j) >= in(o,half_))
00267               {
00268                 peakFound = false;
00269                 break;
00270               }
00271               if (in(o,half_+j) >= in(o,half_))
00272               {
00273                 peakFound = false;
00274                 break;
00275               }
00276             }
00277 
00278 
00279             if (peakFound)
00280             {
00281               out(o, maxIndex) += (in(o, half_)/rms_);
00282               out(o, half_) = in(o,half_)/rms_ + 0.5 * out(o, maxIndex);
00283             }
00284 
00285           }
00286 
00287 
00288           if (triple_ < (peakEnd - peakStart))
00289           {
00290             peakFound = true;
00291 
00292             // peak has to be larger than neighbors
00293             for (int j = 1; j < peakNeighbors; j++)
00294             {
00295               if (in(o,triple_-j) >= in(o,triple_))
00296               {
00297                 peakFound = false;
00298                 break;
00299               }
00300               if (in(o,triple_+j) >= in(o,triple_))
00301               {
00302                 peakFound = false;
00303                 break;
00304               }
00305             }
00306 
00307 
00308             if (peakFound)
00309             {
00310               out(o, maxIndex) += (in(o, triple_)/rms_);
00311               out(o, triple_) = in(o,triple_)/rms_ + 0.5 * out(o, maxIndex);
00312             }
00313 
00314           }
00315 
00316           if (third_ < (peakEnd - peakStart))
00317           {
00318             peakFound = true;
00319 
00320             // peak has to be larger than neighbors
00321 
00322             for (int j = 1; j < peakNeighbors; j++)
00323             {
00324               if (in(o,third_-j) >= in(o,third_))
00325               {
00326                 peakFound = false;
00327                 break;
00328               }
00329               if (in(o,third_+j) >= in(o,triple_))
00330               {
00331                 peakFound = false;
00332                 break;
00333               }
00334             }
00335 
00336             if (peakFound)
00337             {
00338               out(o, maxIndex) += (in(o, third_)/rms_);
00339               out(o, third_) = in(o,third_)/rms_ + 0.5 * out(o, maxIndex);
00340             }
00341 
00342           }
00343         }
00344 
00345         peakFound = true;
00346       }
00347 
00348     }
00349   }
00350 }
00351 
00352 mrs_bool Peaker::doThresholding (mrs_real input, mrs_real peakStrengthRelRms, mrs_real peakStrengthRelMax, mrs_real peakStrengthRelThresh, mrs_real peakStrengthAbs)
00353 {
00354   mrs_real thresh = max (max (max (peakStrengthAbs, peakStrengthRelRms * rms_), peakStrengthRelMax * max_), peakStrengthRelThresh * currThresh_);
00355 
00356 
00357   if (input <= thresh)
00358     return false;
00359   else
00360     return true;
00361 }
00362 
00363 void Peaker::compLpThresh (const mrs_realvec input, mrs_realvec &output)
00364 {
00365   mrs_natural i,
00366               len = input.getCols ();
00367   mrs_real  buf = input(0);
00368   for (i = 0; i < len; ++i)
00369   {
00370     buf         = input(i) + lpCoeff_ * (buf - input(i));
00371     output(i)   = buf;
00372   }
00373   // do reverse filtering to ensure zero group delay
00374   for (i = len-1; i >= 0; i--)
00375   {
00376     buf         = output(i) + lpCoeff_ * (buf - output(i));
00377     output(i)   = buf;
00378   }
00379 }