Marsyas  0.2
/home/gperciva/src/marsyas/src/marsyas/BeatHistoFeatures.cpp
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 
00020 #include "common.h"
00021 #include "BeatHistoFeatures.h"
00022 #include <algorithm>
00023 #include <iterator>
00024 #include <cfloat> 
00025 
00026 using std::ostringstream;
00027 using std::vector;
00028 using std::abs;
00029 
00030 using namespace Marsyas;
00031 
00032 
00033 //#define MTLB_DBG_LOG
00034 
00035 
00036 static void NormInPlace (realvec& vector)
00037 {
00038     mrs_real sum = vector.sum ();
00039     if (sum > 0)
00040         vector  /= sum;
00041     return;
00042 }
00043 
00044 static mrs_real PeriodicCentroid (realvec& vector, mrs_bool isLog = false, mrs_natural startIdx = 200)
00045 {
00046     mrs_natural len = vector.getCols ();
00047     mrs_real res = 0,
00048         norm = 0;
00049 
00050     for (mrs_natural i = startIdx; i < len; i++)
00051     {
00052         mrs_real theta = (isLog)? log(i*1./startIdx)* TWOPI :(i * TWOPI) / startIdx;
00053         res     += vector(i) * (.5*(cos(theta)+1));
00054         norm    += vector(i);
00055     }
00056 
00057     return res/norm;
00058 }
00059 
00060 
00061 static mrs_real PeriodicSpread (realvec& vector, mrs_real centroid, mrs_bool isLog = false, mrs_natural startIdx = 200)
00062 {
00063     mrs_natural len = vector.getCols ();
00064     mrs_real res = 0,
00065         norm = 0;
00066 
00067     for (mrs_natural i = startIdx; i < len; i++)
00068     {
00069         mrs_real theta = (isLog)? log(i*1./startIdx)* TWOPI :(i * TWOPI) / startIdx;
00070         res     += vector(i) * abs(.5*(cos(theta)+1)-centroid);
00071         norm    += vector(i);
00072     }
00073 
00074     return res/norm;
00075 }
00076 
00077 
00078 mrs_real BeatHistoFeatures::NumMax (mrs_realvec& vector)
00079 {
00080 
00081     pkr_->process(vector, pkres_);
00082 
00083     return pkres_.sum();
00084 }
00085 
00086 //static void BeatChroma (realvec& beatChroma, const realvec& beatHistogram, mrs_real beatRes = .25)
00087 //{
00088 //
00089 //  //
00090 //  mrs_natural i,j,k,
00091 //      len = beatHistogram.getCols ();
00092 //  const mrs_real startBpm     = 25;
00093 //  mrs_natural startIdx = (mrs_natural)(startBpm / beatRes + .5);
00094 //
00095 //  beatChroma.stretch (startIdx);      // length equals startIdx
00096 //
00097 //  for (j = 0; j < startIdx; j++)
00098 //  {
00099 //      mrs_real mean = 0;
00100 //      for (k = 0; k < 3; k++)
00101 //      {
00102 //          for (i = -k; i <= k; i++)
00103 //              mean    += beatHistogram((j+startIdx)*k+i);
00104 //
00105 //              beatChroma(j)   += mean/(k+1);
00106 //      }
00107 //  }
00108 //}
00109 
00110 static mrs_real SpectralFlatness (const realvec& beatHistogram, mrs_natural startIdx = 200)
00111 {
00112     mrs_real    res = 0;
00113     mrs_natural len = beatHistogram.getCols ();
00114     //mrs_natural sum = beatHistogram.sum ();
00115 
00116     //beatHistogram /= sum;
00117     for (mrs_natural i = startIdx; i < len; i++)
00118         res     += log(beatHistogram(i)+1e-6);
00119 
00120     return exp(res/(len-startIdx));
00121 }
00122 
00123 static mrs_real Std (const realvec& beatHistogram)
00124 {
00125     return sqrt (beatHistogram.var ());
00126 }
00127 
00128 static mrs_real MaxHps (const realvec& beatHistogram, mrs_natural startIdx = 200)
00129 {
00130     const mrs_natural order = 4;
00131     mrs_natural k,len = beatHistogram.getCols ();
00132     mrs_realvec res = beatHistogram;    // make this a member
00133 
00134     //res.setval(-1e38);
00135 
00136     for (k =2; k < order; k++)
00137     {
00138         for (mrs_natural i = startIdx; i < len; i++)
00139         {
00140             if (k*i >= len)
00141                 break;
00142             res(i)      += log(beatHistogram(k*i)+1e-6);
00143         }
00144     }
00145 
00146     for (k = 0; k < startIdx; k++)
00147         res(k)  = -1e38;
00148 
00149 
00150     return exp (res.maxval ());
00151 }
00152 
00153 static void MaxAcf (mrs_real& max, mrs_real& mean, const realvec& beatHistogram, realvec& res,mrs_natural startSearchAt, mrs_natural stopSearchAt)
00154 {
00155     mrs_natural k,len = beatHistogram.getCols ();
00156     
00157     res.setval(0.); 
00158 
00159     // compute ACF
00160     for (k = startSearchAt; k < stopSearchAt; k++) // this can be optimized
00161     {
00162         mrs_real val    = 0;
00163 
00164         for (mrs_natural i = k; i < len; i++)
00165         {
00166             val += beatHistogram(i) * beatHistogram(i-k);
00167         }
00168         
00169         
00170             res(k)  = val / (len-k);
00171     }
00172 
00173 
00174     //pkr_->process(in, pkres_);
00175 
00176     max = res.maxval ();
00177     
00178     mean = 1e6*res.mean ();
00179 }
00180 
00181 
00182 
00183 
00184 BeatHistoFeatures::BeatHistoFeatures(mrs_string name):MarSystem("BeatHistoFeatures", name)
00185 {
00186     mxr_ = NULL;
00187     pkr_ = NULL;
00188     pkr1_ = NULL;
00189     
00190     addControls();
00191 }
00192 
00193 BeatHistoFeatures::BeatHistoFeatures(const BeatHistoFeatures& a):MarSystem(a)
00194 {
00195     mxr_ = NULL;
00196     pkr_ = NULL;
00197     pkr1_ = NULL;
00198     ctrl_mode_ = getctrl("mrs_string/mode");
00199 }
00200 
00201 BeatHistoFeatures::~BeatHistoFeatures()
00202 {
00203   delete mxr_;
00204   delete pkr_;
00205   delete pkr1_;
00206   
00207 }
00208 
00209 
00210 MarSystem* 
00211 BeatHistoFeatures::clone() const 
00212 {
00213   return new BeatHistoFeatures(*this);
00214 }
00215 
00216 void 
00217 BeatHistoFeatures::addControls()
00218 {
00219     addctrl("mrs_string/mode", "method", ctrl_mode_);
00220 }
00221 
00222 void
00223 BeatHistoFeatures::myUpdate(MarControlPtr sender)
00224 {
00225     (void) sender;
00226     MRSDIAG("BeatHistoFeatures.cpp - BeatHistoFeatures:myUpdate");
00227   
00228     delete mxr_;//[!]
00229     delete pkr_;
00230     delete pkr1_;
00231     
00232     mxr_ = new MaxArgMax("mxr");//[!]
00233     pkr_ = new Peaker("pkr");
00234     pkr1_ = new Peaker("pkr1");
00235     
00236     
00237     setctrl("mrs_natural/onSamples", (mrs_natural)1);   
00238     setctrl("mrs_real/osrate", getctrl("mrs_real/israte"));
00239     
00240     mrs_string mode = ctrl_mode_->to<mrs_string>();
00241      
00242     setctrl("mrs_natural/onObservations", (mrs_natural)18);     // alex 
00243     setctrl("mrs_string/onObsNames", "BeatHisto_Sum,BeatHisto_LowPeakAmp,BeatHisto_MidPeakAmp,BeatHisto_HighPeakAmp,BeatHisto_LowBPM,BeatHisto_MidBPM,BeatHistoHighBPM,BeatHisto_LowMidBPMRatio,BeatHisto_MaxAcr,BeatHisto_MeanACR,BeatHisto_MaxHPS,BeatHisto_Flatness,BeatHisto_Std,BeatHisto_PeriodicCentroid1,BeatHisto_PeriodicCentroi2,BeatHisto_PeriodicSpread1,BeatHisto_PeriodicSpread2,BeatHisto_NumMax");
00244     
00245     flag_.create(getctrl("mrs_natural/inSamples")->to<mrs_natural>());
00246     mxr_->updControl("mrs_natural/inSamples", getctrl("mrs_natural/inSamples"));
00247     mxr_->updControl("mrs_natural/inObservations", getctrl("mrs_natural/inObservations"));
00248     mxr_->updControl("mrs_real/israte", getctrl("mrs_real/israte"));
00249     mxr_->updControl("mrs_natural/nMaximums", 3);
00250 
00251 
00252     pkr_->updControl("mrs_natural/inSamples", getctrl("mrs_natural/inSamples"));
00253     pkr_->updControl("mrs_natural/inObservations", getctrl("mrs_natural/inObservations"));
00254     pkr_->updControl("mrs_real/israte", getctrl("mrs_real/israte"));
00255 
00256 
00257     pkr1_->updControl("mrs_natural/inSamples", getctrl("mrs_natural/inSamples"));
00258     pkr1_->updControl("mrs_natural/inObservations", getctrl("mrs_natural/inObservations"));
00259     pkr1_->updControl("mrs_real/israte", getctrl("mrs_real/israte"));
00260 
00261 
00262     pkr1_->updControl("mrs_natural/peakNeighbors", 40);
00263     pkr1_->updControl("mrs_real/peakSpacing", 0.1);
00264     pkr1_->updControl("mrs_natural/peakStart", 200);
00265     pkr1_->updControl("mrs_natural/peakEnd", 640);
00266 
00267 
00268     pkr_->updControl("mrs_natural/peakNeighbors", 40);
00269     pkr_->updControl("mrs_real/peakSpacing", 0.1);
00270     pkr_->updControl("mrs_natural/peakStart", 200);
00271     pkr_->updControl("mrs_natural/peakEnd", 640);
00272      
00273     pkr_->setctrl("mrs_real/peakStrengthRelMax", 0.5);
00274     pkr_->setctrl("mrs_real/peakStrengthRelThresh", 2.0);
00275     pkr_->setctrl("mrs_real/peakStrengthThreshLpParam", 0.95);
00276     pkr_->setctrl("mrs_natural/peakNeighbors", 4);
00277 
00278 
00279     
00280     mxres_.create(mxr_->getctrl("mrs_natural/onObservations")->to<mrs_natural>(),
00281                   mxr_->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
00282     pkres_.create(pkr_->getctrl("mrs_natural/onObservations")->to<mrs_natural>(),
00283                   pkr_->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
00284 
00285     pkres1_.create(pkr1_->getctrl("mrs_natural/onObservations")->to<mrs_natural>(),
00286                    pkr1_->getctrl("mrs_natural/onSamples")->to<mrs_natural>());
00287     
00288 
00289 }
00290 
00291 mrs_real 
00292 BeatHistoFeatures::sum_nearby(mrs_natural index, mrs_natural radius, mrs_natural size, const realvec& in)
00293 {
00294     mrs_real sum = 0.0;
00295     mrs_natural ix;
00296     for (mrs_natural i = 1; i <= radius; ++i)
00297     {
00298         ix = index - i;
00299         if (0 < ix && ix < size) // Make sure we have a valid index.
00300             sum += in(0, ix);
00301 
00302         ix = index + i;
00303         if (0 < ix && ix < size)
00304             sum += in(0, ix);
00305     }        
00306 
00307     return sum;
00308 }
00309 
00310 // If the bins of <in> surrounding index <factor>*<tmx> have a sum greater
00311 // than <pmax>, set s1, s2, t1, t2. Otherwise do nothing.
00312 void 
00313 BeatHistoFeatures::harm_prob(mrs_real& pmax, mrs_real factor, 
00314                  mrs_real& s1, mrs_natural& t1, 
00315                  mrs_real& s2, mrs_natural& t2, 
00316                  mrs_natural tmx,
00317                  mrs_natural size, 
00318                  const realvec& in)
00319 {
00320 
00321     mrs_natural index = (mrs_natural) floor(factor * tmx + 0.5);    
00322     mrs_real c = (index > 100.0) ? 1.0 : 0.75;
00323     mrs_real a = ((50 < tmx) && (tmx < 100)) ? 1.5 : 0.75;
00324     mrs_real prob = 0.0;
00325 
00326     if (index < size)
00327     {
00328         prob = a * in(0,tmx) + c * in(0, index);
00329 
00330         // Decide how far to look based on index; in the past, values as high
00331         // as 9 have been used if index > 150.
00332         mrs_natural radius = index > 150 ? 6 : 3;
00333         prob += c * sum_nearby(index, radius, size, in);
00334     }
00335 
00336     if (prob > pmax)
00337     {
00338         pmax = prob;
00339         if (tmx < index)
00340         {
00341             s1 = in(0,tmx);
00342             s2 = in(0,index) + sum_nearby(index, 3, size, in);
00343             t1 = tmx+1;
00344         }
00345         else
00346         {
00347             s1 = in(0,index) + sum_nearby(index, 3, size, in);
00348             s2 = in(0,tmx);
00349             t1 = index+1;
00350         }
00351 
00352         t2 = (mrs_natural)(factor * t1);
00353     }  
00354 }
00355 
00356 
00357 
00358 // void 
00359 // BeatHistoFeatures::myProcess(realvec& in, realvec& out)
00360 // {
00361 //   // in.write("histo.plot");
00362 //   //checkFlow(in,out);
00363 //   mrs_natural o,c,t;
00364 //   mrs_real mx = DBL_MIN;
00365 //   mrs_natural tmx  = 0;
00366 //   mrs_real pmax = DBL_MIN;
00367 //   mrs_natural t1 = 0;
00368 //   mrs_natural t2 = 0;
00369 //   mrs_real s1 = 0.0;
00370 //   mrs_real s2 = 0.0;
00371 
00372 //   flag_.setval(0.0);
00373   
00374 //   // c, o, and t are declared in MarSystem.h. 
00375 //   // TODO: Why do we use c < 3 here? (i.e. why 3?)
00376 //   for (c=0; c < 3; ++c)
00377 //   {
00378 //       for (o=0; o < inObservations_; o++)
00379 //       {
00380 //           for (t = 0; t < inSamples_; t++)
00381 //           {
00382 //               if (((in(o,t) > mx) && (flag_(t) == 0.0)) && (40 < t) && (t < 120))
00383 //               {
00384 //                   mx = in(o,t);
00385 //                   tmx = t;
00386 //               }
00387 //           }
00388 //       }
00389 
00390 //       flag_(tmx) = 1.0;
00391 //       mx = DBL_MIN; // reset max
00392 
00393 //       if (tmx < 120.0)
00394 //       {
00395 //           harm_prob(pmax, 2, s1, t1, s2, t2, tmx, inSamples_, in);
00396 //           harm_prob(pmax, 3.0, s1, t1, s2, t2, tmx, inSamples_, in);
00397 //       }
00398 //       else 
00399 //       {
00400 //           harm_prob(pmax, 0.5, s1, t1, s2, t2, tmx, inSamples_, in);
00401 //           harm_prob(pmax, 0.33333, s1, t1, s2, t2, tmx, inSamples_, in);
00402 //       }
00403 //   }
00404   
00405 //   flag_.setval(0.0);
00406   
00407 //   mrs_real sum1 = 0.0;
00408 //   for (t = 40; t < 90; t++)
00409 //     sum1 += in(0,t);
00410 
00411 //   mrs_real sum2 = 0.0;
00412 //   for (t = 90; t < 140; t++)
00413 //     sum2 += in(0,t);
00414 
00415 //   mrs_real sum3 = 0.0;
00416 //   for (t = 40; t < 250; t++)
00417 //     sum3 += in(0,t);
00418     
00419 //   out(0,0) = s1;
00420 //   out(1,0) = t1;
00421 //   out(2,0) = s2;
00422 //   out(3,0) = t2;
00423 //   out(4,0) = (0 == t1) ? 0.0 : (t2 / t1);
00424 //   out(5,0) = sum1;
00425 //   out(6,0) = sum2;
00426 //   out(7,0) = sum3;
00427 // }
00428 
00429 
00430 
00431 
00432 
00433 
00434 void 
00435 BeatHistoFeatures::beatHistoFeatures(realvec& in, realvec& out)
00436 {
00437     
00438     mrs_real sum = 0;
00439   
00440     for (o=0; o < inObservations_; o++)
00441         for (t = 0; t < inSamples_; t++)
00442         {
00443             sum += in(o,t);
00444         }
00445 
00446     
00447     mrs_real result[2];
00448     mrs_natural i,startIdx = 200;
00449     // zero-out below 50BPM 
00450     for (i=0; i < startIdx; i++) 
00451         in(i) = 0;
00452 
00453     for (i = startIdx; i < in.getCols (); i++)
00454         if (in(i) < 0)
00455             in(i) = 0;
00456 
00457 
00458 
00459     
00460 
00461     pkr1_->process(in, pkres1_);
00462     mxr_->process(pkres1_,mxres_);
00463     
00464     
00465     
00466     vector<mrs_real> bpms;
00467     bpms.push_back(mxres_(0,1));
00468     bpms.push_back(mxres_(0,3));
00469     bpms.push_back(mxres_(0,5));
00470     
00471     sort(bpms.begin(), bpms.end());
00472 
00473     out(0,0) = sum; 
00474     for (unsigned int i=0; i<bpms.size(); i++) 
00475         for (unsigned int j =0; j < bpms.size(); j++)
00476         {
00477             if (bpms[i] == mxres_(0,2*j+1))
00478                 out(i+1,0) = mxres_(0,2*j);
00479         }
00480     
00481 
00482     
00483     out(4,0) = bpms[0] /4.0;
00484     out(5,0) = bpms[1] /4.0;
00485     out(6,0) = bpms[2] /4.0;
00486     out(7,0) = out(4,0) / out(5,0);
00487 
00488 
00489     
00490     NormInPlace (in);   
00491 
00492 
00493 
00494 #ifdef MARSYAS_MATLAB
00495 #ifdef MTLB_DBG_LOG
00496     MATLAB_PUT(in, "beathist");
00497     MATLAB_EVAL("figure(1);plot((201:800)/4, beathist(201:800)),grid on");
00498 #endif
00499 #endif
00500 
00501     MaxAcf (result[0], result[1],in, flag_, startIdx, 600);
00502     out(8,0)    = result[0];
00503     out(9,0)    = result[1];
00504     out(10,0)   = MaxHps (in, startIdx);
00505     out(11,0)    = SpectralFlatness (in, startIdx);
00506     out(12,0)   = Std(in);
00507     out(13,0)   = PeriodicCentroid(in, false, startIdx);
00508     out(14,0)   = PeriodicCentroid(in, true, startIdx);
00509     out(15,0)   = PeriodicSpread(in, out(13,0), false, startIdx);
00510     out(16,0)   = PeriodicSpread(in, out(14,0), true, startIdx);
00511     out(17,0)   = NumMax(in);   
00512     
00513 }
00514 
00515 
00516 void 
00517 BeatHistoFeatures::myProcess(realvec& in, realvec& out)
00518 {
00519     mrs_string mode = ctrl_mode_->to<mrs_string>();
00520     beatHistoFeatures(in,out);
00521 }
00522 
00523 
00524 
00525 
00526 
00527 
00528 
00529     
00530 
00531