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