Mercurial > hg > qm-dsp
comparison dsp/tempotracking/TempoTrackV2.cpp @ 479:7e52c034cf62
Untabify, indent, tidy
author | Chris Cannam <cannam@all-day-breakfast.com> |
---|---|
date | Fri, 31 May 2019 10:35:08 +0100 |
parents | 3f7b4b282df1 |
children | bb78ca3fe7de |
comparison
equal
deleted
inserted
replaced
478:c92718cc6ef1 | 479:7e52c034cf62 |
---|---|
23 | 23 |
24 #define EPS 0.0000008 // just some arbitrary small number | 24 #define EPS 0.0000008 // just some arbitrary small number |
25 | 25 |
26 TempoTrackV2::TempoTrackV2(float rate, size_t increment) : | 26 TempoTrackV2::TempoTrackV2(float rate, size_t increment) : |
27 m_rate(rate), m_increment(increment) { } | 27 m_rate(rate), m_increment(increment) { } |
28 | |
28 TempoTrackV2::~TempoTrackV2() { } | 29 TempoTrackV2::~TempoTrackV2() { } |
29 | 30 |
30 void | 31 void |
31 TempoTrackV2::filter_df(d_vec_t &df) | 32 TempoTrackV2::filter_df(d_vec_t &df) |
32 { | 33 { |
33 d_vec_t a(3); | 34 d_vec_t a(3); |
34 d_vec_t b(3); | 35 d_vec_t b(3); |
35 d_vec_t lp_df(df.size()); | 36 d_vec_t lp_df(df.size()); |
36 | 37 |
37 //equivalent in matlab to [b,a] = butter(2,0.4); | 38 //equivalent in matlab to [b,a] = butter(2,0.4); |
38 a[0] = 1.0000; | 39 a[0] = 1.0000; |
39 a[1] = -0.3695; | 40 a[1] = -0.3695; |
40 a[2] = 0.1958; | 41 a[2] = 0.1958; |
47 double out1 = 0.; | 48 double out1 = 0.; |
48 double out2 = 0.; | 49 double out2 = 0.; |
49 | 50 |
50 | 51 |
51 // forwards filtering | 52 // forwards filtering |
52 for (unsigned int i = 0;i < df.size();i++) | 53 for (unsigned int i = 0;i < df.size();i++) { |
53 { | |
54 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2; | 54 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2; |
55 inp2 = inp1; | 55 inp2 = inp1; |
56 inp1 = df[i]; | 56 inp1 = df[i]; |
57 out2 = out1; | 57 out2 = out1; |
58 out1 = lp_df[i]; | 58 out1 = lp_df[i]; |
59 } | 59 } |
60 | 60 |
61 // copy forwards filtering to df... | 61 // copy forwards filtering to df... |
62 // but, time-reversed, ready for backwards filtering | 62 // but, time-reversed, ready for backwards filtering |
63 for (unsigned int i = 0;i < df.size();i++) | 63 for (unsigned int i = 0;i < df.size();i++) { |
64 { | |
65 df[i] = lp_df[df.size()-i-1]; | 64 df[i] = lp_df[df.size()-i-1]; |
66 } | 65 } |
67 | 66 |
68 for (unsigned int i = 0;i < df.size();i++) | 67 for (unsigned int i = 0;i < df.size();i++) { |
69 { | |
70 lp_df[i] = 0.; | 68 lp_df[i] = 0.; |
71 } | 69 } |
72 | 70 |
73 inp1 = 0.; inp2 = 0.; | 71 inp1 = 0.; inp2 = 0.; |
74 out1 = 0.; out2 = 0.; | 72 out1 = 0.; out2 = 0.; |
75 | 73 |
76 // backwards filetering on time-reversed df | 74 // backwards filetering on time-reversed df |
77 for (unsigned int i = 0;i < df.size();i++) | 75 for (unsigned int i = 0;i < df.size();i++) { |
78 { | |
79 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2; | 76 lp_df[i] = b[0]*df[i] + b[1]*inp1 + b[2]*inp2 - a[1]*out1 - a[2]*out2; |
80 inp2 = inp1; | 77 inp2 = inp1; |
81 inp1 = df[i]; | 78 inp1 = df[i]; |
82 out2 = out1; | 79 out2 = out1; |
83 out1 = lp_df[i]; | 80 out1 = lp_df[i]; |
84 } | 81 } |
85 | 82 |
86 // write the re-reversed (i.e. forward) version back to df | 83 // write the re-reversed (i.e. forward) version back to df |
87 for (unsigned int i = 0;i < df.size();i++) | 84 for (unsigned int i = 0;i < df.size();i++) { |
88 { | |
89 df[i] = lp_df[df.size()-i-1]; | 85 df[i] = lp_df[df.size()-i-1]; |
90 } | 86 } |
91 } | 87 } |
92 | 88 |
93 | 89 |
117 // accordingly. | 113 // accordingly. |
118 // note: 60*44100/512 is a magic number | 114 // note: 60*44100/512 is a magic number |
119 // this might (will?) break if a user specifies a different frame rate for the onset detection function | 115 // this might (will?) break if a user specifies a different frame rate for the onset detection function |
120 double rayparam = (60*44100/512)/inputtempo; | 116 double rayparam = (60*44100/512)/inputtempo; |
121 | 117 |
122 // these debug statements can be removed. | |
123 // std::cerr << "inputtempo" << inputtempo << std::endl; | |
124 // std::cerr << "rayparam" << rayparam << std::endl; | |
125 // std::cerr << "constraintempo" << constraintempo << std::endl; | |
126 | |
127 // make rayleigh weighting curve | 118 // make rayleigh weighting curve |
128 d_vec_t wv(wv_len); | 119 d_vec_t wv(wv_len); |
129 | 120 |
130 // check whether or not to use rayleigh weighting (if constraintempo is false) | 121 // check whether or not to use rayleigh weighting (if constraintempo is false) |
131 // or use gaussian weighting it (constraintempo is true) | 122 // or use gaussian weighting it (constraintempo is true) |
132 if (constraintempo) | 123 if (constraintempo) { |
133 { | 124 for (unsigned int i=0; i<wv.size(); i++) { |
134 for (unsigned int i=0; i<wv.size(); i++) | |
135 { | |
136 // MEPD 28/11/12 | 125 // MEPD 28/11/12 |
137 // do a gaussian weighting instead of rayleigh | 126 // do a gaussian weighting instead of rayleigh |
138 wv[i] = exp( (-1.*pow((static_cast<double> (i)-rayparam),2.)) / (2.*pow(rayparam/4.,2.)) ); | 127 wv[i] = exp( (-1.*pow((static_cast<double> (i)-rayparam),2.)) / (2.*pow(rayparam/4.,2.)) ); |
139 } | 128 } |
140 } | 129 } else { |
141 else | 130 for (unsigned int i=0; i<wv.size(); i++) { |
142 { | |
143 for (unsigned int i=0; i<wv.size(); i++) | |
144 { | |
145 // MEPD 28/11/12 | 131 // MEPD 28/11/12 |
146 // standard rayleigh weighting over periodicities | 132 // standard rayleigh weighting over periodicities |
147 wv[i] = (static_cast<double> (i) / pow(rayparam,2.)) * exp((-1.*pow(-static_cast<double> (i),2.)) / (2.*pow(rayparam,2.))); | 133 wv[i] = (static_cast<double> (i) / pow(rayparam,2.)) * exp((-1.*pow(-static_cast<double> (i),2.)) / (2.*pow(rayparam,2.))); |
148 } | 134 } |
149 } | 135 } |
155 // matrix to store output of comb filter bank, increment column of matrix at each frame | 141 // matrix to store output of comb filter bank, increment column of matrix at each frame |
156 d_mat_t rcfmat; | 142 d_mat_t rcfmat; |
157 int col_counter = -1; | 143 int col_counter = -1; |
158 | 144 |
159 // main loop for beat period calculation | 145 // main loop for beat period calculation |
160 for (unsigned int i=0; i+winlen<df.size(); i+=step) | 146 for (unsigned int i=0; i+winlen<df.size(); i+=step) { |
161 { | 147 |
162 // get dfframe | 148 // get dfframe |
163 d_vec_t dfframe(winlen); | 149 d_vec_t dfframe(winlen); |
164 for (unsigned int k=0; k<winlen; k++) | 150 for (unsigned int k=0; k<winlen; k++) { |
165 { | |
166 dfframe[k] = df[i+k]; | 151 dfframe[k] = df[i+k]; |
167 } | 152 } |
168 // get rcf vector for current frame | 153 // get rcf vector for current frame |
169 d_vec_t rcf(wv_len); | 154 d_vec_t rcf(wv_len); |
170 get_rcf(dfframe,wv,rcf); | 155 get_rcf(dfframe,wv,rcf); |
171 | 156 |
172 rcfmat.push_back( d_vec_t() ); // adds a new column | 157 rcfmat.push_back( d_vec_t() ); // adds a new column |
173 col_counter++; | 158 col_counter++; |
174 for (unsigned int j=0; j<rcf.size(); j++) | 159 for (unsigned int j=0; j<rcf.size(); j++) { |
175 { | |
176 rcfmat[col_counter].push_back( rcf[j] ); | 160 rcfmat[col_counter].push_back( rcf[j] ); |
177 } | 161 } |
178 } | 162 } |
179 | 163 |
180 // now call viterbi decoding function | 164 // now call viterbi decoding function |
195 | 179 |
196 MathUtilities::adaptiveThreshold(dfframe); | 180 MathUtilities::adaptiveThreshold(dfframe); |
197 | 181 |
198 d_vec_t acf(dfframe.size()); | 182 d_vec_t acf(dfframe.size()); |
199 | 183 |
200 | 184 for (unsigned int lag=0; lag<dfframe.size(); lag++) { |
201 for (unsigned int lag=0; lag<dfframe.size(); lag++) | |
202 { | |
203 double sum = 0.; | 185 double sum = 0.; |
204 double tmp = 0.; | 186 double tmp = 0.; |
205 | 187 |
206 for (unsigned int n=0; n<(dfframe.size()-lag); n++) | 188 for (unsigned int n=0; n<(dfframe.size()-lag); n++) { |
207 { | |
208 tmp = dfframe[n] * dfframe[n+lag]; | 189 tmp = dfframe[n] * dfframe[n+lag]; |
209 sum += tmp; | 190 sum += tmp; |
210 } | 191 } |
211 acf[lag] = static_cast<double> (sum/ (dfframe.size()-lag)); | 192 acf[lag] = static_cast<double> (sum/ (dfframe.size()-lag)); |
212 } | 193 } |
213 | 194 |
214 // now apply comb filtering | 195 // now apply comb filtering |
215 int numelem = 4; | 196 int numelem = 4; |
216 | 197 |
217 for (unsigned int i = 2;i < rcf.size();i++) // max beat period | 198 for (unsigned int i = 2;i < rcf.size();i++) { // max beat period |
218 { | 199 for (int a = 1;a <= numelem;a++) { // number of comb elements |
219 for (int a = 1;a <= numelem;a++) // number of comb elements | 200 for (int b = 1-a;b <= a-1;b++) { // general state using normalisation of comb elements |
220 { | 201 rcf[i-1] += ( acf[(a*i+b)-1]*wv[i-1] ) / (2.*a-1.); // calculate value for comb filter row |
221 for (int b = 1-a;b <= a-1;b++) // general state using normalisation of comb elements | |
222 { | |
223 rcf[i-1] += ( acf[(a*i+b)-1]*wv[i-1] ) / (2.*a-1.); // calculate value for comb filter row | |
224 } | 202 } |
225 } | 203 } |
226 } | 204 } |
227 | 205 |
228 // apply adaptive threshold to rcf | 206 // apply adaptive threshold to rcf |
229 MathUtilities::adaptiveThreshold(rcf); | 207 MathUtilities::adaptiveThreshold(rcf); |
230 | 208 |
231 double rcfsum =0.; | 209 double rcfsum =0.; |
232 for (unsigned int i=0; i<rcf.size(); i++) | 210 for (unsigned int i=0; i<rcf.size(); i++) { |
233 { | |
234 rcf[i] += EPS ; | 211 rcf[i] += EPS ; |
235 rcfsum += rcf[i]; | 212 rcfsum += rcf[i]; |
236 } | 213 } |
237 | 214 |
238 // normalise rcf to sum to unity | 215 // normalise rcf to sum to unity |
239 for (unsigned int i=0; i<rcf.size(); i++) | 216 for (unsigned int i=0; i<rcf.size(); i++) { |
240 { | |
241 rcf[i] /= (rcfsum + EPS); | 217 rcf[i] /= (rcfsum + EPS); |
242 } | 218 } |
243 } | 219 } |
244 | 220 |
245 void | 221 void |
248 // following Kevin Murphy's Viterbi decoding to get best path of | 224 // following Kevin Murphy's Viterbi decoding to get best path of |
249 // beat periods through rfcmat | 225 // beat periods through rfcmat |
250 | 226 |
251 // make transition matrix | 227 // make transition matrix |
252 d_mat_t tmat; | 228 d_mat_t tmat; |
253 for (unsigned int i=0;i<wv.size();i++) | 229 for (unsigned int i=0;i<wv.size();i++) { |
254 { | |
255 tmat.push_back ( d_vec_t() ); // adds a new column | 230 tmat.push_back ( d_vec_t() ); // adds a new column |
256 for (unsigned int j=0; j<wv.size(); j++) | 231 for (unsigned int j=0; j<wv.size(); j++) { |
257 { | |
258 tmat[i].push_back(0.); // fill with zeros initially | 232 tmat[i].push_back(0.); // fill with zeros initially |
259 } | 233 } |
260 } | 234 } |
261 | 235 |
262 // variance of Gaussians in transition matrix | 236 // variance of Gaussians in transition matrix |
263 // formed of Gaussians on diagonal - implies slow tempo change | 237 // formed of Gaussians on diagonal - implies slow tempo change |
264 double sigma = 8.; | 238 double sigma = 8.; |
265 // don't want really short beat periods, or really long ones | 239 // don't want really short beat periods, or really long ones |
266 for (unsigned int i=20;i <wv.size()-20; i++) | 240 for (unsigned int i=20;i <wv.size()-20; i++) { |
267 { | 241 for (unsigned int j=20; j<wv.size()-20; j++) { |
268 for (unsigned int j=20; j<wv.size()-20; j++) | |
269 { | |
270 double mu = static_cast<double>(i); | 242 double mu = static_cast<double>(i); |
271 tmat[i][j] = exp( (-1.*pow((j-mu),2.)) / (2.*pow(sigma,2.)) ); | 243 tmat[i][j] = exp( (-1.*pow((j-mu),2.)) / (2.*pow(sigma,2.)) ); |
272 } | 244 } |
273 } | 245 } |
274 | 246 |
275 // parameters for Viterbi decoding... this part is taken from | 247 // parameters for Viterbi decoding... this part is taken from |
276 // Murphy's matlab | 248 // Murphy's matlab |
277 | 249 |
278 d_mat_t delta; | 250 d_mat_t delta; |
279 i_mat_t psi; | 251 i_mat_t psi; |
280 for (unsigned int i=0;i <rcfmat.size(); i++) | 252 for (unsigned int i=0;i <rcfmat.size(); i++) { |
281 { | |
282 delta.push_back( d_vec_t()); | 253 delta.push_back( d_vec_t()); |
283 psi.push_back( i_vec_t()); | 254 psi.push_back( i_vec_t()); |
284 for (unsigned int j=0; j<rcfmat[i].size(); j++) | 255 for (unsigned int j=0; j<rcfmat[i].size(); j++) { |
285 { | |
286 delta[i].push_back(0.); // fill with zeros initially | 256 delta[i].push_back(0.); // fill with zeros initially |
287 psi[i].push_back(0); // fill with zeros initially | 257 psi[i].push_back(0); // fill with zeros initially |
288 } | 258 } |
289 } | 259 } |
290 | 260 |
291 | |
292 unsigned int T = delta.size(); | 261 unsigned int T = delta.size(); |
293 | 262 |
294 if (T < 2) return; // can't do anything at all meaningful | 263 if (T < 2) return; // can't do anything at all meaningful |
295 | 264 |
296 unsigned int Q = delta[0].size(); | 265 unsigned int Q = delta[0].size(); |
297 | 266 |
298 // initialize first column of delta | 267 // initialize first column of delta |
299 for (unsigned int j=0; j<Q; j++) | 268 for (unsigned int j=0; j<Q; j++) { |
300 { | |
301 delta[0][j] = wv[j] * rcfmat[0][j]; | 269 delta[0][j] = wv[j] * rcfmat[0][j]; |
302 psi[0][j] = 0; | 270 psi[0][j] = 0; |
303 } | 271 } |
304 | 272 |
305 double deltasum = 0.; | 273 double deltasum = 0.; |
306 for (unsigned int i=0; i<Q; i++) | 274 for (unsigned int i=0; i<Q; i++) { |
307 { | |
308 deltasum += delta[0][i]; | 275 deltasum += delta[0][i]; |
309 } | 276 } |
310 for (unsigned int i=0; i<Q; i++) | 277 for (unsigned int i=0; i<Q; i++) { |
311 { | |
312 delta[0][i] /= (deltasum + EPS); | 278 delta[0][i] /= (deltasum + EPS); |
313 } | 279 } |
314 | |
315 | 280 |
316 for (unsigned int t=1; t<T; t++) | 281 for (unsigned int t=1; t<T; t++) |
317 { | 282 { |
318 d_vec_t tmp_vec(Q); | 283 d_vec_t tmp_vec(Q); |
319 | 284 |
320 for (unsigned int j=0; j<Q; j++) | 285 for (unsigned int j=0; j<Q; j++) { |
321 { | 286 for (unsigned int i=0; i<Q; i++) { |
322 for (unsigned int i=0; i<Q; i++) | |
323 { | |
324 tmp_vec[i] = delta[t-1][i] * tmat[j][i]; | 287 tmp_vec[i] = delta[t-1][i] * tmat[j][i]; |
325 } | 288 } |
326 | 289 |
327 delta[t][j] = get_max_val(tmp_vec); | 290 delta[t][j] = get_max_val(tmp_vec); |
328 | 291 |
331 delta[t][j] *= rcfmat[t][j]; | 294 delta[t][j] *= rcfmat[t][j]; |
332 } | 295 } |
333 | 296 |
334 // normalise current delta column | 297 // normalise current delta column |
335 double deltasum = 0.; | 298 double deltasum = 0.; |
336 for (unsigned int i=0; i<Q; i++) | 299 for (unsigned int i=0; i<Q; i++) { |
337 { | |
338 deltasum += delta[t][i]; | 300 deltasum += delta[t][i]; |
339 } | 301 } |
340 for (unsigned int i=0; i<Q; i++) | 302 for (unsigned int i=0; i<Q; i++) { |
341 { | |
342 delta[t][i] /= (deltasum + EPS); | 303 delta[t][i] /= (deltasum + EPS); |
343 } | 304 } |
344 } | 305 } |
345 | 306 |
346 i_vec_t bestpath(T); | 307 i_vec_t bestpath(T); |
347 d_vec_t tmp_vec(Q); | 308 d_vec_t tmp_vec(Q); |
348 for (unsigned int i=0; i<Q; i++) | 309 for (unsigned int i=0; i<Q; i++) { |
349 { | |
350 tmp_vec[i] = delta[T-1][i]; | 310 tmp_vec[i] = delta[T-1][i]; |
351 } | 311 } |
352 | 312 |
353 // find starting point - best beat period for "last" frame | 313 // find starting point - best beat period for "last" frame |
354 bestpath[T-1] = get_max_ind(tmp_vec); | 314 bestpath[T-1] = get_max_ind(tmp_vec); |
355 | 315 |
356 // backtrace through index of maximum values in psi | 316 // backtrace through index of maximum values in psi |
357 for (unsigned int t=T-2; t>0 ;t--) | 317 for (unsigned int t=T-2; t>0 ;t--) { |
358 { | |
359 bestpath[t] = psi[t+1][bestpath[t+1]]; | 318 bestpath[t] = psi[t+1][bestpath[t+1]]; |
360 } | 319 } |
361 | 320 |
362 // weird but necessary hack -- couldn't get above loop to terminate at t >= 0 | 321 // weird but necessary hack -- couldn't get above loop to terminate at t >= 0 |
363 bestpath[0] = psi[1][bestpath[1]]; | 322 bestpath[0] = psi[1][bestpath[1]]; |
364 | 323 |
365 unsigned int lastind = 0; | 324 unsigned int lastind = 0; |
366 for (unsigned int i=0; i<T; i++) | 325 for (unsigned int i=0; i<T; i++) { |
367 { | |
368 unsigned int step = 128; | 326 unsigned int step = 128; |
369 for (unsigned int j=0; j<step; j++) | 327 for (unsigned int j=0; j<step; j++) { |
370 { | |
371 lastind = i*step+j; | 328 lastind = i*step+j; |
372 beat_period[lastind] = bestpath[i]; | 329 beat_period[lastind] = bestpath[i]; |
373 } | 330 } |
374 // std::cerr << "bestpath[" << i << "] = " << bestpath[i] << " (used for beat_periods " << i*step << " to " << i*step+step-1 << ")" << std::endl; | 331 // std::cerr << "bestpath[" << i << "] = " << bestpath[i] << " (used for beat_periods " << i*step << " to " << i*step+step-1 << ")" << std::endl; |
375 } | 332 } |
376 | 333 |
377 //fill in the last values... | 334 //fill in the last values... |
378 for (unsigned int i=lastind; i<beat_period.size(); i++) | 335 for (unsigned int i=lastind; i<beat_period.size(); i++) { |
379 { | |
380 beat_period[i] = beat_period[lastind]; | 336 beat_period[i] = beat_period[lastind]; |
381 } | 337 } |
382 | 338 |
383 for (unsigned int i = 0; i < beat_period.size(); i++) | 339 for (unsigned int i = 0; i < beat_period.size(); i++) { |
384 { | |
385 tempi.push_back((60. * m_rate / m_increment)/beat_period[i]); | 340 tempi.push_back((60. * m_rate / m_increment)/beat_period[i]); |
386 } | 341 } |
387 } | 342 } |
388 | 343 |
389 double | 344 double |
390 TempoTrackV2::get_max_val(const d_vec_t &df) | 345 TempoTrackV2::get_max_val(const d_vec_t &df) |
391 { | 346 { |
392 double maxval = 0.; | 347 double maxval = 0.; |
393 for (unsigned int i=0; i<df.size(); i++) | 348 for (unsigned int i=0; i<df.size(); i++) { |
394 { | 349 if (maxval < df[i]) { |
395 if (maxval < df[i]) | |
396 { | |
397 maxval = df[i]; | 350 maxval = df[i]; |
398 } | 351 } |
399 } | 352 } |
400 | 353 |
401 return maxval; | 354 return maxval; |
404 int | 357 int |
405 TempoTrackV2::get_max_ind(const d_vec_t &df) | 358 TempoTrackV2::get_max_ind(const d_vec_t &df) |
406 { | 359 { |
407 double maxval = 0.; | 360 double maxval = 0.; |
408 int ind = 0; | 361 int ind = 0; |
409 for (unsigned int i=0; i<df.size(); i++) | 362 for (unsigned int i=0; i<df.size(); i++) { |
410 { | 363 if (maxval < df[i]) { |
411 if (maxval < df[i]) | |
412 { | |
413 maxval = df[i]; | 364 maxval = df[i]; |
414 ind = i; | 365 ind = i; |
415 } | 366 } |
416 } | 367 } |
417 | 368 |
420 | 371 |
421 void | 372 void |
422 TempoTrackV2::normalise_vec(d_vec_t &df) | 373 TempoTrackV2::normalise_vec(d_vec_t &df) |
423 { | 374 { |
424 double sum = 0.; | 375 double sum = 0.; |
425 for (unsigned int i=0; i<df.size(); i++) | 376 for (unsigned int i=0; i<df.size(); i++) { |
426 { | |
427 sum += df[i]; | 377 sum += df[i]; |
428 } | 378 } |
429 | 379 |
430 for (unsigned int i=0; i<df.size(); i++) | 380 for (unsigned int i=0; i<df.size(); i++) { |
431 { | |
432 df[i]/= (sum + EPS); | 381 df[i]/= (sum + EPS); |
433 } | 382 } |
434 } | 383 } |
435 | 384 |
436 // MEPD 28/11/12 | 385 // MEPD 28/11/12 |
446 | 395 |
447 d_vec_t cumscore(df.size()); // store cumulative score | 396 d_vec_t cumscore(df.size()); // store cumulative score |
448 i_vec_t backlink(df.size()); // backlink (stores best beat locations at each time instant) | 397 i_vec_t backlink(df.size()); // backlink (stores best beat locations at each time instant) |
449 d_vec_t localscore(df.size()); // localscore, for now this is the same as the detection function | 398 d_vec_t localscore(df.size()); // localscore, for now this is the same as the detection function |
450 | 399 |
451 for (unsigned int i=0; i<df.size(); i++) | 400 for (unsigned int i=0; i<df.size(); i++) { |
452 { | |
453 localscore[i] = df[i]; | 401 localscore[i] = df[i]; |
454 backlink[i] = -1; | 402 backlink[i] = -1; |
455 } | 403 } |
456 | 404 |
457 //double tightness = 4.; | 405 //double tightness = 4.; |
460 // debug statements that can be removed. | 408 // debug statements that can be removed. |
461 // std::cerr << "alpha" << alpha << std::endl; | 409 // std::cerr << "alpha" << alpha << std::endl; |
462 // std::cerr << "tightness" << tightness << std::endl; | 410 // std::cerr << "tightness" << tightness << std::endl; |
463 | 411 |
464 // main loop | 412 // main loop |
465 for (unsigned int i=0; i<localscore.size(); i++) | 413 for (unsigned int i=0; i<localscore.size(); i++) { |
466 { | 414 |
467 int prange_min = -2*beat_period[i]; | 415 int prange_min = -2*beat_period[i]; |
468 int prange_max = round(-0.5*beat_period[i]); | 416 int prange_max = round(-0.5*beat_period[i]); |
469 | 417 |
470 // transition range | 418 // transition range |
471 d_vec_t txwt (prange_max - prange_min + 1); | 419 d_vec_t txwt (prange_max - prange_min + 1); |
472 d_vec_t scorecands (txwt.size()); | 420 d_vec_t scorecands (txwt.size()); |
473 | 421 |
474 for (unsigned int j=0;j<txwt.size();j++) | 422 for (unsigned int j=0;j<txwt.size();j++) { |
475 { | 423 |
476 double mu = static_cast<double> (beat_period[i]); | 424 double mu = static_cast<double> (beat_period[i]); |
477 txwt[j] = exp( -0.5*pow(tightness * log((round(2*mu)-j)/mu),2)); | 425 txwt[j] = exp( -0.5*pow(tightness * log((round(2*mu)-j)/mu),2)); |
478 | 426 |
479 // IF IN THE ALLOWED RANGE, THEN LOOK AT CUMSCORE[I+PRANGE_MIN+J | 427 // IF IN THE ALLOWED RANGE, THEN LOOK AT CUMSCORE[I+PRANGE_MIN+J |
480 // ELSE LEAVE AT DEFAULT VALUE FROM INITIALISATION: D_VEC_T SCORECANDS (TXWT.SIZE()); | 428 // ELSE LEAVE AT DEFAULT VALUE FROM INITIALISATION: D_VEC_T SCORECANDS (TXWT.SIZE()); |
481 | 429 |
482 int cscore_ind = i+prange_min+j; | 430 int cscore_ind = i+prange_min+j; |
483 if (cscore_ind >= 0) | 431 if (cscore_ind >= 0) { |
484 { | |
485 scorecands[j] = txwt[j] * cumscore[cscore_ind]; | 432 scorecands[j] = txwt[j] * cumscore[cscore_ind]; |
486 } | 433 } |
487 } | 434 } |
488 | 435 |
489 // find max value and index of maximum value | 436 // find max value and index of maximum value |
496 // std::cerr << "backlink[" << i << "] <= " << backlink[i] << std::endl; | 443 // std::cerr << "backlink[" << i << "] <= " << backlink[i] << std::endl; |
497 } | 444 } |
498 | 445 |
499 // STARTING POINT, I.E. LAST BEAT.. PICK A STRONG POINT IN cumscore VECTOR | 446 // STARTING POINT, I.E. LAST BEAT.. PICK A STRONG POINT IN cumscore VECTOR |
500 d_vec_t tmp_vec; | 447 d_vec_t tmp_vec; |
501 for (unsigned int i=cumscore.size() - beat_period[beat_period.size()-1] ; i<cumscore.size(); i++) | 448 for (unsigned int i=cumscore.size() - beat_period[beat_period.size()-1] ; i<cumscore.size(); i++) { |
502 { | |
503 tmp_vec.push_back(cumscore[i]); | 449 tmp_vec.push_back(cumscore[i]); |
504 } | 450 } |
505 | 451 |
506 int startpoint = get_max_ind(tmp_vec) + cumscore.size() - beat_period[beat_period.size()-1] ; | 452 int startpoint = get_max_ind(tmp_vec) + |
453 cumscore.size() - beat_period[beat_period.size()-1] ; | |
507 | 454 |
508 // can happen if no results obtained earlier (e.g. input too short) | 455 // can happen if no results obtained earlier (e.g. input too short) |
509 if (startpoint >= (int)backlink.size()) startpoint = backlink.size()-1; | 456 if (startpoint >= (int)backlink.size()) { |
457 startpoint = backlink.size()-1; | |
458 } | |
510 | 459 |
511 // USE BACKLINK TO GET EACH NEW BEAT (TOWARDS THE BEGINNING OF THE FILE) | 460 // USE BACKLINK TO GET EACH NEW BEAT (TOWARDS THE BEGINNING OF THE FILE) |
512 // BACKTRACKING FROM THE END TO THE BEGINNING.. MAKING SURE NOT TO GO BEFORE SAMPLE 0 | 461 // BACKTRACKING FROM THE END TO THE BEGINNING.. MAKING SURE NOT TO GO BEFORE SAMPLE 0 |
513 i_vec_t ibeats; | 462 i_vec_t ibeats; |
514 ibeats.push_back(startpoint); | 463 ibeats.push_back(startpoint); |
515 // std::cerr << "startpoint = " << startpoint << std::endl; | 464 // std::cerr << "startpoint = " << startpoint << std::endl; |
516 while (backlink[ibeats.back()] > 0) | 465 while (backlink[ibeats.back()] > 0) { |
517 { | |
518 // std::cerr << "backlink[" << ibeats.back() << "] = " << backlink[ibeats.back()] << std::endl; | 466 // std::cerr << "backlink[" << ibeats.back() << "] = " << backlink[ibeats.back()] << std::endl; |
519 int b = ibeats.back(); | 467 int b = ibeats.back(); |
520 if (backlink[b] == b) break; // shouldn't happen... haha | 468 if (backlink[b] == b) break; // shouldn't happen... haha |
521 ibeats.push_back(backlink[b]); | 469 ibeats.push_back(backlink[b]); |
522 } | 470 } |
523 | 471 |
524 // REVERSE SEQUENCE OF IBEATS AND STORE AS BEATS | 472 // REVERSE SEQUENCE OF IBEATS AND STORE AS BEATS |
525 for (unsigned int i=0; i<ibeats.size(); i++) | 473 for (unsigned int i=0; i<ibeats.size(); i++) { |
526 { | |
527 beats.push_back( static_cast<double>(ibeats[ibeats.size()-i-1]) ); | 474 beats.push_back( static_cast<double>(ibeats[ibeats.size()-i-1]) ); |
528 } | 475 } |
529 } | 476 } |
530 | 477 |
531 | 478 |