comparison dsp/tempotracking/TempoTrackV2.cpp @ 70:c3cdb404f807

* build fixes for win32-x-g++
author cannam
date Tue, 02 Jun 2009 09:53:01 +0000
parents d241e7701c0c
children 054c384d860d
comparison
equal deleted inserted replaced
69:6afa0e011f74 70:c3cdb404f807
42 double out1 = 0.; 42 double out1 = 0.;
43 double out2 = 0.; 43 double out2 = 0.;
44 44
45 45
46 // forwards filtering 46 // forwards filtering
47 for (uint i = 0;i < df.size();i++) 47 for (unsigned int i = 0;i < df.size();i++)
48 { 48 {
49 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2; 49 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2;
50 inp2 = inp1; 50 inp2 = inp1;
51 inp1 = df[i]; 51 inp1 = df[i];
52 out2 = out1; 52 out2 = out1;
53 out1 = lp_df[i]; 53 out1 = lp_df[i];
54 } 54 }
55 55
56 // copy forwards filtering to df... 56 // copy forwards filtering to df...
57 // but, time-reversed, ready for backwards filtering 57 // but, time-reversed, ready for backwards filtering
58 for (uint i = 0;i < df.size();i++) 58 for (unsigned int i = 0;i < df.size();i++)
59 { 59 {
60 df[i] = lp_df[df.size()-i-1]; 60 df[i] = lp_df[df.size()-i-1];
61 } 61 }
62 62
63 for (uint i = 0;i < df.size();i++) 63 for (unsigned int i = 0;i < df.size();i++)
64 { 64 {
65 lp_df[i] = 0.; 65 lp_df[i] = 0.;
66 } 66 }
67 67
68 inp1 = 0.; inp2 = 0.; 68 inp1 = 0.; inp2 = 0.;
69 out1 = 0.; out2 = 0.; 69 out1 = 0.; out2 = 0.;
70 70
71 // backwards filetering on time-reversed df 71 // backwards filetering on time-reversed df
72 for (uint i = 0;i < df.size();i++) 72 for (unsigned int i = 0;i < df.size();i++)
73 { 73 {
74 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2; 74 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2;
75 inp2 = inp1; 75 inp2 = inp1;
76 inp1 = df[i]; 76 inp1 = df[i];
77 out2 = out1; 77 out2 = out1;
78 out1 = lp_df[i]; 78 out1 = lp_df[i];
79 } 79 }
80 80
81 // write the re-reversed (i.e. forward) version back to df 81 // write the re-reversed (i.e. forward) version back to df
82 for (uint i = 0;i < df.size();i++) 82 for (unsigned int i = 0;i < df.size();i++)
83 { 83 {
84 df[i] = lp_df[df.size()-i-1]; 84 df[i] = lp_df[df.size()-i-1];
85 } 85 }
86 } 86 }
87 87
94 // calculate the acf, 94 // calculate the acf,
95 // then the rcf.. and then stick the rcfs as columns of a matrix 95 // then the rcf.. and then stick the rcfs as columns of a matrix
96 // then call viterbi decoding with weight vector and transition matrix 96 // then call viterbi decoding with weight vector and transition matrix
97 // and get best path 97 // and get best path
98 98
99 uint wv_len = 128; 99 unsigned int wv_len = 128;
100 double rayparam = 43.; 100 double rayparam = 43.;
101 101
102 // make rayleigh weighting curve 102 // make rayleigh weighting curve
103 d_vec_t wv(wv_len); 103 d_vec_t wv(wv_len);
104 for (uint i=0; i<wv.size(); i++) 104 for (unsigned int i=0; i<wv.size(); i++)
105 { 105 {
106 wv[i] = (static_cast<double> (i) / pow(rayparam,2.)) * exp((-1.*pow(-static_cast<double> (i),2.)) / (2.*pow(rayparam,2.))); 106 wv[i] = (static_cast<double> (i) / pow(rayparam,2.)) * exp((-1.*pow(-static_cast<double> (i),2.)) / (2.*pow(rayparam,2.)));
107 } 107 }
108 108
109 // beat tracking frame size (roughly 6 seconds) and hop (1.5 seconds) 109 // beat tracking frame size (roughly 6 seconds) and hop (1.5 seconds)
110 uint winlen = 512; 110 unsigned int winlen = 512;
111 uint step = 128; 111 unsigned int step = 128;
112 112
113 // matrix to store output of comb filter bank, increment column of matrix at each frame 113 // matrix to store output of comb filter bank, increment column of matrix at each frame
114 d_mat_t rcfmat; 114 d_mat_t rcfmat;
115 int col_counter = -1; 115 int col_counter = -1;
116 116
117 // main loop for beat period calculation 117 // main loop for beat period calculation
118 for (uint i=0; i+winlen<df.size(); i+=step) 118 for (unsigned int i=0; i+winlen<df.size(); i+=step)
119 { 119 {
120 // get dfframe 120 // get dfframe
121 d_vec_t dfframe(winlen); 121 d_vec_t dfframe(winlen);
122 for (uint k=0; k<winlen; k++) 122 for (unsigned int k=0; k<winlen; k++)
123 { 123 {
124 dfframe[k] = df[i+k]; 124 dfframe[k] = df[i+k];
125 } 125 }
126 // get rcf vector for current frame 126 // get rcf vector for current frame
127 d_vec_t rcf(wv_len); 127 d_vec_t rcf(wv_len);
128 get_rcf(dfframe,wv,rcf); 128 get_rcf(dfframe,wv,rcf);
129 129
130 rcfmat.push_back( d_vec_t() ); // adds a new column 130 rcfmat.push_back( d_vec_t() ); // adds a new column
131 col_counter++; 131 col_counter++;
132 for (uint j=0; j<rcf.size(); j++) 132 for (unsigned int j=0; j<rcf.size(); j++)
133 { 133 {
134 rcfmat[col_counter].push_back( rcf[j] ); 134 rcfmat[col_counter].push_back( rcf[j] );
135 } 135 }
136 } 136 }
137 137
154 MathUtilities::adaptiveThreshold(dfframe); 154 MathUtilities::adaptiveThreshold(dfframe);
155 155
156 d_vec_t acf(dfframe.size()); 156 d_vec_t acf(dfframe.size());
157 157
158 158
159 for (uint lag=0; lag<dfframe.size(); lag++) 159 for (unsigned int lag=0; lag<dfframe.size(); lag++)
160 { 160 {
161 double sum = 0.; 161 double sum = 0.;
162 double tmp = 0.; 162 double tmp = 0.;
163 163
164 for (uint n=0; n<(dfframe.size()-lag); n++) 164 for (unsigned int n=0; n<(dfframe.size()-lag); n++)
165 { 165 {
166 tmp = dfframe[n] * dfframe[n+lag]; 166 tmp = dfframe[n] * dfframe[n+lag];
167 sum += tmp; 167 sum += tmp;
168 } 168 }
169 acf[lag] = static_cast<double> (sum/ (dfframe.size()-lag)); 169 acf[lag] = static_cast<double> (sum/ (dfframe.size()-lag));
170 } 170 }
171 171
172 // now apply comb filtering 172 // now apply comb filtering
173 int numelem = 4; 173 int numelem = 4;
174 174
175 for (uint i = 2;i < rcf.size();i++) // max beat period 175 for (unsigned int i = 2;i < rcf.size();i++) // max beat period
176 { 176 {
177 for (int a = 1;a <= numelem;a++) // number of comb elements 177 for (int a = 1;a <= numelem;a++) // number of comb elements
178 { 178 {
179 for (int b = 1-a;b <= a-1;b++) // general state using normalisation of comb elements 179 for (int b = 1-a;b <= a-1;b++) // general state using normalisation of comb elements
180 { 180 {
185 185
186 // apply adaptive threshold to rcf 186 // apply adaptive threshold to rcf
187 MathUtilities::adaptiveThreshold(rcf); 187 MathUtilities::adaptiveThreshold(rcf);
188 188
189 double rcfsum =0.; 189 double rcfsum =0.;
190 for (uint i=0; i<rcf.size(); i++) 190 for (unsigned int i=0; i<rcf.size(); i++)
191 { 191 {
192 rcf[i] += EPS ; 192 rcf[i] += EPS ;
193 rcfsum += rcf[i]; 193 rcfsum += rcf[i];
194 } 194 }
195 195
196 // normalise rcf to sum to unity 196 // normalise rcf to sum to unity
197 for (uint i=0; i<rcf.size(); i++) 197 for (unsigned int i=0; i<rcf.size(); i++)
198 { 198 {
199 rcf[i] /= (rcfsum + EPS); 199 rcf[i] /= (rcfsum + EPS);
200 } 200 }
201 } 201 }
202 202
206 // following Kevin Murphy's Viterbi decoding to get best path of 206 // following Kevin Murphy's Viterbi decoding to get best path of
207 // beat periods through rfcmat 207 // beat periods through rfcmat
208 208
209 // make transition matrix 209 // make transition matrix
210 d_mat_t tmat; 210 d_mat_t tmat;
211 for (uint i=0;i<wv.size();i++) 211 for (unsigned int i=0;i<wv.size();i++)
212 { 212 {
213 tmat.push_back ( d_vec_t() ); // adds a new column 213 tmat.push_back ( d_vec_t() ); // adds a new column
214 for (uint j=0; j<wv.size(); j++) 214 for (unsigned int j=0; j<wv.size(); j++)
215 { 215 {
216 tmat[i].push_back(0.); // fill with zeros initially 216 tmat[i].push_back(0.); // fill with zeros initially
217 } 217 }
218 } 218 }
219 219
220 // variance of Gaussians in transition matrix 220 // variance of Gaussians in transition matrix
221 // formed of Gaussians on diagonal - implies slow tempo change 221 // formed of Gaussians on diagonal - implies slow tempo change
222 double sigma = 8.; 222 double sigma = 8.;
223 // don't want really short beat periods, or really long ones 223 // don't want really short beat periods, or really long ones
224 for (uint i=20;i <wv.size()-20; i++) 224 for (unsigned int i=20;i <wv.size()-20; i++)
225 { 225 {
226 for (uint j=20; j<wv.size()-20; j++) 226 for (unsigned int j=20; j<wv.size()-20; j++)
227 { 227 {
228 double mu = static_cast<double>(i); 228 double mu = static_cast<double>(i);
229 tmat[i][j] = exp( (-1.*pow((j-mu),2.)) / (2.*pow(sigma,2.)) ); 229 tmat[i][j] = exp( (-1.*pow((j-mu),2.)) / (2.*pow(sigma,2.)) );
230 } 230 }
231 } 231 }
233 // parameters for Viterbi decoding... this part is taken from 233 // parameters for Viterbi decoding... this part is taken from
234 // Murphy's matlab 234 // Murphy's matlab
235 235
236 d_mat_t delta; 236 d_mat_t delta;
237 i_mat_t psi; 237 i_mat_t psi;
238 for (uint i=0;i <rcfmat.size(); i++) 238 for (unsigned int i=0;i <rcfmat.size(); i++)
239 { 239 {
240 delta.push_back( d_vec_t()); 240 delta.push_back( d_vec_t());
241 psi.push_back( i_vec_t()); 241 psi.push_back( i_vec_t());
242 for (uint j=0; j<rcfmat[i].size(); j++) 242 for (unsigned int j=0; j<rcfmat[i].size(); j++)
243 { 243 {
244 delta[i].push_back(0.); // fill with zeros initially 244 delta[i].push_back(0.); // fill with zeros initially
245 psi[i].push_back(0); // fill with zeros initially 245 psi[i].push_back(0); // fill with zeros initially
246 } 246 }
247 } 247 }
248 248
249 249
250 uint T = delta.size(); 250 unsigned int T = delta.size();
251 251
252 if (T < 2) return; // can't do anything at all meaningful 252 if (T < 2) return; // can't do anything at all meaningful
253 253
254 uint Q = delta[0].size(); 254 unsigned int Q = delta[0].size();
255 255
256 // initialize first column of delta 256 // initialize first column of delta
257 for (uint j=0; j<Q; j++) 257 for (unsigned int j=0; j<Q; j++)
258 { 258 {
259 delta[0][j] = wv[j] * rcfmat[0][j]; 259 delta[0][j] = wv[j] * rcfmat[0][j];
260 psi[0][j] = 0; 260 psi[0][j] = 0;
261 } 261 }
262 262
263 double deltasum = 0.; 263 double deltasum = 0.;
264 for (uint i=0; i<Q; i++) 264 for (unsigned int i=0; i<Q; i++)
265 { 265 {
266 deltasum += delta[0][i]; 266 deltasum += delta[0][i];
267 } 267 }
268 for (uint i=0; i<Q; i++) 268 for (unsigned int i=0; i<Q; i++)
269 { 269 {
270 delta[0][i] /= (deltasum + EPS); 270 delta[0][i] /= (deltasum + EPS);
271 } 271 }
272 272
273 273
274 for (uint t=1; t<T; t++) 274 for (unsigned int t=1; t<T; t++)
275 { 275 {
276 d_vec_t tmp_vec(Q); 276 d_vec_t tmp_vec(Q);
277 277
278 for (uint j=0; j<Q; j++) 278 for (unsigned int j=0; j<Q; j++)
279 { 279 {
280 for (uint i=0; i<Q; i++) 280 for (unsigned int i=0; i<Q; i++)
281 { 281 {
282 tmp_vec[i] = delta[t-1][i] * tmat[j][i]; 282 tmp_vec[i] = delta[t-1][i] * tmat[j][i];
283 } 283 }
284 284
285 delta[t][j] = get_max_val(tmp_vec); 285 delta[t][j] = get_max_val(tmp_vec);
289 delta[t][j] *= rcfmat[t][j]; 289 delta[t][j] *= rcfmat[t][j];
290 } 290 }
291 291
292 // normalise current delta column 292 // normalise current delta column
293 double deltasum = 0.; 293 double deltasum = 0.;
294 for (uint i=0; i<Q; i++) 294 for (unsigned int i=0; i<Q; i++)
295 { 295 {
296 deltasum += delta[t][i]; 296 deltasum += delta[t][i];
297 } 297 }
298 for (uint i=0; i<Q; i++) 298 for (unsigned int i=0; i<Q; i++)
299 { 299 {
300 delta[t][i] /= (deltasum + EPS); 300 delta[t][i] /= (deltasum + EPS);
301 } 301 }
302 } 302 }
303 303
304 i_vec_t bestpath(T); 304 i_vec_t bestpath(T);
305 d_vec_t tmp_vec(Q); 305 d_vec_t tmp_vec(Q);
306 for (uint i=0; i<Q; i++) 306 for (unsigned int i=0; i<Q; i++)
307 { 307 {
308 tmp_vec[i] = delta[T-1][i]; 308 tmp_vec[i] = delta[T-1][i];
309 } 309 }
310 310
311 // find starting point - best beat period for "last" frame 311 // find starting point - best beat period for "last" frame
312 bestpath[T-1] = get_max_ind(tmp_vec); 312 bestpath[T-1] = get_max_ind(tmp_vec);
313 313
314 // backtrace through index of maximum values in psi 314 // backtrace through index of maximum values in psi
315 for (uint t=T-2; t>0 ;t--) 315 for (unsigned int t=T-2; t>0 ;t--)
316 { 316 {
317 bestpath[t] = psi[t+1][bestpath[t+1]]; 317 bestpath[t] = psi[t+1][bestpath[t+1]];
318 } 318 }
319 319
320 // weird but necessary hack -- couldn't get above loop to terminate at t >= 0 320 // weird but necessary hack -- couldn't get above loop to terminate at t >= 0
321 bestpath[0] = psi[1][bestpath[1]]; 321 bestpath[0] = psi[1][bestpath[1]];
322 322
323 uint lastind = 0; 323 unsigned int lastind = 0;
324 for (uint i=0; i<T; i++) 324 for (unsigned int i=0; i<T; i++)
325 { 325 {
326 uint step = 128; 326 unsigned int step = 128;
327 for (uint j=0; j<step; j++) 327 for (unsigned int j=0; j<step; j++)
328 { 328 {
329 lastind = i*step+j; 329 lastind = i*step+j;
330 beat_period[lastind] = bestpath[i]; 330 beat_period[lastind] = bestpath[i];
331 } 331 }
332 // std::cerr << "bestpath[" << i << "] = " << bestpath[i] << " (used for beat_periods " << i*step << " to " << i*step+step-1 << ")" << std::endl; 332 // std::cerr << "bestpath[" << i << "] = " << bestpath[i] << " (used for beat_periods " << i*step << " to " << i*step+step-1 << ")" << std::endl;
333 } 333 }
334 334
335 //fill in the last values... 335 //fill in the last values...
336 for (uint i=lastind; i<beat_period.size(); i++) 336 for (unsigned int i=lastind; i<beat_period.size(); i++)
337 { 337 {
338 beat_period[i] = beat_period[lastind]; 338 beat_period[i] = beat_period[lastind];
339 } 339 }
340 340
341 for (uint i = 0; i < beat_period.size(); i++) 341 for (unsigned int i = 0; i < beat_period.size(); i++)
342 { 342 {
343 tempi.push_back((60. * m_rate / m_increment)/beat_period[i]); 343 tempi.push_back((60. * m_rate / m_increment)/beat_period[i]);
344 } 344 }
345 } 345 }
346 346
347 double 347 double
348 TempoTrackV2::get_max_val(const d_vec_t &df) 348 TempoTrackV2::get_max_val(const d_vec_t &df)
349 { 349 {
350 double maxval = 0.; 350 double maxval = 0.;
351 for (uint i=0; i<df.size(); i++) 351 for (unsigned int i=0; i<df.size(); i++)
352 { 352 {
353 if (maxval < df[i]) 353 if (maxval < df[i])
354 { 354 {
355 maxval = df[i]; 355 maxval = df[i];
356 } 356 }
362 int 362 int
363 TempoTrackV2::get_max_ind(const d_vec_t &df) 363 TempoTrackV2::get_max_ind(const d_vec_t &df)
364 { 364 {
365 double maxval = 0.; 365 double maxval = 0.;
366 int ind = 0; 366 int ind = 0;
367 for (uint i=0; i<df.size(); i++) 367 for (unsigned int i=0; i<df.size(); i++)
368 { 368 {
369 if (maxval < df[i]) 369 if (maxval < df[i])
370 { 370 {
371 maxval = df[i]; 371 maxval = df[i];
372 ind = i; 372 ind = i;
378 378
379 void 379 void
380 TempoTrackV2::normalise_vec(d_vec_t &df) 380 TempoTrackV2::normalise_vec(d_vec_t &df)
381 { 381 {
382 double sum = 0.; 382 double sum = 0.;
383 for (uint i=0; i<df.size(); i++) 383 for (unsigned int i=0; i<df.size(); i++)
384 { 384 {
385 sum += df[i]; 385 sum += df[i];
386 } 386 }
387 387
388 for (uint i=0; i<df.size(); i++) 388 for (unsigned int i=0; i<df.size(); i++)
389 { 389 {
390 df[i]/= (sum + EPS); 390 df[i]/= (sum + EPS);
391 } 391 }
392 } 392 }
393 393
399 399
400 d_vec_t cumscore(df.size()); // store cumulative score 400 d_vec_t cumscore(df.size()); // store cumulative score
401 i_vec_t backlink(df.size()); // backlink (stores best beat locations at each time instant) 401 i_vec_t backlink(df.size()); // backlink (stores best beat locations at each time instant)
402 d_vec_t localscore(df.size()); // localscore, for now this is the same as the detection function 402 d_vec_t localscore(df.size()); // localscore, for now this is the same as the detection function
403 403
404 for (uint i=0; i<df.size(); i++) 404 for (unsigned int i=0; i<df.size(); i++)
405 { 405 {
406 localscore[i] = df[i]; 406 localscore[i] = df[i];
407 backlink[i] = -1; 407 backlink[i] = -1;
408 } 408 }
409 409
410 double tightness = 4.; 410 double tightness = 4.;
411 double alpha = 0.9; 411 double alpha = 0.9;
412 412
413 // main loop 413 // main loop
414 for (uint i=0; i<localscore.size(); i++) 414 for (unsigned int i=0; i<localscore.size(); i++)
415 { 415 {
416 int prange_min = -2*beat_period[i]; 416 int prange_min = -2*beat_period[i];
417 int prange_max = round(-0.5*beat_period[i]); 417 int prange_max = round(-0.5*beat_period[i]);
418 418
419 // transition range 419 // transition range
420 d_vec_t txwt (prange_max - prange_min + 1); 420 d_vec_t txwt (prange_max - prange_min + 1);
421 d_vec_t scorecands (txwt.size()); 421 d_vec_t scorecands (txwt.size());
422 422
423 for (uint j=0;j<txwt.size();j++) 423 for (unsigned int j=0;j<txwt.size();j++)
424 { 424 {
425 double mu = static_cast<double> (beat_period[i]); 425 double mu = static_cast<double> (beat_period[i]);
426 txwt[j] = exp( -0.5*pow(tightness * log((round(2*mu)-j)/mu),2)); 426 txwt[j] = exp( -0.5*pow(tightness * log((round(2*mu)-j)/mu),2));
427 427
428 // IF IN THE ALLOWED RANGE, THEN LOOK AT CUMSCORE[I+PRANGE_MIN+J 428 // IF IN THE ALLOWED RANGE, THEN LOOK AT CUMSCORE[I+PRANGE_MIN+J
445 // std::cerr << "backlink[" << i << "] <= " << backlink[i] << std::endl; 445 // std::cerr << "backlink[" << i << "] <= " << backlink[i] << std::endl;
446 } 446 }
447 447
448 // STARTING POINT, I.E. LAST BEAT.. PICK A STRONG POINT IN cumscore VECTOR 448 // STARTING POINT, I.E. LAST BEAT.. PICK A STRONG POINT IN cumscore VECTOR
449 d_vec_t tmp_vec; 449 d_vec_t tmp_vec;
450 for (uint i=cumscore.size() - beat_period[beat_period.size()-1] ; i<cumscore.size(); i++) 450 for (unsigned int i=cumscore.size() - beat_period[beat_period.size()-1] ; i<cumscore.size(); i++)
451 { 451 {
452 tmp_vec.push_back(cumscore[i]); 452 tmp_vec.push_back(cumscore[i]);
453 } 453 }
454 454
455 int startpoint = get_max_ind(tmp_vec) + cumscore.size() - beat_period[beat_period.size()-1] ; 455 int startpoint = get_max_ind(tmp_vec) + cumscore.size() - beat_period[beat_period.size()-1] ;
469 if (backlink[b] == b) break; // shouldn't happen... haha 469 if (backlink[b] == b) break; // shouldn't happen... haha
470 ibeats.push_back(backlink[b]); 470 ibeats.push_back(backlink[b]);
471 } 471 }
472 472
473 // REVERSE SEQUENCE OF IBEATS AND STORE AS BEATS 473 // REVERSE SEQUENCE OF IBEATS AND STORE AS BEATS
474 for (uint i=0; i<ibeats.size(); i++) 474 for (unsigned int i=0; i<ibeats.size(); i++)
475 { 475 {
476 beats.push_back( static_cast<double>(ibeats[ibeats.size()-i-1]) ); 476 beats.push_back( static_cast<double>(ibeats[ibeats.size()-i-1]) );
477 } 477 }
478 } 478 }
479 479