cannam@52: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@52: 
cannam@52: /*
cannam@52:     QM DSP Library
cannam@52: 
cannam@52:     Centre for Digital Music, Queen Mary, University of London.
cannam@52:     This file copyright 2008-2009 Matthew Davies and QMUL.
cannam@52:     All rights reserved.
cannam@52: */
cannam@52: 
cannam@52: #include "TempoTrackV2.h"
cannam@52: 
cannam@52: #include <cmath>
cannam@52: #include <cstdlib>
cannam@52: 
cannam@52: 
cannam@52: //#define		FRAMESIZE	512
cannam@52: //#define		BIGFRAMESIZE	1024
cannam@52: #define   TWOPI 6.283185307179586232
cannam@52: #define   EPS 0.0000008 // just some arbitrary small number
cannam@52: 
cannam@52: TempoTrackV2::TempoTrackV2() { }
cannam@52: TempoTrackV2::~TempoTrackV2() { }
cannam@52: 
cannam@52: void
cannam@52: TempoTrackV2::adapt_thresh(d_vec_t &df)
cannam@52: {
cannam@52: 
cannam@52:   d_vec_t smoothed(df.size());
cannam@52: 	
cannam@52: 	int p_post = 7;
cannam@52: 	int p_pre = 8;
cannam@52: 
cannam@52: 	int t = std::min(static_cast<int>(df.size()),p_post);	// what is smaller, p_post of df size. This is to avoid accessing outside of arrays
cannam@52: 
cannam@52: 	// find threshold for first 't' samples, where a full average cannot be computed yet 
cannam@52: 	for (int i = 0;i <= t;i++)
cannam@52: 	{	
cannam@52:             int k = std::min((i+p_pre),static_cast<int>(df.size()));
cannam@52: 		smoothed[i] = mean_array(df,1,k);
cannam@52: 	}
cannam@52: 	// find threshold for bulk of samples across a moving average from [i-p_pre,i+p_post]
cannam@52: 	for (uint i = t+1;i < df.size()-p_post;i++)
cannam@52: 	{
cannam@52: 		smoothed[i] = mean_array(df,i-p_pre,i+p_post);
cannam@52: 	}
cannam@52: 	// for last few samples calculate threshold, again, not enough samples to do as above
cannam@52: 	for (uint i = df.size()-p_post;i < df.size();i++)
cannam@52: 	{
cannam@52:             int k = std::max((static_cast<int> (i) -p_post),1);
cannam@52: 		smoothed[i] = mean_array(df,k,df.size());
cannam@52: 	}
cannam@52: 
cannam@52: 	// subtract the threshold from the detection function and check that it is not less than 0
cannam@52: 	for (uint i = 0;i < df.size();i++)
cannam@52: 	{
cannam@52: 		df[i] -= smoothed[i];
cannam@52: 		if (df[i] < 0)
cannam@52: 		{
cannam@52: 			df[i] = 0;
cannam@52: 		}
cannam@52: 	}
cannam@52: }
cannam@52: 
cannam@52: double
cannam@52: TempoTrackV2::mean_array(const d_vec_t &dfin,int start,int end)
cannam@52: {
cannam@52: 
cannam@52: 	double sum = 0.;
cannam@52: 	
cannam@52: 	// find sum
cannam@52: 	for (int i = start;i < end+1;i++)
cannam@52: 	{
cannam@52: 		sum += dfin[i];
cannam@52: 	}
cannam@52: 
cannam@52: 	return static_cast<double> (sum / (end - start + 1) );	// average and return
cannam@52: }
cannam@52: 
cannam@52: void
cannam@52: TempoTrackV2::filter_df(d_vec_t &df)
cannam@52: {
cannam@52: 
cannam@52: 
cannam@52:   d_vec_t a(3);
cannam@52:   d_vec_t b(3);
cannam@52:   d_vec_t	lp_df(df.size());
cannam@52: 
cannam@52:   //equivalent in matlab to [b,a] = butter(2,0.4);
cannam@52: 	a[0] = 1.0000;
cannam@52: 	a[1] = -0.3695;
cannam@52: 	a[2] = 0.1958;
cannam@52: 	b[0] = 0.2066;
cannam@52: 	b[1] = 0.4131;
cannam@52: 	b[2] = 0.2066;
cannam@52: 	
cannam@52:   double inp1 = 0.;
cannam@52:   double inp2 = 0.;
cannam@52:   double out1 = 0.;
cannam@52:   double out2 = 0.;
cannam@52: 
cannam@52: 
cannam@52:   // forwards filtering
cannam@52: 	for (uint i = 0;i < df.size();i++)
cannam@52: 	{
cannam@52:     lp_df[i] =  b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2;
cannam@52:     inp2 = inp1;
cannam@52:     inp1 = df[i];
cannam@52:     out2 = out1;
cannam@52:     out1 = lp_df[i];
cannam@52: 	}
cannam@52: 
cannam@52: 
cannam@52:   // copy forwards filtering to df...
cannam@52:   // but, time-reversed, ready for backwards filtering
cannam@52: 	for (uint i = 0;i < df.size();i++)
cannam@52: 	{
cannam@52:     df[i] = lp_df[df.size()-i];    
cannam@52:   }
cannam@52: 
cannam@52: 	for (uint i = 0;i < df.size();i++)
cannam@52: 	{
cannam@52:     lp_df[i] = 0.;    
cannam@52:   }
cannam@52: 
cannam@52:   inp1 = 0.; inp2 = 0.;
cannam@52:   out1 = 0.; out2 = 0.;
cannam@52: 
cannam@52:   // backwards filetering on time-reversed df
cannam@52: 	for (uint i = 0;i < df.size();i++)
cannam@52: 	{
cannam@52:     lp_df[i] =  b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2;
cannam@52:     inp2 = inp1;
cannam@52:     inp1 = df[i];
cannam@52:     out2 = out1;
cannam@52:     out1 = lp_df[i];
cannam@52: 	}
cannam@52: 
cannam@52:   // write the re-reversed (i.e. forward) version back to df
cannam@52: 	for (uint i = 0;i < df.size();i++)
cannam@52: 	{
cannam@52:     df[i] = lp_df[df.size()-i];    
cannam@52:   }
cannam@52: 
cannam@52: 
cannam@52: }
cannam@52: 
cannam@52: 
cannam@52: void
cannam@52: TempoTrackV2::calculateBeatPeriod(const d_vec_t &df, d_vec_t &beat_period)
cannam@52: {
cannam@52: 
cannam@52: // to follow matlab.. split into 512 sample frames with a 128 hop size
cannam@52: // calculate the acf,
cannam@52: // then the rcf.. and then stick the rcfs as columns of a matrix
cannam@52: // then call viterbi decoding with weight vector and transition matrix
cannam@52: // and get best path
cannam@52: 
cannam@52:   uint wv_len = 128;
cannam@52:   double rayparam = 43.;
cannam@52: 
cannam@52:   // make rayleigh weighting curve
cannam@52:   d_vec_t wv(wv_len);
cannam@52:   for (uint i=0; i<wv.size(); i++)
cannam@52:   {
cannam@52:     wv[i] = (static_cast<double> (i) / pow(rayparam,2.)) * exp((-1.*pow(-static_cast<double> (i),2.)) / (2.*pow(rayparam,2.)));
cannam@52:   }
cannam@52: 
cannam@52: 
cannam@52:   uint winlen = 512;
cannam@52:   uint step = 128;
cannam@52: 
cannam@52:   d_mat_t rcfmat;
cannam@52:   int col_counter = -1;
cannam@52:   // main loop for beat period calculation
cannam@52:   for (uint i=0; i<(df.size()-winlen); i+=step)
cannam@52:   {
cannam@52:     // get dfframe
cannam@52:     d_vec_t dfframe(winlen);
cannam@52:     for (uint k=0; k<winlen; k++)
cannam@52:     {
cannam@52:       dfframe[k] = df[i+k];
cannam@52:     }
cannam@52:     // get rcf vector for current frame
cannam@52:     d_vec_t rcf(wv_len);    
cannam@52:     get_rcf(dfframe,wv,rcf);
cannam@52:   
cannam@52:     rcfmat.push_back( d_vec_t() ); // adds a new column
cannam@52:     col_counter++;
cannam@52:     for (uint j=0; j<rcf.size(); j++)
cannam@52:     {
cannam@52:       rcfmat[col_counter].push_back( rcf[j] );
cannam@52:     }
cannam@52: 
cannam@52:   }
cannam@52:   
cannam@52:   // now call viterbi decoding function
cannam@52:   viterbi_decode(rcfmat,wv,beat_period);
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52: }
cannam@52: 
cannam@52: 
cannam@52: void
cannam@52: TempoTrackV2::get_rcf(const d_vec_t &dfframe_in, const d_vec_t &wv, d_vec_t &rcf)
cannam@52: {
cannam@52:   // calculate autocorrelation function
cannam@52:   // then rcf
cannam@52:   // just hard code for now... don't really need separate functions to do this
cannam@52: 
cannam@52:   // make acf
cannam@52: 
cannam@52:   d_vec_t dfframe(dfframe_in);
cannam@52: 
cannam@52:   adapt_thresh(dfframe);
cannam@52: 
cannam@52:   d_vec_t acf(dfframe.size());
cannam@52: 
cannam@52: 
cannam@52:   for (uint lag=0; lag<dfframe.size(); lag++)
cannam@52:   {
cannam@52: 
cannam@52:     double sum = 0.;
cannam@52:     double tmp = 0.;
cannam@52: 
cannam@52:     for (uint n=0; n<(dfframe.size()-lag); n++)
cannam@52:     {
cannam@52:       tmp = dfframe[n] * dfframe[n+lag];    
cannam@52:       sum += tmp;
cannam@52:     }
cannam@52:     acf[lag] = static_cast<double> (sum/ (dfframe.size()-lag));
cannam@52:   }
cannam@52: 
cannam@52: 
cannam@52: //  for (uint i=0; i<dfframe.size(); i++)
cannam@52: //  {
cannam@52: //    cout << dfframe[i] << " " << acf[i]  << endl;
cannam@52: //  }
cannam@52: 
cannam@52: //    cout << "~~~~~~~~~~~~~~" << endl;
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52:   // now apply comb filtering
cannam@52: 	int numelem = 4;
cannam@52: 	
cannam@52: //	for (uint i = 1;i < 118;i++) // max beat period
cannam@52: 	for (uint i = 2;i < rcf.size();i++) // max beat period
cannam@52: 	{
cannam@52: 		for (int a = 1;a <= numelem;a++) // number of comb elements
cannam@52: 		{
cannam@52: 			for (int b = 1-a;b <= a-1;b++) // general state using normalisation of comb elements
cannam@52: 			{
cannam@52: 				rcf[i-1] += ( acf[(a*i+b)-1]*wv[i-1] ) / (2.*a-1.);	// calculate value for comb filter row
cannam@52: 			}
cannam@52: 		}
cannam@52: 	}
cannam@52:   
cannam@52:   // apply adaptive threshold to rcf
cannam@52:   adapt_thresh(rcf);
cannam@52:   
cannam@52:   double rcfsum =0.;
cannam@52:   for (uint i=0; i<rcf.size(); i++)
cannam@52:   {
cannam@52:  //  rcf[i] *= acf[i];
cannam@52:     rcf[i] += EPS ;
cannam@52:     rcfsum += rcf[i];
cannam@52:   }
cannam@52: 
cannam@52:   // normalise rcf to sum to unity
cannam@52:   for (uint i=0; i<rcf.size(); i++)
cannam@52:   {
cannam@52:     rcf[i] /= (rcfsum + EPS);
cannam@52:   }
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52: }
cannam@52: 
cannam@52: void
cannam@52: TempoTrackV2::viterbi_decode(const d_mat_t &rcfmat, const d_vec_t &wv, d_vec_t &beat_period)
cannam@52: {
cannam@52: 
cannam@52:   // make transition matrix
cannam@52: 	d_mat_t tmat;
cannam@52: 	for (uint i=0;i<wv.size();i++)
cannam@52: 	{
cannam@52: 		tmat.push_back ( d_vec_t() ); // adds a new column
cannam@52: 		for (uint j=0; j<wv.size(); j++)
cannam@52: 		{	
cannam@52: 			tmat[i].push_back(0.); // fill with zeros initially
cannam@52: 		}
cannam@52: 	}
cannam@52: 
cannam@52:   double sigma = 8.;
cannam@52: 	for (uint i=20;i <wv.size()-20; i++)
cannam@52: 	{
cannam@52: 		for (uint j=20; j<wv.size()-20; j++)
cannam@52: 		{	
cannam@52:       double mu = static_cast<double>(i);
cannam@52: 			tmat[i][j] = exp( (-1.*pow((j-mu),2.)) / (2.*pow(sigma,2.)) );
cannam@52: 		}
cannam@52: 	}
cannam@52: 
cannam@52:   d_mat_t delta;
cannam@52:   i_mat_t psi;
cannam@52: 	for (uint i=0;i <rcfmat.size(); i++)
cannam@52: 	{
cannam@52:     delta.push_back( d_vec_t());
cannam@52:     psi.push_back( i_vec_t());
cannam@52: 		for (uint j=0; j<rcfmat[i].size(); j++)
cannam@52: 		{	
cannam@52: 			delta[i].push_back(0.); // fill with zeros initially
cannam@52: 			psi[i].push_back(0); // fill with zeros initially
cannam@52: 		}
cannam@52: 	}
cannam@52: 
cannam@52: 
cannam@52:   uint T = delta.size();
cannam@52:   uint Q = delta[0].size();
cannam@52: 
cannam@52:   // initialize first column of delta
cannam@52:   for (uint j=0; j<Q; j++)
cannam@52:   {
cannam@52:     delta[0][j] = wv[j] * rcfmat[0][j];
cannam@52:     psi[0][j] = 0;
cannam@52:   }
cannam@52: 
cannam@52:   double deltasum = 0.;
cannam@52:   for (uint i=0; i<Q; i++)
cannam@52:   {
cannam@52:     deltasum += delta[0][i];
cannam@52:   }      
cannam@52:   for (uint i=0; i<Q; i++)
cannam@52:   {
cannam@52:     delta[0][i] /= (deltasum + EPS);
cannam@52:   }      
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52:   for (uint t=1; t<T; t++)
cannam@52:   {
cannam@52:     d_vec_t tmp_vec(Q);
cannam@52: 
cannam@52:     for (uint j=0; j<Q; j++)
cannam@52:     {
cannam@52: 
cannam@52:       for (uint i=0; i<Q; i++)
cannam@52:       {
cannam@52:         tmp_vec[i] = delta[t-1][i] * tmat[j][i];
cannam@52:       }      
cannam@52:    
cannam@52:       delta[t][j] = get_max_val(tmp_vec);    
cannam@52: 
cannam@52:       psi[t][j] = get_max_ind(tmp_vec);
cannam@52:  
cannam@52:       delta[t][j] *= rcfmat[t][j];
cannam@52: 
cannam@52: 
cannam@52:     }
cannam@52: 
cannam@52:     double deltasum = 0.;
cannam@52:     for (uint i=0; i<Q; i++)
cannam@52:     {
cannam@52:       deltasum += delta[t][i];
cannam@52:     }      
cannam@52:     for (uint i=0; i<Q; i++)
cannam@52:     {
cannam@52:       delta[t][i] /= (deltasum + EPS);
cannam@52:     }      
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52: 
cannam@52:   }
cannam@52: 
cannam@52: 
cannam@52: //  ofstream tmatfile;
cannam@52: //  tmatfile.open("/home/matthewd/Desktop/tmat.txt");
cannam@52: 
cannam@52: // 	for (uint i=0;i <delta.size(); i++)
cannam@52: //	{
cannam@52: //		for (uint j=0; j<delta[i].size(); j++)
cannam@52: //		{	
cannam@52: //      tmatfile << rcfmat[i][j] << endl;
cannam@52: //		}
cannam@52: //	}
cannam@52: 
cannam@52: //  tmatfile.close();
cannam@52: 
cannam@52:   i_vec_t bestpath(T);
cannam@52:   d_vec_t tmp_vec(Q);
cannam@52:   for (uint i=0; i<Q; i++)
cannam@52:   {  
cannam@52:     tmp_vec[i] = delta[T-1][i];
cannam@52:   }
cannam@52: 
cannam@52:   
cannam@52:   bestpath[T-1] = get_max_ind(tmp_vec);
cannam@52:  
cannam@52:   for (uint t=T-2; t>0 ;t--)
cannam@52:   {
cannam@52:     bestpath[t] = psi[t+1][bestpath[t+1]];
cannam@52:   }
cannam@52:   // very weird hack!
cannam@52:   bestpath[0] = psi[1][bestpath[1]];
cannam@52: 
cannam@52: //  for (uint i=0; i<bestpath.size(); i++)
cannam@52: //  {
cannam@52: //    cout << bestpath[i] << endl;
cannam@52: //  }
cannam@52: 
cannam@52: 
cannam@52:   uint lastind = 0;
cannam@52:   for (uint i=0; i<T; i++)
cannam@52:   {  
cannam@52:     uint step = 128;
cannam@52:  //   cout << bestpath[i] << " " << i << endl;
cannam@52:     for (uint j=0; j<step; j++)
cannam@52:     {
cannam@52:       lastind = i*step+j;
cannam@52:       beat_period[lastind] = bestpath[i];
cannam@52:       
cannam@52:     }
cannam@52:   }
cannam@52: 
cannam@52:   //fill in the last values...
cannam@52:   for (uint i=lastind; i<beat_period.size(); i++)
cannam@52:   {
cannam@52:     beat_period[i] = beat_period[lastind];
cannam@52:   }
cannam@52:  
cannam@52: 
cannam@52: 
cannam@52: }
cannam@52: 
cannam@52: double
cannam@52: TempoTrackV2::get_max_val(const d_vec_t &df)
cannam@52: {
cannam@52:   double maxval = 0.;
cannam@52:   for (uint i=0; i<df.size(); i++)
cannam@52:   {
cannam@52: 
cannam@52:     if (maxval < df[i])
cannam@52:     {
cannam@52:       maxval = df[i];
cannam@52:     }
cannam@52: 
cannam@52:   }
cannam@52: 
cannam@52:     
cannam@52:   return maxval;
cannam@52: 
cannam@52: }
cannam@52: 
cannam@52: int
cannam@52: TempoTrackV2::get_max_ind(const d_vec_t &df)
cannam@52: {
cannam@52: 
cannam@52:   double maxval = 0.;
cannam@52:   int ind = 0;
cannam@52:   for (uint i=0; i<df.size(); i++)
cannam@52:   {
cannam@52:     if (maxval < df[i])
cannam@52:     {
cannam@52:       maxval = df[i];
cannam@52:       ind = i;
cannam@52:     }
cannam@52: 
cannam@52:   }
cannam@52:   
cannam@52:   return ind;
cannam@52: 
cannam@52: }
cannam@52: 
cannam@52: void
cannam@52: TempoTrackV2::normalise_vec(d_vec_t &df)
cannam@52: {
cannam@52:   double sum = 0.;
cannam@52:   for (uint i=0; i<df.size(); i++)
cannam@52:   {
cannam@52:     sum += df[i];
cannam@52:   }
cannam@52: 
cannam@52:   for (uint i=0; i<df.size(); i++)
cannam@52:   {
cannam@52:     df[i]/= (sum + EPS);
cannam@52:   }
cannam@52: 
cannam@52: 
cannam@52: }
cannam@52: 
cannam@52: void
cannam@52: TempoTrackV2::calculateBeats(const d_vec_t &df, const d_vec_t &beat_period,
cannam@52:                              d_vec_t &beats)
cannam@52: {
cannam@52: 
cannam@52:   d_vec_t cumscore(df.size());
cannam@52:   i_vec_t backlink(df.size());
cannam@52:   d_vec_t localscore(df.size());
cannam@52: 
cannam@52:   // WHEN I FIGURE OUT HOW, I'LL WANT TO DO SOME FILTERING ON THIS... 
cannam@52:   for (uint i=0; i<df.size(); i++)
cannam@52:   {
cannam@52:     localscore[i] = df[i];
cannam@52:     backlink[i] = -1;
cannam@52:   }
cannam@52: 
cannam@52:   double tightness = 4.;
cannam@52:   double alpha = 0.9;
cannam@52: 
cannam@52:   // main loop
cannam@52:   for (uint i=3*beat_period[0]; i<localscore.size(); i++)
cannam@52:   {
cannam@52:     int prange_min = -2*beat_period[i];
cannam@52:     int prange_max = round(-0.5*beat_period[i]);
cannam@52: 
cannam@52:     d_vec_t txwt (prange_max - prange_min + 1);
cannam@52:     d_vec_t scorecands (txwt.size());
cannam@52: 
cannam@52:     for (uint j=0;j<txwt.size();j++)
cannam@52:     {
cannam@52:       double mu = static_cast<double> (beat_period[i]);
cannam@52:       txwt[j] = exp( -0.5*pow(tightness * log((round(2*mu)-j)/mu),2));
cannam@52: 
cannam@52:       scorecands[j] = txwt[j] * cumscore[i+prange_min+j];
cannam@52:     }
cannam@52: 
cannam@52:     double vv = get_max_val(scorecands);
cannam@52:     int xx = get_max_ind(scorecands);
cannam@52: 
cannam@52:     cumscore[i] = alpha*vv + (1.-alpha)*localscore[i];
cannam@52: 
cannam@52:     backlink[i] = i+prange_min+xx;
cannam@52: 
cannam@52:   }
cannam@52: 
cannam@52: 
cannam@52:   d_vec_t tmp_vec;
cannam@52:   for (uint i=cumscore.size() - beat_period[beat_period.size()-1] ; i<cumscore.size(); i++)
cannam@52:   {
cannam@52:     tmp_vec.push_back(cumscore[i]);
cannam@52:   }  
cannam@52: 
cannam@52:   int startpoint = get_max_ind(tmp_vec) + cumscore.size() - beat_period[beat_period.size()-1] ;
cannam@52: 
cannam@52:   i_vec_t ibeats;
cannam@52:   ibeats.push_back(startpoint);
cannam@52:   while (backlink[ibeats.back()] > 3*beat_period[0])
cannam@52:   {
cannam@52:     ibeats.push_back(backlink[ibeats.back()]);
cannam@52:   }
cannam@52:   
cannam@52: 
cannam@52:   for (uint i=0; i<ibeats.size(); i++)
cannam@52:   { 
cannam@52: 
cannam@52:     beats.push_back( static_cast<double>(ibeats[i]) );
cannam@52: 
cannam@52:  //   cout << ibeats[i] << " "  << beats[i] <<endl;
cannam@52:   }
cannam@52: }
cannam@52: 
cannam@52: