Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/PhaseLock.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2006 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 "PhaseLock.h"
00020 #include "../common_source.h"
00021 
00022 #include <marsyas/FileName.h>
00023 //#include "../common_source.h"
00024 
00025 #define MINIMUMREAL 0.000001 //(0.000001 minimum float recognized)
00026 #define NA -10000.0 //undefined value flag
00027 
00028 using namespace std;
00029 using namespace Marsyas;
00030 
00031 PhaseLock::PhaseLock(mrs_string name):MarSystem("PhaseLock", name)
00032 {
00033   addControls();
00034   timeElapsed_ = 0;
00035   lastGTBeatTime_ = -1;
00036   inductionFinished_ = false;
00037   gtAfter2ndBeat_ = false;
00038   triggerInduction_ = false;
00039   lastGTIBI_ = 0.0;
00040 }
00041 
00042 PhaseLock::PhaseLock(const PhaseLock& a) : MarSystem(a)
00043 {
00044   // For any MarControlPtr in a MarSystem
00045   // it is necessary to perform this getctrl
00046   // in the copy constructor in order for cloning to work
00047   ctrl_beatHypotheses_ = getctrl("mrs_realvec/beatHypotheses");
00048   ctrl_inductionTime_ = getctrl("mrs_natural/inductionTime");
00049   ctrl_nrPeriodHyps_ = getctrl("mrs_natural/nrPeriodHyps");
00050   ctrl_nrPhasesPerPeriod_ = getctrl("mrs_natural/nrPhasesPerPeriod");
00051   ctrl_scoreFunc_ = getctrl("mrs_string/scoreFunc");
00052   ctrl_gtBeatsFile_ = getctrl("mrs_string/gtBeatsFile");
00053   ctrl_hopSize_ = getctrl("mrs_natural/hopSize");
00054   ctrl_srcFs_ = getctrl("mrs_real/srcFs");
00055   ctrl_mode_ = getctrl("mrs_string/mode");
00056   ctrl_backtrace_ = getctrl("mrs_bool/backtrace");
00057   ctrl_tickCount_ = getctrl("mrs_natural/tickCount");
00058   ctrl_innerMargin_ = getctrl("mrs_real/innerMargin");
00059   ctrl_lftOutterMargin_ = getctrl("mrs_real/lftOutterMargin");
00060   ctrl_rgtOutterMargin_ = getctrl("mrs_real/rgtOutterMargin");
00061   ctrl_corFactor_ = getctrl("mrs_real/corFactor");
00062   ctrl_maxPeriod_ = getctrl("mrs_natural/maxPeriod");
00063   ctrl_minPeriod_ = getctrl("mrs_natural/minPeriod");
00064   ctrl_adjustment_ = getctrl("mrs_natural/adjustment");
00065   ctrl_dumbInduction_ = getctrl("mrs_bool/dumbInduction");
00066   ctrl_inductionOut_ = getctrl("mrs_string/inductionOut");
00067   ctrl_triggerInduction_ = getctrl("mrs_bool/triggerInduction");
00068   ctrl_curBestScore_ = getctrl("mrs_real/curBestScore");
00069   ctrl_triggerBestScoreFactor_ = getctrl("mrs_real/triggerBestScoreFactor");
00070 
00071   inductionFinished_ = a.inductionFinished_;
00072   gtInitPhase_ = a.gtInitPhase_;
00073   gtScore_ = a.gtScore_;
00074   gtInitPeriod_ = a.gtInitPeriod_;
00075   gtLastPeriod_ = a.gtLastPeriod_;
00076   srcFs_ = a.srcFs_;
00077   hopSize_ = a.hopSize_;
00078   adjustment_ = a.adjustment_;
00079   gtAfter2ndBeat_ = a.gtAfter2ndBeat_;
00080   triggerInduction_ = a.triggerInduction_;
00081   mode_ = a.mode_;
00082   curBestScore_ = a.curBestScore_;
00083   lastGTBeatTime_ = a.lastGTBeatTime_;
00084   triggerBestScoreFactor_ = a.triggerBestScoreFactor_;
00085   lastGTIBI_ = a.lastGTIBI_;
00086 }
00087 
00088 PhaseLock::~PhaseLock()
00089 {
00090 }
00091 
00092 MarSystem*
00093 PhaseLock::clone() const
00094 {
00095   return new PhaseLock(*this);
00096 }
00097 
00098 void
00099 PhaseLock::addControls()
00100 {
00101   //Add specific controls needed by this MarSystem.
00102   addctrl("mrs_realvec/beatHypotheses", realvec(), ctrl_beatHypotheses_);
00103   addctrl("mrs_natural/inductionTime", -1, ctrl_inductionTime_);
00104   setctrlState("mrs_natural/inductionTime", true);
00105   addctrl("mrs_natural/nrPeriodHyps", 6, ctrl_nrPeriodHyps_);
00106   setctrlState("mrs_natural/nrPeriodHyps", true);
00107   addctrl("mrs_natural/nrPhasesPerPeriod", 20, ctrl_nrPhasesPerPeriod_);
00108   setctrlState("mrs_natural/nrPhasesPerPeriod", true);
00109   addctrl("mrs_string/scoreFunc", "regular", ctrl_scoreFunc_);
00110   setctrlState("mrs_string/scoreFunc", true);
00111   addctrl("mrs_natural/hopSize", -1, ctrl_hopSize_);
00112   setctrlState("mrs_natural/hopSize", true);
00113   addctrl("mrs_real/srcFs", -1.0, ctrl_srcFs_);
00114   setctrlState("mrs_real/srcFs", true);
00115   addctrl("mrs_string/gtBeatsFile", "input.txt", ctrl_gtBeatsFile_);
00116   addctrl("mrs_string/mode", "regular", ctrl_mode_);
00117   setctrlState("mrs_string/mode", true);
00118   addctrl("mrs_bool/backtrace", false, ctrl_backtrace_);
00119   setctrlState("mrs_bool/backtrace", true);
00120   addctrl("mrs_natural/tickCount", 0, ctrl_tickCount_);
00121   addctrl("mrs_real/innerMargin", 3.0, ctrl_innerMargin_);
00122   setctrlState("mrs_real/innerMargin", true);
00123   addctrl("mrs_real/lftOutterMargin", 0.2, ctrl_lftOutterMargin_);
00124   setctrlState("mrs_real/lftOutterMargin", true);
00125   addctrl("mrs_real/rgtOutterMargin", 0.4, ctrl_rgtOutterMargin_);
00126   setctrlState("mrs_real/rgtOutterMargin", true);
00127   addctrl("mrs_real/corFactor", 0.5, ctrl_corFactor_);
00128   setctrlState("mrs_real/corFactor", true);
00129   addctrl("mrs_natural/maxPeriod", -1, ctrl_maxPeriod_);
00130   setctrlState("mrs_natural/maxPeriod", true);
00131   addctrl("mrs_natural/minPeriod", -1, ctrl_minPeriod_);
00132   setctrlState("mrs_natural/minPeriod", true);
00133   addctrl("mrs_natural/adjustment", 0, ctrl_adjustment_);
00134   setctrlState("mrs_natural/adjustment", true);
00135   addctrl("mrs_bool/dumbInduction", false, ctrl_dumbInduction_);
00136   addctrl("mrs_string/inductionOut", "-1", ctrl_inductionOut_);
00137   setctrlState("mrs_string/inductionOut", true);
00138   addctrl("mrs_bool/triggerInduction", false, ctrl_triggerInduction_);
00139   setctrlState("mrs_bool/triggerInduction", true);
00140   addctrl("mrs_real/curBestScore", NA, ctrl_curBestScore_);
00141   setctrlState("mrs_real/curBestScore", true);
00142   addctrl("mrs_real/triggerBestScoreFactor", 0.95, ctrl_triggerBestScoreFactor_);
00143 }
00144 
00145 void
00146 PhaseLock::myUpdate(MarControlPtr sender)
00147 {
00148   (void) sender;  //suppress warning of unused parameter(s)
00149   MRSDIAG("PhaseLock.cpp - PhaseLock:myUpdate");
00150 
00151   inductionTime_ = ctrl_inductionTime_->to<mrs_natural>();
00152   nrPeriodHyps_ = ctrl_nrPeriodHyps_->to<mrs_natural>();
00153   nrPhasesPerPeriod_ = ctrl_nrPhasesPerPeriod_->to<mrs_natural>();
00154   scoreFunc_ = ctrl_scoreFunc_->to<mrs_string>();
00155   mode_ = ctrl_mode_->to<mrs_string>();
00156   backtrace_ = ctrl_backtrace_->to<mrs_bool>();
00157   innerMargin_ = (mrs_natural) ctrl_innerMargin_->to<mrs_real>();
00158   lftOutterMargin_ = ctrl_lftOutterMargin_->to<mrs_real>();
00159   rgtOutterMargin_ = ctrl_rgtOutterMargin_->to<mrs_real>();
00160   corFactor_ = ctrl_corFactor_->to<mrs_real>();
00161   maxPeriod_ = ctrl_maxPeriod_->to<mrs_natural>();
00162   minPeriod_ = ctrl_minPeriod_->to<mrs_natural>();
00163   hopSize_ = ctrl_hopSize_->to<mrs_natural>();
00164   srcFs_ = ctrl_srcFs_->to<mrs_real>();
00165   adjustment_ = (mrs_real) ctrl_adjustment_->to<mrs_natural>();
00166   inductionOut_ = ctrl_inductionOut_->to<mrs_string>();
00167   triggerInduction_ = ctrl_triggerInduction_->to<mrs_bool>();
00168   curBestScore_ = ctrl_curBestScore_->to<mrs_real>();
00169   triggerBestScoreFactor_ = ctrl_triggerBestScoreFactor_->to<mrs_real>();
00170 
00171   ctrl_onSamples_->setValue(3, NOUPDATE);
00172   ctrl_onObservations_->setValue(nrPeriodHyps_, NOUPDATE);
00173   ctrl_osrate_->setValue(ctrl_israte_, NOUPDATE);
00174 
00175   nInitHyp_ = nrPeriodHyps_ * nrPhasesPerPeriod_; //Nr. of initial hypothesis
00176 
00177   hypSignals_.create(nInitHyp_, inSamples_); //N hypothesis signals of induction timming size.
00178   firstBeatPoint_.create(nInitHyp_); //for accumulating reversed phases (one for each hypothesis pair)
00179   trackingScore_.create(nInitHyp_);
00180   beatCount_.create(nInitHyp_);
00181   maxLocalTrackingScore_.create(nrPeriodHyps_);
00182   maxLocalTrackingScoreInd_.create(nrPeriodHyps_);
00183   initPhases_.create(nrPeriodHyps_);
00184   lastPhases_.create(nrPeriodHyps_);
00185   initPeriods_.create(nrPeriodHyps_);
00186   lastPeriods_.create(nrPeriodHyps_);
00187   metricalSalience_.create(nrPeriodHyps_);
00188   rawScore_.create(nrPeriodHyps_);
00189   rawScoreNorm_.create(nrPeriodHyps_);
00190   metricalRelScore_.create(nrPeriodHyps_);
00191   scoreNorm_.create(nrPeriodHyps_);
00192 
00193   // reset maxLoclTrackingScores and indexes to allow negative scores
00194   for (int i = 0; i < nrPeriodHyps_; i++)
00195   {
00196     maxLocalTrackingScore_(i) = NA;
00197     maxLocalTrackingScoreInd_(i) = -1;
00198   }
00199 }
00200 
00201 mrs_realvec
00202 PhaseLock::GTInitialization(realvec& in, realvec& out, mrs_natural gtInitPhase, mrs_natural gtInitPeriod)
00203 {
00204   (void) out;
00205   //MATLAB_PUT(in, "Flux_FlowThrued");
00206 
00207   //cout << "InitPhase: " << gtInitPhase << "(" << (((gtInitPhase * hopSize_)-hopSize_/2) / srcFs_) << "); initPeriod: " << gtInitPeriod << endl;
00208 
00209   mrs_realvec initHypothesis(5); //backtracedPhase, initPeriod, lastPhase, lastPeriod, initScore
00210   for(int i = 0; i < initHypothesis.getCols(); i++)
00211     initHypothesis(i) = 0; //clean array
00212 
00213   //Calculate score inside induction window (simulate tracking behaviour) given the ground-truth (two beats)
00214   //adjust initPhase around innerMargin
00215   mrs_real localPeakAmp = in(gtInitPhase);
00216   mrs_natural localPeak = gtInitPhase;
00217   //cout << "gtInitPhase: " << gtInitPhase << "; value: " << localPeakAmp << endl;
00218   for (int t = gtInitPhase-innerMargin_; t <= gtInitPhase+innerMargin_; t++)
00219   {
00220     //Limit analysis to induction window
00221     if(t >= (inSamples_-1-inductionTime_))
00222     {
00223       if(in(t) > localPeakAmp)
00224       {
00225         localPeakAmp = in(t);
00226         localPeak = t;
00227       }
00228     }
00229   }
00230   //cout << "2- adjInitPhase: " << localPeak << "; adjValue: " << localPeakAmp << " period: " << gtInitPeriod << endl;
00231   initHypothesis(0) = localPeak; //keep first beat point (adjusted to maximum wihtin innerMargin)
00232   initHypothesis(1) = gtInitPeriod; //keep initPeriod
00233   initHypothesis(4) += (localPeakAmp * (gtInitPeriod/(mrs_real)maxPeriod_)); //assume perfect initial beat prediction
00234 
00235   outterWinLft_ = (mrs_natural) ceil(gtInitPeriod * lftOutterMargin_); //(% of IBI)
00236   outterWinRgt_ = (mrs_natural) ceil(gtInitPeriod * rgtOutterMargin_); //(% of IBI)
00237   mrs_natural beatPoint = gtInitPhase;
00238   mrs_real scoreFraction = 0.0;
00239   mrs_natural error = 0;
00240   mrs_real errorIn = MINIMUMREAL; //to avoid divide by zero
00241   mrs_natural period = gtInitPeriod;
00242   mrs_natural periodTmp = 0;
00243   //mrs_natural it=1; //used for cout debugging
00244   do
00245   {
00246     //cout << "3- beatPoint: " << beatPoint << " period: " << period << endl;
00247 
00248     periodTmp = period + ((mrs_natural) ((errorIn*corFactor_) + ((errorIn/abs(errorIn)) * 0.5)));
00249     //if the updated period don't respect the limits then keeps the former value
00250     if(periodTmp >= minPeriod_ && periodTmp <= maxPeriod_)
00251       period = periodTmp;
00252 
00253     //(for safechecking)
00254     if((beatPoint + (period + (mrs_natural) ((errorIn*corFactor_) + ((errorIn/abs(errorIn)) * 0.5)))) >= inSamples_)
00255     {
00256       initHypothesis(2) = beatPoint; //Keep current lastPhase
00257       initHypothesis(3) = period; //Keep current lastPeriod
00258       break;
00259     }
00260 
00261     beatPoint += (period + (mrs_natural) ((errorIn*corFactor_) + ((errorIn/abs(errorIn)) * 0.5)));
00262     errorIn = MINIMUMREAL; //to avoid divide by zero
00263 
00264     initHypothesis(2) = beatPoint; //Keep current lastPhase
00265     initHypothesis(3) = period; //Keep current lastPeriod
00266     //Check maximum around outter window
00267     localPeakAmp = in(beatPoint);
00268     localPeak = beatPoint;
00269     //cout << "predicition: " << localPeak << "; predValue: " << localPeakAmp << "; period: " << period << endl;
00270     for (int t = beatPoint-outterWinLft_; t <= beatPoint+outterWinRgt_; t++)
00271     {
00272       if(t > (inSamples_-1-inductionTime_) && t < inSamples_)
00273       {
00274         //cout << "t: " << t << "; in(t): " << in(t) << endl;
00275         if(in(t) > localPeakAmp)
00276         {
00277           localPeakAmp = in(t);
00278           localPeak = t;
00279 
00280           //cout << "Max: " << localPeak << "; MaxV: " << localPeakAmp << endl;
00281         }
00282       }
00283     }
00284 
00285     //Check if maximum is inside inner window
00286     //calc positive regular score and adjust hypothesis
00287     error = localPeak-beatPoint;
00288     if(localPeak >= beatPoint-innerMargin_ && localPeak <= beatPoint+innerMargin_)
00289     {
00290       //only consider error in innerWindow for adjustment
00291       errorIn = error;
00292       scoreFraction = (mrs_real) ((mrs_real)abs(error) / outterWinRgt_);
00293       initHypothesis(4) += (((1 - scoreFraction) * localPeakAmp) * (period/(mrs_real)maxPeriod_)); //update score
00294 
00295       if(errorIn == 0)
00296         errorIn = MINIMUMREAL;
00297     }
00298     //if outside inner window calc negative regular score
00299     else
00300     {
00301       scoreFraction = (mrs_real) ((mrs_real)abs(error)/ outterWinRgt_);
00302       //cout << "SC: " << initHypothesis(4) << "; V: " << localPeakAmp << "; sf: " << scoreFraction << endl;
00303       initHypothesis(4) += ((-1 * scoreFraction * localPeakAmp) * (period/(mrs_real)maxPeriod_));  //update score
00304     }
00305 
00306     //cout << "-it:" << it << "; Error: " << error << "; Max: " << localPeak << "; MaxV: " << localPeakAmp
00307     //  << "; Sf: " << scoreFraction << "(" << initHypothesis(4) << "); BeatPoint: " << beatPoint << "; period: " << period << endl;
00308     //it++;
00309 
00310     //don't evaluate beat if the outterWin around it surprasses the induction window
00311   } while(beatPoint < inSamples_);
00312 
00313   return initHypothesis;
00314 }
00315 
00316 //read beat-times from groundtruth file
00317 mrs_realvec
00318 PhaseLock::readGTFile(mrs_string gtFilePath)
00319 {
00320   mrs_realvec gtHypotheses(4); //save gt hypotheses retrieved from gt file
00321   mrs_real gtBeatTimeInit1 = NA, gtBeatTimeInit2 = NA, gtBeatTimeLast1 = NA, gtBeatTimeLast2 = NA;
00322 
00323   //cout << "\nINDUCTION TIME: " << timeElapsed_ << "(" << (((timeElapsed_* hopSize_)-hopSize_/2) / srcFs_) << ") [" << inductionTime_ << "]" << endl;
00324   inStream_.open(gtFilePath.c_str());
00325 
00326   getline (inStream_, line_);
00327 
00328   //for beat groundtruth files (in line, separated by spaces):
00329   mrs_natural pos0, pos1, pos2;
00330   pos0 = (mrs_natural) line_.find_first_of(" ", 0); //initial delimiter
00331   pos1 = (mrs_natural) line_.find_first_of(" ", pos0+1); //second delimiter
00332   pos2 = (mrs_natural) line_.find_first_of(" ", pos1+1); //third delimiter
00333 
00334   if(pos1 >= 0) //if beat times column by column (on first line)
00335   {
00336     //check gtfile for eof at first run -> save last beat time
00337     if(lastGTBeatTime_ < 0)
00338     {
00339       std::istringstream iss(line_);
00340       char c[10]; //big enough array
00341       while (iss >> c) //space ("") delimiter
00342       {
00343         //last c contains last value in file
00344         lastGTBeatTime_ = atof(c); //save end value
00345       }
00346       iss.clear();
00347       //cout << "last GT beat time: " << lastGTBeatTime_ << endl;
00348     }
00349 
00350     //INIT PHASE & PERIOD
00351     //consider initial search point as inductionTime before the current time point
00352     mrs_real initBeatPoint = (((mrs_real)((timeElapsed_-inductionTime_) * hopSize_)) - adjustment_) / srcFs_;
00353     //cout << "InitSearchPoint: " << initBeatPoint << endl;
00354     //if more than two beatTimes given in the annotation file
00355     //discart initial beatTime (due to inconsistencies in the beggining of some annotation files)
00356     mrs_real lastGTBeatTimeInit1 = 0.0;
00357     //mrs_real lastGTBeatTimeLast1 = 0.0;
00358     do
00359     {
00360       //cout << "gtBeatTimeInit1: " << gtBeatTimeInit1 << "; gtBeatTimeInit2:"
00361       //        << gtBeatTimeInit2 << "; initBeatPoint: " << initBeatPoint << endl;
00362 
00363       //gtBeatTimeInit1 == lastGTBeatTime => reached end of gt file =>
00364       //=> propagate values to current position based on last given hypotheses in the GT file
00365       if((gtBeatTimeInit1 == lastGTBeatTime_) && (gtBeatTimeInit1 < initBeatPoint))
00366       {
00367         mrs_real initGTIBI = abs(gtBeatTimeInit1 - lastGTBeatTimeInit1);
00368         mrs_natural k = (mrs_natural) ceil((initBeatPoint+MINIMUMREAL - gtBeatTimeInit1) / initGTIBI);
00369         gtBeatTimeInit1 = gtBeatTimeInit1 + (k * initGTIBI);
00370         gtBeatTimeInit2 = gtBeatTimeInit2 + initGTIBI;
00371 
00372         break;
00373       }
00374       else
00375       {
00376         //gtBeatTimeInit1 > initBeatPoint => reached gt beat time corresponding to current time point => ok
00377         if(gtBeatTimeInit1 > initBeatPoint)
00378           break;
00379         else
00380         {
00381           lastGTBeatTimeInit1 = gtBeatTimeInit1; //save last gtBeatTimeInit1
00382 
00383           gtBeatTimeInit1 = strtod(line_.substr(pos0+1, pos1).c_str(), NULL);
00384           gtBeatTimeInit2 = strtod(line_.substr(pos1+1, pos2).c_str(), NULL);
00385 
00386           //keep on looking
00387           pos0 = (mrs_natural) line_.find_first_of(" ", pos0+1); //initial delimiter
00388           pos1 = (mrs_natural) line_.find_first_of(" ", pos0+1); //second delimiter
00389         }
00390       }
00391 
00392     } while(gtBeatTimeInit1 <= initBeatPoint);
00393 
00394     //cout << "INIT-> gtBeatTime1: " << gtBeatTimeInit1 << "(" << pos0 << "); gtBeatTime2: " << gtBeatTimeInit2
00395     //  << "(" << pos1 << "); initBP: " << initBeatPoint << endl;
00396 
00397     //LAST PHASE & PERIOD
00398     mrs_real lastBeatPoint = (((mrs_real)(timeElapsed_ * hopSize_)) - adjustment_) / srcFs_;
00399     if((gtBeatTimeInit2 < lastGTBeatTime_) && (gtBeatTimeInit2 > initBeatPoint)) //if still not reached end of file
00400     {
00401       do
00402       {
00403         //if still looking but beat2 already reached end of file =>
00404         //assign beat1 as last beat in file (beat2) and beat2 one period (based on last ibi) above
00405         if((gtBeatTimeLast2 == lastGTBeatTime_) && (gtBeatTimeLast1 < lastBeatPoint))
00406         {
00407           //cout << "it2: beat2 reached eof: gtBeatTimeLast2: " << gtBeatTimeLast2 << "; gtBeatTimeLast1: " << gtBeatTimeLast1
00408           //    << "; lastGTBeatTime: " << lastGTBeatTime_ << "; lastBeatPoint: " << lastBeatPoint
00409           //    << "; ibi: " << lastGTIBI_ << endl;
00410 
00411           lastGTIBI_ = abs(gtBeatTimeLast2 - gtBeatTimeLast1);
00412           //lastGTBeatTimeLast1 = gtBeatTimeLast1;
00413           gtBeatTimeLast1 = gtBeatTimeLast2;
00414           gtBeatTimeLast2 = gtBeatTimeLast1 + lastGTIBI_;
00415         }
00416 
00417         //gtBeatTimeLast1 == lastGTBeatTime => reached end of gt file =>
00418         //=> propagate values to current position based on last given hypotheses in the GT file
00419         if((gtBeatTimeLast1 == lastGTBeatTime_) && (gtBeatTimeLast1 < lastBeatPoint))
00420         {
00421           //cout << "it2: beat1 reached eof: gtBeatTimeLast1: " << gtBeatTimeLast1 << "(" << lastGTBeatTimeLast1 << ")"
00422           //    << "; lastGTBeatTime: " << lastGTBeatTime_ << "; lastBeatPoint: " << lastBeatPoint << "; ibi: " << lastGTIBI_ << endl;
00423 
00424           mrs_natural k = (mrs_natural) ceil((lastBeatPoint+MINIMUMREAL - gtBeatTimeLast1) / lastGTIBI_);
00425           gtBeatTimeLast1 = gtBeatTimeLast1 + (k * lastGTIBI_);
00426           gtBeatTimeLast2 = gtBeatTimeLast1 + lastGTIBI_;
00427 
00428           break;
00429         }
00430         else
00431         {
00432           //gtBeatTimeLast1 > initBeatPoint => reached gt beat time corresponding to current time point => ok
00433           if(gtBeatTimeLast1 > lastBeatPoint)
00434             break;
00435           else
00436           {
00437             //lastGTBeatTimeLast1 = gtBeatTimeLast1; //save last gtBeatTimeLast1
00438 
00439             //keep on looking
00440             pos0 = (mrs_natural) line_.find_first_of(" ", pos0+1); //initial delimiter
00441             pos1 = (mrs_natural) line_.find_first_of(" ", pos0+1); //second delimiter
00442 
00443             gtBeatTimeLast1 = strtod(line_.substr(pos0+1, pos1).c_str(), NULL);
00444             gtBeatTimeLast2 = strtod(line_.substr(pos1+1, pos2).c_str(), NULL);
00445 
00446             if(gtBeatTimeLast2 > gtBeatTimeLast1)
00447               lastGTIBI_ = abs(gtBeatTimeLast2 - gtBeatTimeLast1);
00448 
00449             //cout << "gtBeatTimeLast1: " << gtBeatTimeLast1 << "; gtBeatTimeLast2:"
00450             //  << gtBeatTimeLast2 << "; lastBeatPoint: " << lastBeatPoint << "; ibi: " << lastGTIBI_ << endl;
00451           }
00452         }
00453 
00454       } while(gtBeatTimeLast1 <= lastBeatPoint);
00455     }
00456     else //if already reached end of gt file on the first request (beat-time at the begining of this induction window)
00457     {
00458       //calculate gtBeatTimeLast1 and gtBeatTimeLast2 from init values
00459       cerr << "Reached end of ground-truth file...Last GT values propagated from the last hypotheses given by the GT file!" << endl;
00460 
00461       mrs_natural k = (mrs_natural) ceil((lastBeatPoint+MINIMUMREAL - gtBeatTimeInit1) / lastGTIBI_);
00462       gtBeatTimeLast1 = gtBeatTimeInit1 + (k * lastGTIBI_);
00463       gtBeatTimeLast2 = gtBeatTimeLast1 + lastGTIBI_;
00464     }
00465 
00466     //cout << "LAST-> gtBeatTime1: " << gtBeatTimeLast1 << "; gtBeatTime2: " << gtBeatTimeLast2 << endl;
00467 
00468     //cout <<  "gtInitBeatTime1: " << gtBeatTimeInit1 << "; gtInitBeatTime2: " << gtBeatTimeInit2
00469     //  << "; gtLastBeatTime1: " << gtBeatTimeLast1 << "; gtLastBeatTime2: " << gtBeatTimeLast2
00470     //  << "; initBeat: " << initBeatPoint << "; lastBeat: " << lastBeatPoint <<"; lastGtBeat: "
00471     //  << lastGTBeatTime_ << endl;
00472   }
00473 
00474   gtHypotheses(0) = gtBeatTimeInit1;
00475   gtHypotheses(1) = gtBeatTimeInit2;
00476   gtHypotheses(2) = gtBeatTimeLast1;
00477   gtHypotheses(3) = gtBeatTimeLast2;
00478 
00479   return gtHypotheses;
00480 }
00481 
00482 //testing trigger ground-truth without calculating pseudo-tracking starting from gt values
00483 void
00484 PhaseLock::handleGTHypotheses(realvec& in, realvec& out,mrs_string gtFilePath, mrs_realvec gtHypotheses)
00485 {
00486   mrs_real gtBeatTimeInit1 = gtHypotheses(0);
00487   mrs_real gtBeatTimeInit2 = gtHypotheses(1);
00488   mrs_real gtBeatTimeLast1 = gtHypotheses(2);
00489   mrs_real gtBeatTimeLast2 = gtHypotheses(3);
00490 
00491   //gt init period and phase in frames
00492   gtInitPeriod_ = (mrs_natural) (((abs(gtBeatTimeInit2 - gtBeatTimeInit1) * srcFs_) / hopSize_) + 0.5);
00493   gtInitPhase_ = (mrs_natural) ceil(((gtBeatTimeInit1 * srcFs_) + adjustment_) / hopSize_);
00494   //gtInitPhase_ = (mrs_natural) ceil(((gtBeatTimeIni1 * srcFs_)) / hopSize_);
00495 
00496   //gt last period and phase in frames
00497   gtLastPeriod_ = (mrs_natural) (((abs(gtBeatTimeLast2 - gtBeatTimeLast1) * srcFs_) / hopSize_) + 0.5);
00498   gtLastPhase_ = (mrs_natural) ceil(((gtBeatTimeLast1 * srcFs_) + adjustment_) / hopSize_);
00499   //gtLastPhase_ = (mrs_natural) ceil(((gtBeatTimeLast1 * srcFs_)) / hopSize_);
00500 
00501   //cout << "iniBeat1: " << gtBeatTimeInit1 << "; iniBeat2: " << gtBeatTimeInit2 << "; lasBeat1: " << gtBeatTimeLast1
00502   //    << "; lasBeat2: "<< gtBeatTimeLast2 << "; gtIniPh: " << gtInitPhase_ << "; gtIniP: " << gtInitPeriod_
00503   //    << "(" << (60/(gtBeatTimeInit2 - gtBeatTimeInit1)) << "BPM)(" << (gtBeatTimeInit2 - gtBeatTimeInit1) << "); gtLastPh: "
00504   //    << gtLastPhase_ << "; gtLastP: " << gtLastPeriod_ << "(" << (60/(gtBeatTimeLast2 - gtBeatTimeLast1)) << "BPM)" << endl;
00505 
00506   //ADJUST PHASE TO INDUCTION (ACCUMULATOR) WINDOW
00507   mrs_natural phaseAccAdjust = (inSamples_ - 1 - timeElapsed_) + gtInitPhase_;
00508   //cout << "; gtPAdj: " << phaseAccAdjust << "t: " << timeElapsed_ << "; inS: " << inSamples_ << endl;
00509 
00510   //calculate k for reversing phase till minimum (by subtracting period multiples)
00511   mrs_natural k = (mrs_natural) (((mrs_real)((inSamples_-1-inductionTime_)-phaseAccAdjust)/(mrs_real)gtInitPeriod_)+MINIMUMREAL);
00512   k += 1; //start on 2nd beat prediction (due to unconsistencies on the very beggining of the accumulator)
00513   phaseAccAdjust = ((mrs_natural) phaseAccAdjust) + ((mrs_natural) (k * gtInitPeriod_));
00514 
00515   if(strcmp(mode_.c_str(), "2b2") == 0 || strcmp(mode_.c_str(), "2b") == 0) //adjust values to induction window
00516   {
00517     mrs_realvec gtInitHypothesis = GTInitialization(in, out, phaseAccAdjust, gtInitPeriod_);
00518 
00519     if(curBestScore_ == NA) //if initial induction window
00520     {
00521       gtScore_ = gtInitHypothesis(4);
00522       //cout << "InitGtScore: " << gtScore_ << endl;
00523     }
00524     else //if trigger induction
00525     {
00526       if(curBestScore_ > 0) //to avoid best score inversion
00527         gtScore_ = curBestScore_ * triggerBestScoreFactor_;
00528       else
00529         gtScore_ = curBestScore_ / triggerBestScoreFactor_;
00530 
00531       //cout << "t: " << timeElapsed_ << "; curBest: " << curBestScore_ << "; gtScore: " << gtScore_ << endl;
00532     }
00533 
00534     if(strcmp(mode_.c_str(), "2b") == 0) //adjust phases and periods
00535     {
00536       gtInitPeriod_ = (mrs_natural) gtInitHypothesis(1); //initPeriod
00537       gtLastPeriod_ = (mrs_natural) gtInitHypothesis(3); //lastPeriod
00538       gtInitPhase_ = (mrs_natural) gtInitHypothesis(0) - (inSamples_-1 - timeElapsed_);
00539       gtLastPhase_ = (mrs_natural) gtInitHypothesis(2) - (inSamples_-1 - timeElapsed_) + gtLastPeriod_;
00540 
00541       //cout << "real gtInitPhase: " << gtInitPhase_ << "; gtLastPhase: " << gtLastPhase_ << endl;
00542 
00543       //for safechecking (really unlikely to happen) if lastPhases still don't surprasses the induction window then sum other period
00544       if(gtLastPhase_ <= timeElapsed_)
00545         gtLastPhase_ += gtLastPeriod_;
00546     }
00547 
00548     //cout << "iniPh: " << gtInitPhase_ << "(" << (((gtInitPhase_ * hopSize_) - adjustment_) / srcFs_) << "); iniPer: " << gtInitPeriod_
00549     //  << "; lasPh: " << gtLastPhase_ << "(" << (((gtLastPhase_ * hopSize_) - adjustment_) / srcFs_) << "); lasPer: " << gtLastPeriod_
00550     //  << "; gtSc: " << gtScore_ << endl;
00551 
00552     cerr << "\nInduction replaced by IBI given by two beats from ground-truth file at: " << gtFilePath << endl;
00553   }
00554 
00555   //simulate tracking under induction windows, given first beat (or any beat)
00556   //retrieve initial score, backtraced phase (subtracting period multiples), and phase just after the induction
00557   //(ground-truth initial phase and period as calculated through ACF (normal induction))
00558   if(strcmp(mode_.c_str(), "1b") == 0 || strcmp(mode_.c_str(), "1b1") == 0)
00559   {
00560     //matrix with all the N generated beat hypotheses in the induction stage
00561     beatHypotheses_ = ctrl_beatHypotheses_->to<mrs_realvec>();
00562 
00563     mrs_realvec lastBeatPoint(nInitHyp_);
00564     mrs_realvec lastLocalPeriod(nInitHyp_);
00565     for(int i = 0; i < nrPeriodHyps_; i++)
00566     {
00567       mrs_natural phase = phaseAccAdjust;
00568       mrs_natural period = (mrs_natural) beatHypotheses_(i*nrPhasesPerPeriod_, 0);
00569 
00570       initPeriods_(i) = period; //keep each init period hypothesis
00571       if(period > 1) //(phase is necessarily valid)
00572       {
00573         mrs_realvec gtInitHypothesis = GTInitialization(in, out, phase, period);
00574         lastPeriods_(i) = (mrs_natural) gtInitHypothesis(3); //keep each lastPeriod
00575 
00576         if(strcmp(mode_.c_str(), "1b") == 0) //adjusted phases to induction window
00577         {
00578           initPhases_(i) = (mrs_natural) gtInitHypothesis(0) - (inSamples_-1 - timeElapsed_);
00579           lastPhases_(i) = (mrs_natural) gtInitHypothesis(2) - (inSamples_-1 - timeElapsed_) + lastPeriods_(i);
00580 
00581           //for safechecking (really unlikely to happen) if lastPhases still don't surprasses the induction window then sum other period
00582           if(lastPhases_(i) <= timeElapsed_)
00583             lastPhases_(i) += lastPeriods_(i);
00584         }
00585 
00586         trackingScore_(i) = gtInitHypothesis(4);
00587 
00588         rawScore_(i) = trackingScore_(i); //rawScore equals simulated trackingScore
00589         maxLocalTrackingScore_(i) = trackingScore_(i); //every hypothesis is maximum local because there's only one phase per period
00590 
00591         //cout << "i-" << i << "; iniPhase: " << initPhases_(i) << "; initPeriod: " << initPeriods_(i)
00592         //  << "; lastPhase: " << lastPhases_(i) << "; lastPeriod: " << lastPeriods_(i)
00593         //  << "; rawScore: " << rawScore_(i) << endl;
00594       }
00595     }
00596 
00597     mrs_real maxMetricalRelScore = NA; //to allow best negative score
00598     mrs_natural avgPeriod = 0;
00599     mrs_real avgLocalTrackingScore_ = 0.0;
00600     mrs_real maxGlobalTrackingScore_ = NA; //to allow best negative score
00601     //mrs_natural maxMetricalRelScoreInd = -1;
00602     for(int i = 0; i < nrPeriodHyps_; i++)
00603     {
00604       if(initPeriods_(i) > MINIMUMREAL) //if the given period_ > 0 (= valid period_) (initial phase, being ground-truth, is necessarily valid
00605       {
00606         metricalRelScore_(i) = calcRelationalScore(i, rawScore_);
00607 
00608         if(metricalRelScore_(i) > maxMetricalRelScore)
00609         {
00610           maxMetricalRelScore = metricalRelScore_(i);
00611           //maxMetricalRelScoreInd = i;
00612         }
00613 
00614         if(maxLocalTrackingScore_(i) > maxGlobalTrackingScore_)
00615           maxGlobalTrackingScore_ = maxLocalTrackingScore_(i);
00616 
00617         avgLocalTrackingScore_ += maxLocalTrackingScore_(i);
00618         avgPeriod += (mrs_natural) initPeriods_(i);
00619 
00620         //cout << "MetrRelScore" << i << ":" << metricalRelScore_(i) << "; iPeriod:" << initPeriods_(i)
00621         //  << "; iPhase:" << initPhases_(i) << "(" << beatHypotheses_((mrs_natural)maxLocalTrackingScoreInd_(i), 1)
00622         //  << "); MaxMetrRelScore:" << maxMetricalRelScore << "(" << maxMetricalRelScoreInd << ")" << endl;
00623 
00624       }
00625     }
00626 
00627     for(int i = 0; i < nrPeriodHyps_; i++)
00628     {
00629       if(strcmp(mode_.c_str(), "1b") == 0) //adjusted phase to induction window
00630       {
00631         if(backtrace_)
00632         {
00633           //Period:
00634           out(i, 0) = initPeriods_(i);
00635           //Phase:
00636           out(i, 1) = initPhases_(i);
00637 
00638           if(i == 0)
00639           {
00640             if(gtAfter2ndBeat_)
00641               cerr << "\nInitial phase backtraced from second beat in ground-truth file at: " << gtFilePath << endl;
00642             else
00643               cerr << "\nInitial phase backtraced from first beat in ground-truth file at: " << gtFilePath << endl;
00644 
00645             cerr << "Initial phase: " << (((initPhases_(i) * hopSize_) - adjustment_) / srcFs_) << "s" << endl;
00646           }
00647         }
00648         else //if no backtrace
00649         {
00650           //Period:
00651           out(i, 0) = lastPeriods_(i);
00652           //Phase:
00653           out(i, 1) = lastPhases_(i);
00654 
00655           if(i == 0)
00656           {
00657             if(gtAfter2ndBeat_)
00658               cerr << "\nInitial phase adjusted from second beat in ground-truth file at: " << gtFilePath << endl;
00659             else
00660               cerr << "\nInitial phase adjusted from first beat in ground-truth file at: " << gtFilePath << endl;
00661 
00662             cerr << "Initial phase: " << (((lastPhases_(i) * hopSize_) - adjustment_) / srcFs_) << "s" << endl;
00663           }
00664         }
00665       }
00666       else if(strcmp(mode_.c_str(), "1b1") == 0) //phase as given by gt file
00667       {
00668         if(backtrace_)
00669         {
00670           //Period:
00671           out(i, 0) = initPeriods_(i);
00672           //Phase:
00673           out(i, 1) = gtInitPhase_; //keep first phase as extracted from gt file
00674 
00675           if(i == 0)
00676           {
00677             if(gtAfter2ndBeat_)
00678               cerr << "\nInitial phase replaced by second beat in ground-truth file at: " << gtFilePath << endl;
00679             else
00680               cerr << "\nInitial phase replaced by first beat in ground-truth file at: " << gtFilePath << endl;
00681 
00682             cerr << "Initial phase: " << (((gtInitPhase_ * hopSize_) - adjustment_) / srcFs_) << "s" << endl;
00683           }
00684         }
00685         else //if no backtrace (causal)
00686         {
00687           //Period:
00688           out(i, 0) = lastPeriods_(i);
00689           //Phase:
00690           out(i, 1) = gtLastPhase_; //keep first phase as extracted from gt file
00691 
00692           if(i == 0)
00693           {
00694             cerr << "\nInitial phase replaced by first beat in ground-truth file at: " << gtFilePath << endl;
00695             cerr << "Initial phase: " << (((gtLastPhase_ * hopSize_) - adjustment_) / srcFs_) << "s" << endl;
00696           }
00697         }
00698       }
00699 
00700       //After-Induction Score Calulcation: ===========================================================
00701       mrs_real finalScore; //final score
00702       if(!dumbInduction_) //if not in dumbInduction consider metrical relations
00703       {
00704         mrs_real metricalRelFraction = 0.0;
00705         metricalRelFraction = (metricalRelScore_(i) / maxMetricalRelScore);
00706 
00707         //to avoid negative score inversions
00708         if(metricalRelScore_(i) < 0 && maxMetricalRelScore > 0 && maxGlobalTrackingScore_ > 0
00709             && abs(metricalRelScore_(i)) > maxMetricalRelScore)
00710           metricalRelFraction = -1;
00711         if(metricalRelScore_(i) < 0 && maxMetricalRelScore < 0 && maxGlobalTrackingScore_ > 0)
00712           metricalRelFraction = (maxMetricalRelScore / metricalRelScore_(i));
00713 
00714 
00715         if(curBestScore_ == NA) //initially consider a proportion of the maximum tracking score for this induction window
00716           finalScore = metricalRelFraction * maxGlobalTrackingScore_;
00717         else //at each trigger consider a proportion of the current best score
00718         {
00719           if(curBestScore_ > 0) //to avoid best score inversion
00720             finalScore = (metricalRelFraction * curBestScore_) * triggerBestScoreFactor_;
00721           else
00722             finalScore = (metricalRelFraction * curBestScore_) / triggerBestScoreFactor_;
00723 
00724           //cout << "relFrac: " << metricalRelFraction << "; curBestScore: " << curBestScore_ << "; finalScore: " << finalScore << endl;
00725         }
00726 
00727         //cout << i << ": " << finalScore << "; MetrScore: " << metricalRelScore_(i) << "; MaxMetrRelScore: "
00728         //  << maxMetricalRelScore << "; MaxCorrScore: " << maxGlobalTrackingScore_ << "; initPeriod: " << initPeriods_(i)
00729         //  << "; initPhase: " << initPhases_(i) << "; lastPeriod: " << lastPeriods_(i) << "; lastPhase: "
00730         //  << lastPhases_(i) << endl;
00731 
00732         //cout << "scoreFunc:" << scoreFunc_ << "; avgMaxSum:" << avgMaxSum << "; avgPeriod:" << avgPeriod << endl;
00733       }
00734       else //if dumbInduction (either manual or needed-no salient periodicities) don't consider metrical relations
00735       {
00736         //cout << "DUMB INDUCTION" << endl;
00737 
00738         if(curBestScore_ == NA) //initially consider a proportion of the maximum tracking score for this induction window
00739           finalScore = rawScore_(i);
00740         else //at each trigger consider a proportion of the current best score
00741         {
00742           //final score = a fraction of the curBestScore given by the proportion of this agent's induction score
00743           //to the maximum induction score, and by a user-defined triggerBestScoreFactor
00744 
00745           //to avoid rawScoreFraction negative inversions
00746           mrs_real rawScoreFraction;
00747           //(to avoid divide by 0 if hypothesis exists)
00748           if(initPeriods_(i) > MINIMUMREAL && initPhases_(i) > MINIMUMREAL)
00749           {
00750             if(rawScore_(i) == 0.0) rawScore_(i) = MINIMUMREAL;
00751             if(maxGlobalTrackingScore_ == 0.0) maxGlobalTrackingScore_ = MINIMUMREAL;
00752           }
00753 
00754           //if rawScore negative and maxRawScore positive && rawScore or maxRawScore < 1
00755           if((rawScore_(i) < 0 && maxGlobalTrackingScore_ > 0) && (fabs(rawScore_(i)) < 1 || maxGlobalTrackingScore_ < 1))
00756             rawScoreFraction = rawScore_(i) * maxGlobalTrackingScore_;
00757           else
00758           {
00759             if(fabs(rawScore_(i)) < maxGlobalTrackingScore_) //if abs(rawScore) < maxRawScore
00760               rawScoreFraction = rawScore_(i) / maxGlobalTrackingScore_;
00761             else //if abs(rawScore) > maxRawScore
00762               rawScoreFraction = maxGlobalTrackingScore_ / rawScore_(i);
00763           }
00764 
00765           if(curBestScore_ > 0) //to avoid best score inversion
00766             finalScore = (rawScoreFraction * curBestScore_) * triggerBestScoreFactor_;
00767           else
00768             finalScore = (rawScoreFraction * curBestScore_) / triggerBestScoreFactor_;
00769 
00770           //cout << "rawSc: " << rawScore_(i) << "; maxRawSc: " << maxGlobalTrackingScore_ << "; rawFr: " << rawScoreFraction
00771           //    << "; curBestScore: " << curBestScore_ << "; finalScore: " << finalScore << "; triggerScFact: "
00772           //    << triggerBestScoreFactor_ << endl;
00773         }
00774 
00775         //reset dumb induction
00776         ctrl_dumbInduction_->setValue(false);
00777       }
00778 
00779       //aditional initial score normalization dependent of score function (due to their relative weights)
00780       if(strcmp(scoreFunc_.c_str(), "correlation") == 0 || strcmp(scoreFunc_.c_str(), "squareCorr") == 0 )
00781         finalScore *= 5; //("correlation" and "squareCorr" are, in average, 5times more reactive than "regular")
00782 
00783       out(i, 2) = finalScore;
00784     }
00785   }
00786 
00787   //cout << "Beat1: " << gtBeatTime1_ << "(" << gtPhase_ << "-" << gtInitPhase_ << "); Beat2: "
00788   //    << gtBeatTime2_ << "; Period: " << gtInitPeriod_ << "("
00789   //    << (((gtBeatTime2_ - gtBeatTime1_) * srcFs_) / hopSize_) << "); Score:" << gtScore_ << endl;
00790 
00791   inStream_.close();
00792 }
00793 
00794 mrs_natural
00795 PhaseLock::metricalRelation(mrs_real period1, mrs_real period2)
00796 {
00797   mrs_real periodA = (period1 > period2 ? period1 : period2);
00798   mrs_real periodB = (period1 > period2 ? period2 : period1);
00799   mrs_real multiple = periodA / periodB;
00800   mrs_real tolerance = 0.15; //4/periodB;
00801 
00802   mrs_natural floorMultiple = (mrs_natural) floor(multiple);
00803   mrs_natural ceilMultiple = (mrs_natural) ceil(multiple);
00804   mrs_natural rel = 0;
00805 
00806   //if(print)
00807   //cout << "Per1:" << period1 << "; Per2:" << period2 << "; Multiple:" << multiple
00808   //<< "(" << (multiple - tolerance) << ":" << (multiple + tolerance) << ")" << endl;
00809 
00810   //if periods are integerly related (with a given tolerance)
00811   if(multiple - tolerance <= floorMultiple)
00812   {
00813     if(floorMultiple == 2 || floorMultiple == 4)
00814       rel = 6 - floorMultiple; //then retrieve integer relation
00815     else if(floorMultiple == 3)
00816       rel = 3; //then retrieve integer relation
00817     else if(floorMultiple >= 5 && floorMultiple <= 8)
00818       rel = 1;
00819     else
00820       rel = 0;
00821   }
00822 
00823   else if(multiple + tolerance >= ceilMultiple)
00824   {
00825     if(ceilMultiple == 2 || ceilMultiple == 4)
00826       rel = 6 - ceilMultiple; //then retrieve integer relation
00827     else if(ceilMultiple == 3)
00828       rel = 3; //then retrieve integer relation
00829     else if(ceilMultiple >= 5 && ceilMultiple <= 8)
00830       rel = 1;
00831     else rel = 0;
00832   }
00833 
00834   else rel = 0; //else retrieve 0 (periods non metricaly related)
00835 
00836   return rel;
00837 }
00838 
00839 mrs_real
00840 PhaseLock::calcRelationalScore(mrs_natural i, mrs_realvec rawScoreVec)
00841 {
00842   mrs_real score = 0.0;
00843   score = 2*5*rawScoreVec(i); //self-relation = 10*
00844 
00845   mrs_natural metricalRel;
00846   for(int j = 0; j < nrPeriodHyps_; j++)
00847   {
00848     if(j != i && initPeriods_(i) > MINIMUMREAL && initPeriods_(j) > MINIMUMREAL)
00849     {
00850       metricalRel = metricalRelation((mrs_natural)initPeriods_(i), (mrs_natural)initPeriods_(j));
00851       score += rawScoreVec(j) * metricalRel;
00852     }
00853   }
00854 
00855   return score;
00856 }
00857 
00858 
00859 void
00860 PhaseLock::regularFunc(realvec& in, realvec& out)
00861 {
00862   //matrix with all the N generated beat hypotheses in the induction stage
00863   beatHypotheses_ = ctrl_beatHypotheses_->to<mrs_realvec>();
00864 
00865   //MATLAB_PUT(beatHypotheses_, "BeatHypotheses");
00866   //MATLAB_PUT(in, "Flux_FlowThrued2");
00867 
00868   //Build N hypotheses signals (phase + k*periods):
00869   mrs_realvec lastBeatPoint(nInitHyp_);
00870   mrs_realvec lastLocalPeriod(nInitHyp_);
00871   for(int h = 0; h < nInitHyp_; h++)
00872   {
00873     mrs_natural phase = (mrs_natural) beatHypotheses_(h, 1);
00874     mrs_natural period = (mrs_natural) beatHypotheses_(h, 0);
00875     if(phase > 0 && period > 1)
00876     {
00877       //if(timeElapsed_ == 1228)
00878       //cout << "h: " << h << "; ph: " << phase << "; per: " << period << endl;
00879 
00880       //calculate k for reversing phase till minimum (by subtracting period multiples)
00881       mrs_natural k = (mrs_natural) (((mrs_real)((inSamples_-1-inductionTime_)-phase)/(mrs_real)period)+MINIMUMREAL);
00882 
00883       if((phase + ((mrs_natural) ((k+1) * period))) < inSamples_) //to avoid surpassing accumulator
00884         k += 1; //start on 2nd beat prediction (due to unconsistencies on the very beggining of the accumulator)
00885 
00886       firstBeatPoint_(h) = phase + ((mrs_natural) (k * period));
00887 
00888       //cout << "H: " << h << "-1BP: " << firstBeatPoint_(h) << endl;
00889 
00890       //adjust initPhase around innerMargin
00891       beatCount_(h) = 1;
00892       mrs_real localPeakAmp = in((mrs_natural)firstBeatPoint_(h));
00893       mrs_natural localPeak = (mrs_natural) firstBeatPoint_(h);
00894       //mrs_natural errorTmp = 0;
00895       for (int t = (mrs_natural)firstBeatPoint_(h)-innerMargin_; t <= (mrs_natural)firstBeatPoint_(h)+innerMargin_; t++)
00896       {
00897         //Limit analysis to induction window
00898         if(t >= (inSamples_-1-inductionTime_))
00899         {
00900           if(in(t) > localPeakAmp)
00901           {
00902             //cout << "localPeak: " << localPeak << "; Amp: " << localPeakAmp << "; newPeak: " << t << "; newAmp: " << in(t) << endl;
00903             localPeakAmp = in(t);
00904             localPeak = t;
00905           }
00906         }
00907       }
00908       //errorTmp = (mrs_natural) firstBeatPoint_(h) - localPeak;
00909       //cout << "INIT_ERROR: " << errorTmp << endl;
00910       //cout << "Adj: " << localPeak << "; REALinitPhase: " << localPeak - (inSamples_-1 - timeElapsed_) << endl;
00911 
00912       trackingScore_(h) += localPeakAmp * ((mrs_real)period / (mrs_real)maxPeriod_); //Assume perfect initial beat prediction
00913 
00914       //if(timeElapsed_ == 1228)
00915       //cout << "TrackingScore1: " << trackingScore_(h) << " [" << localPeakAmp * ((mrs_real)period / (mrs_real)maxPeriod_) << "]" << endl;
00916 
00917       //cout << "h:" << h << "-k:" << k << "(" << (((mrs_real)((inSamples_-1-inductionTime_)-phase)/(mrs_real)period)+MINIMUMREAL) << ")"
00918       //    << ";period: " << beatHypotheses_(h, 0) << ";phase: " << firstBeatPoint_(h) << "(" << (mrs_natural) beatHypotheses_(h, 1) << ")"
00919       //    << ";adPhase: " << localPeak << ";S: " << trackingScore_(h) << endl;
00920 
00921       firstBeatPoint_(h) = localPeak; //Keep first beat point (adjusted to maximum wihtin innerMargin)
00922 
00923       //cout << "FirstBeat: " << firstBeatPoint_(h) << endl;
00924       hypSignals_(h, (mrs_natural) firstBeatPoint_(h)) = 1.0;
00925 
00926       outterWinLft_ = (mrs_natural) ceil(period * lftOutterMargin_); //(% of IBI)
00927       outterWinRgt_ = (mrs_natural) ceil(period * rgtOutterMargin_); //(% of IBI)
00928       mrs_natural beatPoint = (mrs_natural) firstBeatPoint_(h);
00929       mrs_real scoreFraction = 0.0;
00930       mrs_natural error = 0;
00931       mrs_real errorIn = MINIMUMREAL; //to avoid divide by zero
00932       mrs_natural periodTmp = 0;
00933       //mrs_natural it=1; //used for cout debugging
00934       do
00935       {
00936         //cout << "Period: " << period << endl;
00937         //cout << "Hypothesis (" << h << "): Phase: " << phase << " period_: " << period << "; corFactor: " << corFactor_ << endl;
00938         periodTmp = period + ((mrs_natural) ((errorIn*corFactor_) + ((errorIn/abs(errorIn)) * 0.5)));
00939 
00940         //cout << "NewPeriod: " << periodTmp << "; errorIn: " << errorIn << "; adj: " << ((errorIn*corFactor_) + ((errorIn/abs(errorIn)) * 0.5)) << endl;
00941         //if the updated period don't respect the limits then keeps the former value
00942         if(periodTmp >= minPeriod_ && periodTmp <= maxPeriod_)
00943           period = periodTmp;
00944 
00945         //(for safechecking)
00946         if((beatPoint + (period + (mrs_natural) ((errorIn*corFactor_) + ((errorIn/abs(errorIn)) * 0.5)))) >= inSamples_)
00947         {
00948           lastBeatPoint(h) = beatPoint; //Keep current lastPhase
00949           lastLocalPeriod(h) = period; //Keep current lastPeriod
00950           break;
00951         }
00952 
00953         beatPoint += (period + (mrs_natural) ((errorIn*corFactor_) + ((errorIn/abs(errorIn)) * 0.5)));
00954         errorIn = MINIMUMREAL; //to avoid divide by zero
00955 
00956         //cout << "curBeatPoint: " << beatPoint << endl;
00957 
00958         lastBeatPoint(h) = beatPoint;
00959         lastLocalPeriod(h) = period;
00960         //Check maximum around outter window
00961         localPeakAmp = in(beatPoint);
00962         localPeak = beatPoint;
00963         for (int t = beatPoint-outterWinLft_; t <= beatPoint+outterWinRgt_; t++)
00964         {
00965           if(t > (inSamples_-1-inductionTime_) && t < inSamples_)
00966           {
00967             if(in(t) > localPeakAmp)
00968             {
00969               localPeakAmp = in(t);
00970               localPeak = t;
00971             }
00972           }
00973         }
00974 
00975         //Check if maximum is inside inner window
00976         //calc positive regular score and adjust hypothesis
00977         error = localPeak-beatPoint;
00978         if(localPeak >= beatPoint-innerMargin_ && localPeak <= beatPoint+innerMargin_)
00979         {
00980           //only consider error in innerWindow for adjustment
00981           errorIn = error;
00982           scoreFraction = (mrs_real) ((mrs_real)abs(error) / outterWinRgt_);
00983           trackingScore_(h) += (((1 - scoreFraction) * localPeakAmp) * ((mrs_real)period / (mrs_real)maxPeriod_));
00984 
00985           if(errorIn == 0)
00986             errorIn = MINIMUMREAL;
00987 
00988           //cout << "SCORE_FRACTION: " << (((1 - scoreFraction) * localPeakAmp) * ((mrs_real)period / (mrs_real)maxPeriod_)) << endl;
00989         }
00990         //if outside inner window calc negative regular score
00991         else
00992         {
00993           scoreFraction = (mrs_real) ((mrs_real)abs(error)/ outterWinRgt_);
00994           trackingScore_(h) += ((-1 * scoreFraction * localPeakAmp) * ((mrs_real)period / (mrs_real)maxPeriod_));
00995 
00996           //cout << "SCORE_FRACTION: " << ((-1 * scoreFraction * localPeakAmp) * ((mrs_real)period / (mrs_real)maxPeriod_)) << endl;
00997         }
00998 
00999         hypSignals_(h, beatPoint) = 1.0;
01000         beatCount_(h)++;
01001 
01002         //if(timeElapsed_ == 1228)
01003         //cout << "H:" << h << "-it:" << it << "; Error: " << error << "; Max: " << localPeakAmp
01004         //  << "; S: " << trackingScore_(h) << "; BeatPoint: " << beatPoint << "; period: " << period << endl;
01005         //it++;
01006         //don't evaluate beat if the outterWin around it surprasses the induction window
01007       } while(beatPoint < inSamples_);
01008     }
01009   }
01010 
01011   //MATLAB_PUT(beatHypotheses_, "BeatHypotheses");
01012 
01013   //Retrieve best M (nrPeriodHyps_) {period_, phase} pairs, by period:
01014   for(int i = 0; i < nrPeriodHyps_; i++)
01015   {
01016     //cout << "i- " << i << "; period: " << beatHypotheses_(i*nrPhasesPerPeriod_, 0) << endl;
01017     //if the given period_ > 0 (= valid period_)
01018     if(beatHypotheses_(i*nrPhasesPerPeriod_, 0) > 0)
01019     {
01020       for(int j = i*nrPhasesPerPeriod_; j < (i*nrPhasesPerPeriod_)+nrPhasesPerPeriod_; j++)
01021       {
01022         // if given phase > 0 (= valid phase)
01023         if(beatHypotheses_(j, 1) > 0)
01024         {
01025           if(trackingScore_(j) > maxLocalTrackingScore_(i))
01026           {
01027             maxLocalTrackingScore_(i) = trackingScore_(j);
01028             maxLocalTrackingScoreInd_(i) = j;
01029           }
01030         }
01031       }
01032 
01033       initPeriods_(i) = beatHypotheses_((mrs_natural)maxLocalTrackingScoreInd_(i), 0); //keep each init period_ hypothesis
01034       lastPeriods_(i) = lastLocalPeriod((mrs_natural)maxLocalTrackingScoreInd_(i));
01035       //initPhases_(i) = firstBeatPoint_((mrs_natural)maxLocalTrackingScoreInd_(i)) - (inSamples_-1 - inductionTime_);
01036       //lastPhases_(i) = lastBeatPoint((mrs_natural)maxLocalTrackingScoreInd_(i)) - (inSamples_-1 - inductionTime_) + lastPeriods_(i);
01037       initPhases_(i) = firstBeatPoint_((mrs_natural)maxLocalTrackingScoreInd_(i)) - (inSamples_-1 - timeElapsed_);
01038       lastPhases_(i) = lastBeatPoint((mrs_natural)maxLocalTrackingScoreInd_(i)) - (inSamples_-1 - timeElapsed_) + lastPeriods_(i);
01039       //for safechecking (really unlikely to happen) if lastPhases still don't surprasses the induction window then sum other period
01040       if(lastPhases_(i) <= timeElapsed_)
01041         lastPhases_(i) += lastPeriods_(i);
01042 
01043       //metricalSalience_(i) = beatHypotheses_((mrs_natural)maxLocalTrackingScoreInd_(i), 2) * 10000; //(x10000 for real value)
01044       //rawScore_(i) = pow(metricalSalience_(i),2) * maxLocalTrackingScore_(i) / sqrt(beatCount((mrs_natural)maxLocalTrackingScoreInd_(i)));
01045       rawScore_(i) =  maxLocalTrackingScore_(i);// * (initPeriods_(i) / (mrs_real) maxPeriod_);
01046       metricalSalience_(i) = beatHypotheses_((mrs_natural)maxLocalTrackingScoreInd_(i), 2);
01047       //rawScore_(i) = metricalSalience_(i) * (initPeriods_(i) / 105);
01048 
01049       //cout << i << "(" << (mrs_natural)maxLocalTrackingScoreInd_(i) << ")- corrScore:" << maxLocalTrackingScore_(i) << "; Phase: "
01050       //    << beatHypotheses_((mrs_natural)maxLocalTrackingScoreInd_(i), 1) << "; initPeriod: " << initPeriods_(i) << "; initPhase:"
01051       //    << initPhases_(i) << "; lastPeriod: " << lastPeriods_(i) << "; lastPhases: " << lastPhases_(i) << endl;
01052     }
01053   }
01054 
01055   mrs_real maxMetricalRelScore = NA; //to allow negative scores
01056   mrs_natural avgPeriod = 0;
01057   mrs_real avgLocalTrackingScore_ = 0.0;
01058   mrs_real maxGlobalTrackingScore_ = NA; //to allow negative scores
01059   mrs_natural maxMetricalRelScoreInd = -1;
01060   for(int i = 0; i < nrPeriodHyps_; i++)
01061   {
01062     //if the given period_ > 0 && given phase > 0 (= valid period_ && valid phase)
01063     if(initPeriods_(i) > MINIMUMREAL && initPhases_(i) > MINIMUMREAL)
01064     {
01065       metricalRelScore_(i) = calcRelationalScore(i, rawScore_);
01066 
01067       if(metricalRelScore_(i) > maxMetricalRelScore)
01068       {
01069         maxMetricalRelScore = metricalRelScore_(i);
01070         maxMetricalRelScoreInd = i;
01071       }
01072 
01073       if(maxLocalTrackingScore_(i) > maxGlobalTrackingScore_)
01074         maxGlobalTrackingScore_ = maxLocalTrackingScore_(i);
01075 
01076       avgLocalTrackingScore_ += maxLocalTrackingScore_(i);
01077       avgPeriod += (mrs_natural) initPeriods_(i);
01078 
01079       //cout << "MetrRelScore" << i << ":" << metricalRelScore_(i) << "; iPeriod:" << initPeriods_(i)
01080       //    << "; iPhase:" << initPhases_(i) << "(" << beatHypotheses_((mrs_natural)maxLocalTrackingScoreInd_(i), 1)
01081       //    << "); MaxMetrRelScore:" << maxMetricalRelScore << "(" << maxMetricalRelScoreInd << ")" << endl;
01082     }
01083   }
01084 
01085   //if requested output induction best period hypothesis (as estimated by ACF and after adjustment in induction)
01086   if(strcmp(inductionOut_.c_str(), "-1") != 0 && !dumbInduction_)
01087   {
01088     ostringstream ossInitPeriod, ossLastPeriod;
01089     fstream outStream;
01090     mrs_natural indInitPeriod = (mrs_natural) (((60.0 / initPeriods_(maxMetricalRelScoreInd))
01091                                 * (srcFs_ / hopSize_)) + 0.5); // (+ 0.5 for round integer)
01092     mrs_natural indLastPeriod = (mrs_natural) (((60.0 / lastPeriods_(maxMetricalRelScoreInd))
01093                                 * (srcFs_ / hopSize_)) + 0.5); // (+ 0.5 for round integer)
01094 
01095     //write best initPeriod calculation (from ACF peaks)
01096     ossInitPeriod << inductionOut_ << "_indInitTempo.txt";
01097     outStream.open(ossInitPeriod.str().c_str(), ios::out|ios::trunc);
01098     outStream << indInitPeriod;
01099     outStream.close();
01100 
01101     //write best lastPeriod calculation (after adjusting the initial best period till the end of induction)
01102     ossLastPeriod << inductionOut_ << "_indLastTempo.txt";
01103     outStream.open(ossLastPeriod.str().c_str(), ios::out|ios::trunc);
01104     outStream << indLastPeriod;
01105     outStream.close();
01106   }
01107 
01108   avgLocalTrackingScore_ /= nrPeriodHyps_;
01109   avgPeriod /= nrPeriodHyps_;
01110 
01111   for(int i = 0; i < nrPeriodHyps_; i++)
01112   {
01113     if(backtrace_)
01114     {
01115       //Period:
01116       out(i, 0) = initPeriods_(i);
01117       //Phase:
01118       out(i, 1) = initPhases_(i);
01119     }
01120     else
01121     {
01122       //Period:
01123       out(i, 0) = lastPeriods_(i);
01124       //Phase:
01125       out(i, 1) = lastPhases_(i);
01126     }
01127 
01128     //After-Induction Score Calulcation: ===========================================================
01129     mrs_real finalScore; //final score
01130     if(!dumbInduction_) //if not in dumbInduction consider metrical relations
01131     {
01132       mrs_real metricalRelFraction = 0.0;
01133       metricalRelFraction = (metricalRelScore_(i) / maxMetricalRelScore);
01134 
01135       //to avoid negative score inversions
01136       if(metricalRelScore_(i) < 0 && maxMetricalRelScore > 0 && maxGlobalTrackingScore_ > 0
01137           && abs(metricalRelScore_(i)) > maxMetricalRelScore)
01138         metricalRelFraction = -1;
01139       if(metricalRelScore_(i) < 0 && maxMetricalRelScore < 0 && maxGlobalTrackingScore_ > 0)
01140         metricalRelFraction = (maxMetricalRelScore / metricalRelScore_(i));
01141 
01142       //initially consider a proportion of the maximum tracking score for this induction window
01143       if(curBestScore_ == NA)
01144         finalScore = metricalRelFraction * maxGlobalTrackingScore_;
01145       else //at each trigger consider a proportion of the current best score
01146       {
01147         if(curBestScore_ > 0) //to avoid best score inversion
01148           finalScore = (metricalRelFraction * curBestScore_) * triggerBestScoreFactor_;
01149         else
01150           finalScore = (metricalRelFraction * curBestScore_) / triggerBestScoreFactor_;
01151 
01152         //cout << "relFrac: " << metricalRelFraction << "; curBestScore: " << curBestScore_ << "; finalScore: " << finalScore
01153         //  << "; triggerBestScoreFactor: " << triggerBestScoreFactor_ << endl;
01154       }
01155 
01156       //cout << i << "-> " << finalScore << "; MetrScore: " << metricalRelScore_(i) << "; MaxMetrRelScore: "
01157       //    << maxMetricalRelScore << "; MaxCorrScore: " << maxGlobalTrackingScore_ << "; initPeriod: " << initPeriods_(i)
01158       //    << "; initPhase: " << initPhases_(i) << "; lastPeriod: " << lastPeriods_(i) << "; lastPhase: "
01159       //    << lastPhases_(i) << endl;
01160 
01161       //cout << "scoreFunc:" << scoreFunc_ << "; avgMaxSum:" << avgMaxSum << "; avgPeriod:" << avgPeriod << endl;
01162     }
01163     else //if dumbInduction don't consider metrical relations
01164     {
01165       //cout << "DUMB INDUCTION" << endl;
01166 
01167       if(curBestScore_ == NA) //initially consider a proportion of the maximum tracking score for this induction window
01168         finalScore = rawScore_(i);
01169       else //at each trigger consider a proportion of the current best score
01170       {
01171         //final score = a fraction of the curBestScore given by the proportion of this agent's induction score
01172         //to the maximum induction score, and by a user-defined triggerBestScoreFactor
01173 
01174         //to avoid rawScoreFraction negative inversions
01175         mrs_real rawScoreFraction;
01176         //(to avoid divide by 0 if hypothesis exists)
01177         if(initPeriods_(i) > MINIMUMREAL && initPhases_(i) > MINIMUMREAL)
01178         {
01179           if(rawScore_(i) == 0.0) rawScore_(i) = MINIMUMREAL;
01180           if(maxGlobalTrackingScore_ == 0.0) maxGlobalTrackingScore_ = MINIMUMREAL;
01181         }
01182 
01183         //if rawScore negative and maxRawScore positive && rawScore or maxRawScore < 1
01184         if((rawScore_(i) < 0 && maxGlobalTrackingScore_ > 0) && (fabs(rawScore_(i)) < 1 || maxGlobalTrackingScore_ < 1))
01185           rawScoreFraction = rawScore_(i) * maxGlobalTrackingScore_;
01186         else
01187         {
01188           if(fabs(rawScore_(i)) < maxGlobalTrackingScore_) //if abs(rawScore) < maxRawScore
01189             rawScoreFraction = rawScore_(i) / maxGlobalTrackingScore_;
01190           else //if abs(rawScore) > maxRawScore
01191             rawScoreFraction = maxGlobalTrackingScore_ / rawScore_(i);
01192         }
01193 
01194         if(curBestScore_ > 0) //to avoid best score inversion
01195           finalScore = (rawScoreFraction * curBestScore_) * triggerBestScoreFactor_;
01196         else
01197           finalScore = (rawScoreFraction * curBestScore_) / triggerBestScoreFactor_;
01198 
01199         //cout << "rawSc: " << rawScore_(i) << "; maxRawSc: " << maxGlobalTrackingScore_ << "; rawFr: " << rawScoreFraction
01200         //  << "; curBestScore: " << curBestScore_ << "; finalScore: " << finalScore << "; triggerScFact: "
01201         //  << triggerBestScoreFactor_ << endl;
01202       }
01203 
01204       //reset dumb induction
01205       ctrl_dumbInduction_->setValue(false);
01206     }
01207 
01208     //aditional initial score normalization dependent of score function (due to their relative weights)
01209     if(strcmp(scoreFunc_.c_str(), "correlation") == 0 || strcmp(scoreFunc_.c_str(), "squareCorr") == 0 )
01210       finalScore *= 5; //("correlation" and "squareCorr" are, in average, 5times more reactive than "regular")
01211 
01212     out(i, 2) = finalScore;
01213 
01214     //cout << "i- " << i << "FINALSCORE: " << finalScore << endl;
01215   }
01216 
01217   //MATLAB_PUT(in, "Flux_FlowThrued2");
01218   //MATLAB_PUT(hypSignals_, "HypSignals");
01219   //MATLAB_PUT(beatHypotheses_, "BeatHypotheses");
01220   //MATLAB_PUT(out, "FinalHypotheses");
01221 }
01222 
01223 void
01224 PhaseLock::forceInitPeriods(mrs_string mode)
01225 {
01226   (void) mode; // seems to be done with mode_ instead?
01227   cerr << "\nInitial period(s) given by ground-truth file at: " << ctrl_gtBeatsFile_->to<mrs_string>() << endl;
01228 
01229   beatHypotheses_ = ctrl_beatHypotheses_->to<mrs_realvec>();
01230   if(strcmp(mode_.c_str(), "p") == 0)
01231   {
01232     //just pass gt period
01233     if(backtrace_) //if in non-causal or backtrace mode
01234     {
01235       for(int h = 0; h < nInitHyp_; h++) //give init gt period
01236         beatHypotheses_(h, 0) = (mrs_real) gtInitPeriod_;
01237 
01238       cerr << "Period as ibi of given first 2 beats: ";
01239       cerr << ((60.0/gtInitPeriod_)*(srcFs_ / hopSize_)) << " (BPMs)" << endl;
01240     }
01241     else
01242     {
01243       for(int h = 0; h < nInitHyp_; h++) //give last gt period
01244         beatHypotheses_(h, 0) = (mrs_real) gtLastPeriod_;
01245 
01246       cerr << "Period as ibi of given last 2 beats: ";
01247       cerr << ((60.0/gtLastPeriod_)*(srcFs_ / hopSize_)) << " (BPMs)" << endl;
01248     }
01249 
01250     nrPeriodHyps_ = 1;
01251   }
01252   else if(strcmp(mode_.c_str(), "p_mr") == 0)
01253   {
01254     nrPeriodHyps_ = 5;
01255     nInitHyp_ = nrPeriodHyps_ * nrPhasesPerPeriod_; //Nr. of initial hypothesis
01256     mrs_natural p = 0, ph = 0;
01257 
01258     mrs_real gtPeriod;
01259     if(backtrace_) //noncausal or backtrace mode
01260     {
01261       gtPeriod = gtInitPeriod_;
01262       cerr << "Periods as ibi of given first 2 beats + others metrical related: ";
01263     }
01264     else //causal mode
01265     {
01266       gtPeriod = gtLastPeriod_;
01267       cerr << "Periods as ibi of given last 2 beats + others metrical related: ";
01268     }
01269 
01270     double periods[] = {
01271       (double) gtPeriod,
01272       (double) (mrs_natural) (gtPeriod * 2.0),
01273       (double) (mrs_natural) (gtPeriod * 0.5),
01274       (double) (mrs_natural) (gtPeriod * 3.0),
01275       (double) (mrs_natural) (gtPeriod * 0.333)
01276     };
01277     //assume gt period + (4) others metrical related (2x, 1/2x, 3x, 1/3x)
01278     for(int i = 0; i < nrPeriodHyps_; i++)
01279     {
01280       for(int h = i*nrPhasesPerPeriod_; h < (i*nrPhasesPerPeriod_)+nrPhasesPerPeriod_; h++)
01281       {
01282         beatHypotheses_(h, 0) = periods[p];
01283         beatHypotheses_(h, 1) = beatHypotheses_(ph, 1);
01284         //cout << "i-" << i <<"; h-" << h << "; p-" << p << ": " << ((60.0/periods[p])*(srcFs_ / hopSize_)) << "; ph: " << beatHypotheses_(h, 1) << endl;
01285         ph++;
01286       }
01287       p++;
01288       ph = 0;
01289     }
01290 
01291     //redefine tempo ranges to account for imposed values
01292     if(periods[3] > maxPeriod_) maxPeriod_ = ((mrs_natural) (periods[3] + 0.5));
01293     if(periods[4] < minPeriod_) minPeriod_ = ((mrs_natural) (periods[4] + 0.5));
01294     updControl(ctrl_maxPeriod_, maxPeriod_);
01295     updControl(ctrl_minPeriod_, minPeriod_);
01296 
01297     cerr << ((60.0/periods[0])*(srcFs_ / hopSize_)) << "; " << ((60.0/periods[1])*(srcFs_ / hopSize_)) << "; " <<
01298          ((60.0/periods[2])*(srcFs_ / hopSize_)) << "; " << ((60.0/periods[3])*(srcFs_ / hopSize_)) << "; " <<
01299          ((60.0/periods[4])*(srcFs_ / hopSize_)) << " (BPMs)" << endl;
01300   }
01301   else if(strcmp(mode_.c_str(), "p_nr") == 0)
01302   {
01303     nrPeriodHyps_ = 5;
01304     nInitHyp_ = nrPeriodHyps_ * nrPhasesPerPeriod_; //Nr. of initial hypothesis
01305     mrs_natural p = 0, ph = 0;
01306 
01307     mrs_real gtPeriod;
01308     if(backtrace_) //noncausal or backtrace mode
01309     {
01310       gtPeriod = gtInitPeriod_;
01311       cerr << "Periods as ibi of given first 2 beats + others non-related: ";
01312     }
01313     else //causal mode
01314     {
01315       gtPeriod = gtLastPeriod_;
01316       cerr << "Periods as ibi of given last 2 beats + others non-related: ";
01317     }
01318 
01319     double periods[] = {
01320       (double) gtPeriod,
01321       (double) (mrs_natural) (gtInitPeriod_ * 1.8),
01322       (double) (mrs_natural) (gtPeriod * 1.2),
01323       (double) (mrs_natural) (gtPeriod * 2.3),
01324       (double) (mrs_natural) (gtPeriod * 0.7)
01325     };
01326     //assume gt period + (4) others metrical related (2x, 1/2x, 3x, 1/3x)
01327     for(int i = 0; i < nrPeriodHyps_; i++)
01328     {
01329       for(int h = i*nrPhasesPerPeriod_; h < (i*nrPhasesPerPeriod_)+nrPhasesPerPeriod_; h++)
01330       {
01331         beatHypotheses_(h, 0) = periods[p];
01332         beatHypotheses_(h, 1) = beatHypotheses_(ph, 1);
01333         //cout << "i-" << i <<"; h-" << h << "; p-" << p << ": " << ((60.0/periods[p])*(srcFs_ / hopSize_)) << "; ph: " << beatHypotheses_(h, 1) << endl;
01334         ph++;
01335       }
01336       p++;
01337       ph = 0;
01338     }
01339 
01340     //redefine tempo ranges to account for imposed values
01341     if(periods[3] > maxPeriod_) maxPeriod_ = ((mrs_natural) (periods[3] + 0.5));
01342     if(periods[4] < minPeriod_) minPeriod_ = ((mrs_natural) (periods[4] + 0.5));
01343     updControl(ctrl_maxPeriod_, maxPeriod_);
01344     updControl(ctrl_minPeriod_, minPeriod_);
01345 
01346     cerr << ((60.0/periods[0])*(srcFs_ / hopSize_)) << "; " << ((60.0/periods[1])*(srcFs_ / hopSize_)) << "; " <<
01347          ((60.0/periods[2])*(srcFs_ / hopSize_)) << "; " << ((60.0/periods[3])*(srcFs_ / hopSize_)) << "; " <<
01348          ((60.0/periods[4])*(srcFs_ / hopSize_)) << " (BPMs)" << endl;
01349   }
01350 
01351   updControl(ctrl_beatHypotheses_, beatHypotheses_);
01352 }
01353 
01354 void
01355 PhaseLock::myProcess(realvec& in, realvec& out)
01356 {
01357   mrs_natural o,t;
01358   //timeElapsed_ is constantly updated with the referee's next time frame
01359   timeElapsed_ = ctrl_tickCount_->to<mrs_natural>();
01360 
01361   //cout << "PLock: " << timeElapsed_ << "; bestScoreFact: " << triggerBestScoreFactor_ << endl;
01362 
01363   //cout << "PLock: " << timeElapsed_ << "; Ind: " << inductionTime_ << "; maxPer: " << maxPeriod_ << "; minPer: " << minPeriod_ << endl;
01364 
01365   //Output only defined just after induction time
01366   //until then output is undefined...
01367   for (o=0; o < onObservations_; o++)
01368   {
01369     for (t = 0; t < onSamples_; t++)
01370     {
01371       out(o,t) = -1.0;
01372     }
01373   }
01374 
01375   //After the given induction stage, generate and return the correspondent beat hypotheses
01376   //[occuring in the beginning of the analysis or if triggered in "trigger induction" mode]
01377   //Calculate the N (nrPeriodHyps_) best (period_, phase) pairs + respective scores:
01378   triggerInduction_ = ctrl_triggerInduction_->to<mrs_bool>();
01379   if(triggerInduction_)
01380   {
01381     dumbInduction_ = ctrl_dumbInduction_->to<mrs_bool>();
01382 
01383     //cout << "TRIGGER @ " << timeElapsed_ << endl;
01384 
01385     cerr << "\nRequested Induction in \"" << mode_ << "\" mode at: " << (((timeElapsed_* hopSize_)-hopSize_/2) / srcFs_) << "s" << endl; //(" << timeElapsed_ << ")" << endl;
01386 
01387     //maxScore_ = calcGTNormScore(in);
01388     if(strcmp(mode_.c_str(), "2b") == 0 || strcmp(mode_.c_str(), "2b2") == 0) //gt period+phase (adjusted)
01389     {
01390       mrs_realvec gtHypotheses = readGTFile(ctrl_gtBeatsFile_->to<mrs_string>()); //process gt file
01391       handleGTHypotheses(in, out, ctrl_gtBeatsFile_->to<mrs_string>(), gtHypotheses);
01392 
01393       if(backtrace_)
01394       {
01395         //Period:
01396         out(0, 0) = gtInitPeriod_; //keep period as ibi from first 2 beats in gt file
01397         //Phase:
01398         out(0, 1) = gtInitPhase_; //first/second beat in gt file
01399 
01400         if(strcmp(mode_.c_str(), "2b") == 0)
01401         {
01402           if(gtAfter2ndBeat_)
01403             cerr << "Initial phase backtraced from second beat of given ground-truth file: ";
01404           else
01405             cerr << "Initial phase backtraced from first beat of given ground-truth file: ";
01406         }
01407         else if(strcmp(mode_.c_str(), "2b2") == 0)
01408         {
01409           if(gtAfter2ndBeat_)
01410             cerr << "Initial phase as second beat of given ground-truth file: ";
01411           else
01412             cerr << "Initial phase as first beat of given ground-truth file: ";
01413         }
01414 
01415         cerr << (((gtInitPhase_ * hopSize_) - adjustment_) / srcFs_) << "s" << endl;
01416         cerr << "Ground-truth period: " << ((60.0/gtInitPeriod_)*(srcFs_ / hopSize_)) << " (BPMs)" << endl;
01417       }
01418       else //no backtrace (causal mode)
01419       {
01420         //Period:
01421         out(0, 0) = gtLastPeriod_; //period adjusted from ibi given by first 2 beats in gt file
01422         //Phase:
01423         out(0, 1) = gtLastPhase_; //adjusted from first/second beat in gt file
01424 
01425         if(strcmp(mode_.c_str(), "2b") == 0)
01426         {
01427           if(gtAfter2ndBeat_)
01428             cerr << "Initial phase adjusted from second beat after induction time, of given ground-truth file: ";
01429           else
01430             cerr << "Initial phase adjusted from first beat after induction time, of given ground-truth file: ";
01431 
01432           cerr << (((gtLastPhase_ * hopSize_) - adjustment_) / srcFs_) << "s" << endl;
01433           cerr << "Ground-truth period (adjusted): " << ((60.0/gtLastPeriod_)*(srcFs_ / hopSize_)) << " (BPMs)" << endl;
01434         }
01435         else if(strcmp(mode_.c_str(), "2b2") == 0)
01436         {
01437           cerr << "Initial phase as first beat after induction time, of given ground-truth file: ";
01438           cerr << (((gtLastPhase_ * hopSize_) - adjustment_) / srcFs_) << "s" << "[" << gtLastPhase_ << "]" << endl;
01439           cerr << "Ground-truth period: " << ((60.0/gtLastPeriod_)*(srcFs_ / hopSize_)) << " (BPMs)" << endl;
01440         }
01441         //cout << "CALC Initial phase: " << gtInitPhase_ << "(" << (((gtInitPhase_ * hopSize_) - adjustment_) / srcFs_) << ")" << endl;
01442       }
01443 
01444       //score
01445       out(0, 2) = gtScore_; //initialScore
01446 
01447       //cout << "gtLPe: " << gtLastPeriod_ << "(" << out(0,0) << ")" << "gtLPh: " << gtLastPhase_ << "(" << out(0,1) << ")" <<
01448       //    "gtSc: " << gtScore_ << "(" << out(0,3) << ")" << endl;
01449     }
01450     else if(strcmp(mode_.c_str(), "1b") == 0
01451             || strcmp(mode_.c_str(), "1b1") == 0 ) //gt phase
01452     {
01453       mrs_realvec gtHypotheses = readGTFile(ctrl_gtBeatsFile_->to<mrs_string>()); //process gt file
01454       handleGTHypotheses(in, out, ctrl_gtBeatsFile_->to<mrs_string>(), gtHypotheses);
01455     }
01456     else if((strcmp(mode_.c_str(), "p") == 0) || (strcmp(mode_.c_str(), "p_mr") == 0)
01457             || (strcmp(mode_.c_str(), "p_nr") == 0))
01458     {
01459       mrs_realvec gtHypotheses = readGTFile(ctrl_gtBeatsFile_->to<mrs_string>()); //process gt file
01460       handleGTHypotheses(in, out, ctrl_gtBeatsFile_->to<mrs_string>(), gtHypotheses);
01461 
01462       //force gt period (given by gt file - ibi of first 2 beats) in beatHypotheses vector
01463       forceInitPeriods(mode_);
01464 
01465       //run as regular induction
01466       regularFunc(in, out);
01467     }
01468     else if(strcmp(mode_.c_str(), "regular") == 0)
01469       regularFunc(in , out);
01470 
01471     cerr << "===================FINISH INDUCTION=====================" << endl;
01472 
01473     /*
01474     ostringstream oss;
01475     fstream outStream;
01476     oss << "phaselock_test.txt";
01477     for(int i = 0; i < in.getCols(); i++)
01478     {
01479         if(i == 0)
01480         {
01481             outStream.open(oss.str().c_str(), ios::out|ios::trunc);
01482             outStream << in(i) << endl;
01483             outStream.close();
01484         }
01485         else
01486         {
01487             outStream.open(oss.str().c_str(), ios::out|ios::app);
01488             outStream << in(i) << endl;
01489             outStream.close();
01490         }
01491     }
01492     */
01493   }
01494 
01495   //MATLAB_PUT(in, "Flux_FlowThrued");
01496   //MATLAB_PUT(out, "PhaseLockOut");
01497   //MATLAB_EVAL("hold on;");
01498   //MATLAB_EVAL("plot(Flux_FlowThrued), g");
01499   //MATLAB_EVAL("FluxFlowTS = [FluxFlowTS, Flux_FlowThrued];");
01500 }