Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/BeatHistoFeatures.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 
00020 #include "../common_source.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;  //suppress warning of unused parameter(s)
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 (mrs_natural o=0; o < inObservations_; o++)
00441     for (mrs_natural 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