c@277: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@277: 
c@277: /*
c@277:     QM DSP Library
c@277: 
c@277:     Centre for Digital Music, Queen Mary, University of London.
c@277:     This file copyright 2008-2009 Matthew Davies and QMUL.
c@309: 
c@309:     This program is free software; you can redistribute it and/or
c@309:     modify it under the terms of the GNU General Public License as
c@309:     published by the Free Software Foundation; either version 2 of the
c@309:     License, or (at your option) any later version.  See the file
c@309:     COPYING included with this distribution for more information.
c@277: */
c@277: 
c@277: #include "TempoTrackV2.h"
c@277: 
c@277: #include <cmath>
c@277: #include <cstdlib>
c@278: #include <iostream>
c@277: 
c@279: #include "maths/MathUtilities.h"
c@277: 
cannam@493: using std::vector;
cannam@493: 
c@277: #define   EPS 0.0000008 // just some arbitrary small number
c@277: 
cannam@493: TempoTrackV2::TempoTrackV2(float rate, int increment) :
cannam@493:     m_rate(rate), m_increment(increment) {
cannam@493: }
cannam@479: 
c@277: TempoTrackV2::~TempoTrackV2() { }
c@277: 
c@277: void
c@277: TempoTrackV2::filter_df(d_vec_t &df)
c@277: {
cannam@501:     int df_len = int(df.size());
cannam@501:     
c@278:     d_vec_t a(3);
c@278:     d_vec_t b(3);
cannam@501:     d_vec_t lp_df(df_len);
c@277: 
c@278:     //equivalent in matlab to [b,a] = butter(2,0.4);
c@278:     a[0] = 1.0000;
c@278:     a[1] = -0.3695;
c@278:     a[2] = 0.1958;
c@278:     b[0] = 0.2066;
c@278:     b[1] = 0.4131;
c@278:     b[2] = 0.2066;
luis@327: 
c@278:     double inp1 = 0.;
c@278:     double inp2 = 0.;
c@278:     double out1 = 0.;
c@278:     double out2 = 0.;
c@277: 
c@277: 
c@278:     // forwards filtering
cannam@501:     for (int i = 0; i < df_len; i++) {
c@278:         lp_df[i] =  b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2;
c@278:         inp2 = inp1;
c@278:         inp1 = df[i];
c@278:         out2 = out1;
c@278:         out1 = lp_df[i];
c@278:     }
c@277: 
c@278:     // copy forwards filtering to df...
c@278:     // but, time-reversed, ready for backwards filtering
cannam@501:     for (int i = 0; i < df_len; i++) {
cannam@501:         df[i] = lp_df[df_len - i - 1];
c@278:     }
c@277: 
cannam@501:     for (int i = 0; i < df_len; i++) {
luis@327:         lp_df[i] = 0.;
c@278:     }
c@277: 
c@278:     inp1 = 0.; inp2 = 0.;
c@278:     out1 = 0.; out2 = 0.;
c@277: 
cannam@479:     // backwards filetering on time-reversed df
cannam@501:     for (int i = 0; i < df_len; i++) {
c@278:         lp_df[i] =  b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2;
c@278:         inp2 = inp1;
c@278:         inp1 = df[i];
c@278:         out2 = out1;
c@278:         out1 = lp_df[i];
c@278:     }
c@277: 
cannam@479:     // write the re-reversed (i.e. forward) version back to df
cannam@501:     for (int i = 0; i < df_len; i++) {
cannam@501:         df[i] = lp_df[df_len - i - 1];
c@278:     }
c@277: }
c@277: 
c@277: 
luis@327: // MEPD 28/11/12
luis@327: // This function now allows for a user to specify an inputtempo (in BPM)
luis@327: // and a flag "constraintempo" which replaces the general rayleigh weighting for periodicities
luis@327: // with a gaussian which is centered around the input tempo
luis@327: // Note, if inputtempo = 120 and constraintempo = false, then functionality is
luis@327: // as it was before
c@277: void
c@304: TempoTrackV2::calculateBeatPeriod(const vector<double> &df,
c@304:                                   vector<double> &beat_period,
luis@327:                                   vector<double> &tempi,
luis@327:                                   double inputtempo, bool constraintempo)
c@277: {
c@278:     // to follow matlab.. split into 512 sample frames with a 128 hop size
c@278:     // calculate the acf,
c@278:     // then the rcf.. and then stick the rcfs as columns of a matrix
c@278:     // then call viterbi decoding with weight vector and transition matrix
c@278:     // and get best path
c@277: 
cannam@501:     int wv_len = 128;
luis@327: 
luis@327:     // MEPD 28/11/12
luis@327:     // the default value of inputtempo in the beat tracking plugin is 120
luis@327:     // so if the user specifies a different inputtempo, the rayparam will be updated
luis@327:     // accordingly.
luis@327:     // note: 60*44100/512 is a magic number
luis@327:     // this might (will?) break if a user specifies a different frame rate for the onset detection function
luis@327:     double rayparam = (60*44100/512)/inputtempo;
luis@327: 
c@278:     // make rayleigh weighting curve
c@278:     d_vec_t wv(wv_len);
luis@327: 
luis@327:     // check whether or not to use rayleigh weighting (if constraintempo is false)
luis@327:     // or use gaussian weighting it (constraintempo is true)
cannam@479:     if (constraintempo) {
cannam@501:         for (int i = 0; i < wv_len; i++) {
luis@327:             // MEPD 28/11/12
luis@327:             // do a gaussian weighting instead of rayleigh
cannam@501:             wv[i] = exp( (-1.*pow((double(i)-rayparam),2.)) / (2.*pow(rayparam/4.,2.)) );
luis@327:         }
cannam@479:     } else {
cannam@501:         for (int i = 0; i < wv_len; i++) {
luis@327:             // MEPD 28/11/12
luis@327:             // standard rayleigh weighting over periodicities
cannam@501:             wv[i] = (double(i) / pow(rayparam,2.)) * exp((-1.*pow(-double(i),2.)) / (2.*pow(rayparam,2.)));
luis@327:         }
c@277:     }
c@277: 
c@278:     // beat tracking frame size (roughly 6 seconds) and hop (1.5 seconds)
cannam@501:     int winlen = 512;
cannam@501:     int step = 128;
c@278: 
c@278:     // matrix to store output of comb filter bank, increment column of matrix at each frame
c@278:     d_mat_t rcfmat;
c@278:     int col_counter = -1;
cannam@501:     int df_len = int(df.size());
c@278: 
c@278:     // main loop for beat period calculation
cannam@501:     for (int i = 0; i+winlen < df_len; i+=step) {
cannam@479:         
c@278:         // get dfframe
c@278:         d_vec_t dfframe(winlen);
cannam@501:         for (int k=0; k < winlen; k++) {
c@278:             dfframe[k] = df[i+k];
c@278:         }
c@278:         // get rcf vector for current frame
luis@327:         d_vec_t rcf(wv_len);
c@278:         get_rcf(dfframe,wv,rcf);
luis@327: 
c@278:         rcfmat.push_back( d_vec_t() ); // adds a new column
c@278:         col_counter++;
cannam@501:         for (int j = 0; j < wv_len; j++) {
c@278:             rcfmat[col_counter].push_back( rcf[j] );
c@278:         }
c@278:     }
luis@327: 
c@278:     // now call viterbi decoding function
c@278:     viterbi_decode(rcfmat,wv,beat_period,tempi);
c@277: }
c@277: 
c@277: 
c@277: void
c@277: TempoTrackV2::get_rcf(const d_vec_t &dfframe_in, const d_vec_t &wv, d_vec_t &rcf)
c@277: {
c@278:     // calculate autocorrelation function
c@278:     // then rcf
c@278:     // just hard code for now... don't really need separate functions to do this
c@277: 
c@278:     // make acf
c@277: 
c@278:     d_vec_t dfframe(dfframe_in);
c@277: 
c@279:     MathUtilities::adaptiveThreshold(dfframe);
c@277: 
cannam@501:     int dfframe_len = int(dfframe.size());
cannam@501:     int rcf_len = int(rcf.size());
cannam@501:     
cannam@501:     d_vec_t acf(dfframe_len);
c@277: 
cannam@501:     for (int lag = 0; lag < dfframe_len; lag++) {
c@278:         double sum = 0.;
c@278:         double tmp = 0.;
c@277: 
cannam@501:         for (int n = 0; n < (dfframe_len - lag); n++) {
cannam@501:             tmp = dfframe[n] * dfframe[n + lag];
c@278:             sum += tmp;
c@278:         }
cannam@501:         acf[lag] = double(sum/ (dfframe_len - lag));
c@278:     }
c@277: 
c@278:     // now apply comb filtering
c@278:     int numelem = 4;
luis@327: 
cannam@501:     for (int i = 2; i < rcf_len; i++) { // max beat period
cannam@501:         for (int a = 1; a <= numelem; a++) { // number of comb elements
cannam@501:             for (int b = 1-a; b <= a-1; b++) { // general state using normalisation of comb elements
cannam@479:                 rcf[i-1] += ( acf[(a*i+b)-1]*wv[i-1] ) / (2.*a-1.);     // calculate value for comb filter row
c@278:             }
c@278:         }
c@278:     }
luis@327: 
c@278:     // apply adaptive threshold to rcf
c@279:     MathUtilities::adaptiveThreshold(rcf);
luis@327: 
c@278:     double rcfsum =0.;
cannam@501:     for (int i = 0; i < rcf_len; i++) {
c@278:         rcf[i] += EPS ;
c@278:         rcfsum += rcf[i];
c@278:     }
c@277: 
c@278:     // normalise rcf to sum to unity
cannam@501:     for (int i = 0; i < rcf_len; i++) {
c@278:         rcf[i] /= (rcfsum + EPS);
c@277:     }
c@277: }
c@277: 
c@277: void
c@278: TempoTrackV2::viterbi_decode(const d_mat_t &rcfmat, const d_vec_t &wv, d_vec_t &beat_period, d_vec_t &tempi)
c@277: {
c@278:     // following Kevin Murphy's Viterbi decoding to get best path of
c@278:     // beat periods through rfcmat
cannam@501:     
cannam@501:     int wv_len = int(wv.size());
cannam@501:     
c@278:     // make transition matrix
c@278:     d_mat_t tmat;
cannam@501:     for (int i = 0; i < wv_len; i++) {
c@278:         tmat.push_back ( d_vec_t() ); // adds a new column
cannam@501:         for (int j = 0; j < wv_len; j++) {
c@278:             tmat[i].push_back(0.); // fill with zeros initially
c@278:         }
c@278:     }
luis@327: 
c@278:     // variance of Gaussians in transition matrix
c@278:     // formed of Gaussians on diagonal - implies slow tempo change
c@278:     double sigma = 8.;
c@278:     // don't want really short beat periods, or really long ones
cannam@501:     for (int i = 20; i  < wv_len - 20; i++) {
cannam@501:         for (int j = 20; j < wv_len - 20; j++) {
cannam@501:             double mu = double(i);
c@278:             tmat[i][j] = exp( (-1.*pow((j-mu),2.)) / (2.*pow(sigma,2.)) );
c@278:         }
c@278:     }
c@277: 
c@278:     // parameters for Viterbi decoding... this part is taken from
c@278:     // Murphy's matlab
c@277: 
c@278:     d_mat_t delta;
c@278:     i_mat_t psi;
cannam@501:     for (int i = 0; i < int(rcfmat.size()); i++) {
cannam@501:         delta.push_back(d_vec_t());
cannam@501:         psi.push_back(i_vec_t());
cannam@501:         for (int j = 0; j < int(rcfmat[i].size()); j++) {
c@278:             delta[i].push_back(0.); // fill with zeros initially
c@278:             psi[i].push_back(0); // fill with zeros initially
c@278:         }
c@278:     }
c@277: 
cannam@502:     int T = int(delta.size());
c@281: 
c@281:     if (T < 2) return; // can't do anything at all meaningful
c@281: 
cannam@502:     int Q = int(delta[0].size());
c@277: 
c@278:     // initialize first column of delta
cannam@501:     for (int j = 0; j < Q; j++) {
c@278:         delta[0][j] = wv[j] * rcfmat[0][j];
c@278:         psi[0][j] = 0;
c@277:     }
luis@327: 
c@277:     double deltasum = 0.;
cannam@501:     for (int i = 0; i < Q; i++) {
c@278:         deltasum += delta[0][i];
luis@327:     }
cannam@501:     for (int i = 0; i < Q; i++) {
c@278:         delta[0][i] /= (deltasum + EPS);
luis@327:     }
c@277: 
cannam@501:     for (int t=1; t < T; t++)
c@278:     {
c@278:         d_vec_t tmp_vec(Q);
c@277: 
cannam@501:         for (int j = 0; j < Q; j++) {
cannam@501:             for (int i = 0; i < Q; i++) {
c@278:                 tmp_vec[i] = delta[t-1][i] * tmat[j][i];
luis@327:             }
luis@327: 
luis@327:             delta[t][j] = get_max_val(tmp_vec);
c@277: 
c@278:             psi[t][j] = get_max_ind(tmp_vec);
luis@327: 
c@278:             delta[t][j] *= rcfmat[t][j];
c@278:         }
c@277: 
c@278:         // normalise current delta column
c@278:         double deltasum = 0.;
cannam@501:         for (int i = 0; i < Q; i++) {
c@278:             deltasum += delta[t][i];
luis@327:         }
cannam@501:         for (int i = 0; i < Q; i++) {
c@278:             delta[t][i] /= (deltasum + EPS);
luis@327:         }
c@278:     }
c@277: 
c@278:     i_vec_t bestpath(T);
c@278:     d_vec_t tmp_vec(Q);
cannam@501:     for (int i = 0; i < Q; i++) {
c@278:         tmp_vec[i] = delta[T-1][i];
c@278:     }
c@277: 
c@278:     // find starting point - best beat period for "last" frame
c@278:     bestpath[T-1] = get_max_ind(tmp_vec);
luis@327: 
c@278:     // backtrace through index of maximum values in psi
cannam@501:     for (int t=T-2; t>0 ;t--) {
c@278:         bestpath[t] = psi[t+1][bestpath[t+1]];
c@278:     }
c@277: 
c@278:     // weird but necessary hack -- couldn't get above loop to terminate at t >= 0
c@278:     bestpath[0] = psi[1][bestpath[1]];
c@277: 
cannam@501:     int lastind = 0;
cannam@501:     for (int i = 0; i < T; i++) {
cannam@501:         int step = 128;
cannam@501:         for (int j = 0; j < step; j++) {
c@278:             lastind = i*step+j;
c@278:             beat_period[lastind] = bestpath[i];
c@278:         }
c@282: //        std::cerr << "bestpath[" << i << "] = " << bestpath[i] << " (used for beat_periods " << i*step << " to " << i*step+step-1 << ")" << std::endl;
c@278:     }
c@277: 
cannam@501:     // fill in the last values...
cannam@501:     for (int i = lastind; i < int(beat_period.size()); i++) {
c@278:         beat_period[i] = beat_period[lastind];
c@278:     }
c@277: 
cannam@501:     for (int i = 0; i < int(beat_period.size()); i++) {
c@279:         tempi.push_back((60. * m_rate / m_increment)/beat_period[i]);
c@277:     }
c@277: }
c@277: 
c@277: double
c@277: TempoTrackV2::get_max_val(const d_vec_t &df)
c@277: {
c@278:     double maxval = 0.;
cannam@501:     int df_len = int(df.size());
cannam@501:     
cannam@501:     for (int i = 0; i < df_len; i++) {
cannam@479:         if (maxval < df[i]) {
c@278:             maxval = df[i];
c@278:         }
c@277:     }
luis@327: 
c@278:     return maxval;
c@277: }
c@277: 
c@277: int
c@277: TempoTrackV2::get_max_ind(const d_vec_t &df)
c@277: {
c@278:     double maxval = 0.;
c@278:     int ind = 0;
cannam@501:     int df_len = int(df.size());
cannam@501:     
cannam@501:     for (int i = 0; i < df_len; i++) {
cannam@479:         if (maxval < df[i]) {
c@278:             maxval = df[i];
c@278:             ind = i;
c@278:         }
c@277:     }
luis@327: 
c@278:     return ind;
c@277: }
c@277: 
c@277: void
c@277: TempoTrackV2::normalise_vec(d_vec_t &df)
c@277: {
c@278:     double sum = 0.;
cannam@501:     int df_len = int(df.size());
cannam@501:     
cannam@501:     for (int i = 0; i < df_len; i++) {
c@278:         sum += df[i];
c@278:     }
luis@327: 
cannam@501:     for (int i = 0; i < df_len; i++) {
c@278:         df[i]/= (sum + EPS);
c@278:     }
c@277: }
c@277: 
luis@327: // MEPD 28/11/12
luis@327: // this function has been updated to allow the "alpha" and "tightness" parameters
luis@327: // of the dynamic program to be set by the user
luis@327: // the default value of alpha = 0.9 and tightness = 4
c@277: void
c@304: TempoTrackV2::calculateBeats(const vector<double> &df,
c@304:                              const vector<double> &beat_period,
luis@327:                              vector<double> &beats, double alpha, double tightness)
c@277: {
c@281:     if (df.empty() || beat_period.empty()) return;
c@281: 
cannam@501:     int df_len = int(df.size());
c@277: 
cannam@501:     d_vec_t cumscore(df_len); // store cumulative score
cannam@501:     i_vec_t backlink(df_len); // backlink (stores best beat locations at each time instant)
cannam@501:     d_vec_t localscore(df_len); // localscore, for now this is the same as the detection function
cannam@501: 
cannam@501:     for (int i = 0; i < df_len; i++) {
c@278:         localscore[i] = df[i];
c@278:         backlink[i] = -1;
c@277:     }
c@277: 
luis@327:     //double tightness = 4.;
luis@327:     //double alpha = 0.9;
luis@327:     // MEPD 28/11/12
luis@327:     // debug statements that can be removed.
c@330: //    std::cerr << "alpha" << alpha << std::endl;
c@330: //    std::cerr << "tightness" << tightness << std::endl;
c@277: 
c@278:     // main loop
cannam@501:     for (int i = 0; i < df_len; i++) {
cannam@479:         
c@278:         int prange_min = -2*beat_period[i];
c@278:         int prange_max = round(-0.5*beat_period[i]);
c@277: 
c@278:         // transition range
cannam@501:         int txwt_len = prange_max - prange_min + 1;
cannam@501:         d_vec_t txwt (txwt_len);
cannam@501:         d_vec_t scorecands (txwt_len);
c@277: 
cannam@501:         for (int j = 0; j < txwt_len; j++) {
cannam@479:             
cannam@501:             double mu = double(beat_period[i]);
c@278:             txwt[j] = exp( -0.5*pow(tightness * log((round(2*mu)-j)/mu),2));
c@277: 
c@278:             // IF IN THE ALLOWED RANGE, THEN LOOK AT CUMSCORE[I+PRANGE_MIN+J
c@278:             // ELSE LEAVE AT DEFAULT VALUE FROM INITIALISATION:  D_VEC_T SCORECANDS (TXWT.SIZE());
c@277: 
cannam@501:             int cscore_ind = i + prange_min + j;
cannam@479:             if (cscore_ind >= 0) {
c@278:                 scorecands[j] = txwt[j] * cumscore[cscore_ind];
c@278:             }
c@278:         }
c@277: 
c@278:         // find max value and index of maximum value
c@278:         double vv = get_max_val(scorecands);
c@278:         int xx = get_max_ind(scorecands);
c@277: 
c@278:         cumscore[i] = alpha*vv + (1.-alpha)*localscore[i];
c@278:         backlink[i] = i+prange_min+xx;
c@280: 
c@282: //        std::cerr << "backlink[" << i << "] <= " << backlink[i] << std::endl;
c@278:     }
c@278: 
c@278:     // STARTING POINT, I.E. LAST BEAT.. PICK A STRONG POINT IN cumscore VECTOR
c@278:     d_vec_t tmp_vec;
cannam@501:     for (int i = df_len - beat_period[beat_period.size()-1] ; i < df_len; i++) {
c@278:         tmp_vec.push_back(cumscore[i]);
luis@327:     }
c@278: 
cannam@479:     int startpoint = get_max_ind(tmp_vec) +
cannam@501:         df_len - beat_period[beat_period.size()-1] ;
c@278: 
c@281:     // can happen if no results obtained earlier (e.g. input too short)
cannam@501:     if (startpoint >= int(backlink.size())) {
cannam@501:         startpoint = int(backlink.size()) - 1;
cannam@479:     }
c@281: 
c@278:     // USE BACKLINK TO GET EACH NEW BEAT (TOWARDS THE BEGINNING OF THE FILE)
c@278:     //  BACKTRACKING FROM THE END TO THE BEGINNING.. MAKING SURE NOT TO GO BEFORE SAMPLE 0
c@278:     i_vec_t ibeats;
c@278:     ibeats.push_back(startpoint);
c@282: //    std::cerr << "startpoint = " << startpoint << std::endl;
cannam@479:     while (backlink[ibeats.back()] > 0) {
c@282: //        std::cerr << "backlink[" << ibeats.back() << "] = " << backlink[ibeats.back()] << std::endl;
c@281:         int b = ibeats.back();
c@281:         if (backlink[b] == b) break; // shouldn't happen... haha
c@281:         ibeats.push_back(backlink[b]);
c@278:     }
luis@327: 
c@278:     // REVERSE SEQUENCE OF IBEATS AND STORE AS BEATS
cannam@501:     for (int i = 0; i < int(ibeats.size()); i++) {
cannam@501:         beats.push_back(double(ibeats[ibeats.size() - i - 1]));
c@278:     }
c@277: }
c@277: 
c@277: