Marsyas
0.6.0-alpha
|
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