comparison src/OnsetDetectionFunction.cpp @ 72:f4d9410f187e

flow: Merged <release> '1.0.0' to <master> ('master').
author Adam Stark <adamstark@users.noreply.github.com>
date Tue, 08 Jul 2014 12:32:27 +0100
parents b387d8327729
children 5eeabb24d677
comparison
equal deleted inserted replaced
44:3049937d6ef1 72:f4d9410f187e
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 */ 19 */
20 //======================================================================= 20 //=======================================================================
21 21
22 #include <math.h> 22 #include <math.h>
23 #include <iostream>
24 #include "OnsetDetectionFunction.h" 23 #include "OnsetDetectionFunction.h"
25 using namespace std; 24
26 25
27 //------------------------------------------------------------------------------- 26 //=======================================================================
28 // Constructor 27 OnsetDetectionFunction::OnsetDetectionFunction(int hopSize_,int frameSize_) : onsetDetectionFunctionType(ComplexSpectralDifferenceHWR), windowType(HanningWindow)
29 OnsetDetectionFunction :: OnsetDetectionFunction() 28 {
29 // indicate that we have not initialised yet
30 initialised = false;
31
32 // set pi
33 pi = 3.14159265358979;
34
35 // initialise with arguments to constructor
36 initialise(hopSize_,frameSize_,ComplexSpectralDifferenceHWR,HanningWindow);
37 }
38
39 //=======================================================================
40 OnsetDetectionFunction::OnsetDetectionFunction(int hopSize_,int frameSize_,int onsetDetectionFunctionType_,int windowType_) : onsetDetectionFunctionType(ComplexSpectralDifferenceHWR), windowType(HanningWindow)
30 { 41 {
31 // indicate that we have not initialised yet 42 // indicate that we have not initialised yet
32 initialised = 0; 43 initialised = false;
33 44
34 // set pi 45 // set pi
35 pi = 3.14159265358979; 46 pi = 3.14159265358979;
36 47
37 // initialise with hopsize = 512, framesize = 1024, complex spectral difference DF and hanning window
38 initialise(512,1024,6,1);
39 }
40
41
42 //-------------------------------------------------------------------------------
43 // Constructor (with arguments)
44 OnsetDetectionFunction :: OnsetDetectionFunction(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
45 {
46 // indicate that we have not initialised yet
47 initialised = 0;
48
49 // set pi
50 pi = 3.14159265358979;
51
52 // initialise with arguments to constructor 48 // initialise with arguments to constructor
53 initialise(arg_hsize,arg_fsize,arg_df_type,arg_win_type); 49 initialise(hopSize_,frameSize_,onsetDetectionFunctionType_,windowType_);
54 } 50 }
55 51
56 52
57 //-------------------------------------------------------------------------------------- 53 //=======================================================================
58 // Destructor 54 OnsetDetectionFunction::~OnsetDetectionFunction()
59 OnsetDetectionFunction :: ~OnsetDetectionFunction() 55 {
60 { 56 if (initialised)
61 // destroy fft plan 57 {
62 fftw_destroy_plan(p); 58 // destroy fft plan
63 fftw_free(in); 59 fftw_destroy_plan(p);
64 fftw_free(out); 60 fftw_free(complexIn);
65 61 fftw_free(complexOut);
66 // deallocate memory 62 }
67 delete [] frame; 63 }
68 frame = NULL; 64
69 delete [] window; 65 //=======================================================================
70 window = NULL; 66 void OnsetDetectionFunction::initialise(int hopSize_,int frameSize_)
71 delete [] wframe; 67 {
72 wframe = NULL; 68 // use the already initialised onset detection function and window type and
73 delete [] mag; 69 // pass the new frame and hop size to the main initialisation function
74 mag = NULL; 70 initialise(hopSize_, frameSize_, onsetDetectionFunctionType, windowType);
75 delete [] mag_old; 71 }
76 mag_old = NULL; 72
77 delete [] phase; 73 //=======================================================================
78 phase = NULL; 74 void OnsetDetectionFunction::initialise(int hopSize_,int frameSize_,int onsetDetectionFunctionType_,int windowType_)
79 delete [] phase_old; 75 {
80 phase_old = NULL; 76 if (initialised) // if we have already initialised FFT plan
81 delete [] phase_old_2; 77 {
82 phase_old_2 = NULL;
83 }
84
85 //-------------------------------------------------------------------------------
86 // Initialisation
87 void OnsetDetectionFunction :: initialise(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
88 {
89 if (initialised == 1) // if we have already initialised some buffers and an FFT plan
90 {
91 //////////////////////////////////
92 // TIDY UP FIRST - If initialise is called after the class has been initialised
93 // then we want to free up memory and cancel existing FFT plans
94
95 // destroy fft plan 78 // destroy fft plan
96 fftw_destroy_plan(p); 79 fftw_destroy_plan(p);
97 fftw_free(in); 80 fftw_free(complexIn);
98 fftw_free(out); 81 fftw_free(complexOut);
99 82
100 83 }
101 // deallocate memory 84
102 delete [] frame; 85 hopSize = hopSize_; // set hopsize
103 frame = NULL; 86 frameSize = frameSize_; // set framesize
104 delete [] window; 87
105 window = NULL; 88 onsetDetectionFunctionType = onsetDetectionFunctionType_; // set detection function type
106 delete [] wframe; 89 windowType = windowType_; // set window type
107 wframe = NULL;
108 delete [] mag;
109 mag = NULL;
110 delete [] mag_old;
111 mag_old = NULL;
112 delete [] phase;
113 phase = NULL;
114 delete [] phase_old;
115 phase_old = NULL;
116 delete [] phase_old_2;
117 phase_old_2 = NULL;
118
119 ////// END TIDY UP ///////////////
120 //////////////////////////////////
121 }
122
123 hopsize = arg_hsize; // set hopsize
124 framesize = arg_fsize; // set framesize
125
126 df_type = arg_df_type; // set detection function type
127 90
128 // initialise buffers 91 // initialise buffers
129 frame = new double[framesize]; 92 frame.resize(frameSize);
130 window = new double[framesize]; 93 window.resize(frameSize);
131 wframe = new double[framesize]; 94 magSpec.resize(frameSize);
132 95 prevMagSpec.resize(frameSize);
133 mag = new double[framesize]; 96 phase.resize(frameSize);
134 mag_old = new double[framesize]; 97 prevPhase.resize(frameSize);
135 98 prevPhase2.resize(frameSize);
136 phase = new double[framesize];
137 phase_old = new double[framesize];
138 phase_old_2 = new double[framesize];
139 99
140 100
141 // set the window to the specified type 101 // set the window to the specified type
142 switch (arg_win_type){ 102 switch (windowType){
143 case 0: 103 case RectangularWindow:
144 set_win_rectangular(); // Rectangular window 104 calculateRectangularWindow(); // Rectangular window
145 break; 105 break;
146 case 1: 106 case HanningWindow:
147 set_win_hanning(); // Hanning Window 107 calculateHanningWindow(); // Hanning Window
148 break; 108 break;
149 case 2: 109 case HammingWindow:
150 set_win_hamming(); // Hamming Window 110 calclulateHammingWindow(); // Hamming Window
151 break; 111 break;
152 case 3: 112 case BlackmanWindow:
153 set_win_blackman(); // Blackman Window 113 calculateBlackmanWindow(); // Blackman Window
154 break; 114 break;
155 case 4: 115 case TukeyWindow:
156 set_win_tukey(); // Tukey Window 116 calculateTukeyWindow(); // Tukey Window
157 break; 117 break;
158 default: 118 default:
159 set_win_hanning(); // DEFAULT: Hanning Window 119 calculateHanningWindow(); // DEFAULT: Hanning Window
160 } 120 }
161
162
163
164 121
165 // initialise previous magnitude spectrum to zero 122 // initialise previous magnitude spectrum to zero
166 for (int i = 0;i < framesize;i++) 123 for (int i = 0;i < frameSize;i++)
167 { 124 {
168 mag_old[i] = 0.0; 125 prevMagSpec[i] = 0.0;
169 phase_old[i] = 0.0; 126 prevPhase[i] = 0.0;
170 phase_old_2[i] = 0.0; 127 prevPhase2[i] = 0.0;
171 frame[i] = 0.0; 128 frame[i] = 0.0;
172 } 129 }
173 130
174 energy_sum_old = 0.0; // initialise previous energy sum value to zero 131 prevEnergySum = 0.0; // initialise previous energy sum value to zero
175 132
176 /* Init fft */ 133 /* Init fft */
177 in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize); // complex array to hold fft data 134 complexIn = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * frameSize); // complex array to hold fft data
178 out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize); // complex array to hold fft data 135 complexOut = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * frameSize); // complex array to hold fft data
179 p = fftw_plan_dft_1d(framesize, in, out, FFTW_FORWARD, FFTW_ESTIMATE); // FFT plan initialisation 136 p = fftw_plan_dft_1d(frameSize, complexIn, complexOut, FFTW_FORWARD, FFTW_ESTIMATE); // FFT plan initialisation
180 137
181 initialised = 1; 138 initialised = true;
182 } 139 }
183 140
184 //-------------------------------------------------------------------------------------- 141 //=======================================================================
185 // set the detection function type to that specified by the argument 142 void OnsetDetectionFunction :: setOnsetDetectionFunctionType(int onsetDetectionFunctionType_)
186 void OnsetDetectionFunction :: set_df_type(int arg_df_type) 143 {
187 { 144 onsetDetectionFunctionType = onsetDetectionFunctionType_; // set detection function type
188 df_type = arg_df_type; // set detection function type 145 }
189 } 146
190 147 //=======================================================================
191 148 double OnsetDetectionFunction :: calculateOnsetDetectionFunctionSample(double *buffer)
192 //--------------------------------------------------------------------------------------
193 // calculates a single detection function sample from a single audio frame.
194 double OnsetDetectionFunction :: getDFsample(double inputbuffer[])
195 { 149 {
196 double df_sample; 150 double odfSample;
197 151
198 // shift audio samples back in frame by hop size 152 // shift audio samples back in frame by hop size
199 for (int i = 0; i < (framesize-hopsize);i++) 153 for (int i = 0; i < (frameSize-hopSize);i++)
200 { 154 {
201 frame[i] = frame[i+hopsize]; 155 frame[i] = frame[i+hopSize];
202 } 156 }
203 157
204 // add new samples to frame from input buffer 158 // add new samples to frame from input buffer
205 int j = 0; 159 int j = 0;
206 for (int i = (framesize-hopsize);i < framesize;i++) 160 for (int i = (frameSize-hopSize);i < frameSize;i++)
207 { 161 {
208 frame[i] = inputbuffer[j]; 162 frame[i] = buffer[j];
209 j++; 163 j++;
210 } 164 }
211 165
212 switch (df_type){ 166 switch (onsetDetectionFunctionType){
213 case 0: 167 case EnergyEnvelope:
214 df_sample = energy_envelope(); // calculate energy envelope detection function sample 168 {
215 break; 169 // calculate energy envelope detection function sample
216 case 1: 170 odfSample = energyEnvelope();
217 df_sample = energy_difference(); // calculate half-wave rectified energy difference detection function sample 171 break;
218 break; 172 }
219 case 2: 173 case EnergyDifference:
220 df_sample = spectral_difference(); // calculate spectral difference detection function sample 174 {
221 break; 175 // calculate half-wave rectified energy difference detection function sample
222 case 3: 176 odfSample = energyDifference();
223 df_sample = spectral_difference_hwr(); // calculate spectral difference detection function sample (half wave rectified) 177 break;
224 break; 178 }
225 case 4: 179 case SpectralDifference:
226 df_sample = phase_deviation(); // calculate phase deviation detection function sample (half wave rectified) 180 {
227 break; 181 // calculate spectral difference detection function sample
228 case 5: 182 odfSample = spectralDifference();
229 df_sample = complex_spectral_difference(); // calcualte complex spectral difference detection function sample 183 break;
230 break; 184 }
231 case 6: 185 case SpectralDifferenceHWR:
232 df_sample = complex_spectral_difference_hwr(); // calcualte complex spectral difference detection function sample (half-wave rectified) 186 {
233 break; 187 // calculate spectral difference detection function sample (half wave rectified)
234 case 7: 188 odfSample = spectralDifferenceHWR();
235 df_sample = high_frequency_content(); // calculate high frequency content detection function sample 189 break;
236 break; 190 }
237 case 8: 191 case PhaseDeviation:
238 df_sample = high_frequency_spectral_difference(); // calculate high frequency spectral difference detection function sample 192 {
239 break; 193 // calculate phase deviation detection function sample (half wave rectified)
240 case 9: 194 odfSample = phaseDeviation();
241 df_sample = high_frequency_spectral_difference_hwr(); // calculate high frequency spectral difference detection function (half-wave rectified) 195 break;
242 break; 196 }
197 case ComplexSpectralDifference:
198 {
199 // calcualte complex spectral difference detection function sample
200 odfSample = complexSpectralDifference();
201 break;
202 }
203 case ComplexSpectralDifferenceHWR:
204 {
205 // calcualte complex spectral difference detection function sample (half-wave rectified)
206 odfSample = complexSpectralDifferenceHWR();
207 break;
208 }
209 case HighFrequencyContent:
210 {
211 // calculate high frequency content detection function sample
212 odfSample = highFrequencyContent();
213 break;
214 }
215 case HighFrequencySpectralDifference:
216 {
217 // calculate high frequency spectral difference detection function sample
218 odfSample = highFrequencySpectralDifference();
219 break;
220 }
221 case HighFrequencySpectralDifferenceHWR:
222 {
223 // calculate high frequency spectral difference detection function (half-wave rectified)
224 odfSample = highFrequencySpectralDifferenceHWR();
225 break;
226 }
243 default: 227 default:
244 df_sample = 1.0; 228 {
245 } 229 odfSample = 1.0;
246 230 }
247 return df_sample; 231 }
248 } 232
249 233 return odfSample;
250 234 }
251 //-------------------------------------------------------------------------------------- 235
252 // performs the fft, storing the complex result in 'out' 236
253 void OnsetDetectionFunction :: perform_FFT() 237 //=======================================================================
254 { 238 void OnsetDetectionFunction :: performFFT()
255 int fsize2 = (framesize/2); 239 {
240 int fsize2 = (frameSize/2);
256 241
257 // window frame and copy to complex array, swapping the first and second half of the signal 242 // window frame and copy to complex array, swapping the first and second half of the signal
258 for (int i = 0;i < fsize2;i++) 243 for (int i = 0;i < fsize2;i++)
259 { 244 {
260 in[i][0] = frame[i+fsize2] * window[i+fsize2]; 245 complexIn[i][0] = frame[i+fsize2] * window[i+fsize2];
261 in[i][1] = 0.0; 246 complexIn[i][1] = 0.0;
262 in[i+fsize2][0] = frame[i] * window[i]; 247 complexIn[i+fsize2][0] = frame[i] * window[i];
263 in[i+fsize2][1] = 0.0; 248 complexIn[i+fsize2][1] = 0.0;
264 } 249 }
265 250
266 // perform the fft 251 // perform the fft
267 fftw_execute(p); 252 fftw_execute(p);
268 } 253 }
269 254
270 //////////////////////////////////////////////////////////////////////////////////////////////// 255 ////////////////////////////////////////////////////////////////////////////////////////////////
271 //////////////////////////////////////////////////////////////////////////////////////////////// 256 ////////////////////////////////////////////////////////////////////////////////////////////////
272 ////////////////////////////// Methods for Detection Functions ///////////////////////////////// 257 ////////////////////////////// Methods for Detection Functions /////////////////////////////////
273 258
274 //-------------------------------------------------------------------------------------- 259 //=======================================================================
275 // calculates an energy envelope detection function sample 260 double OnsetDetectionFunction :: energyEnvelope()
276 double OnsetDetectionFunction :: energy_envelope()
277 { 261 {
278 double sum; 262 double sum;
279 263
280 sum = 0; // initialise sum 264 sum = 0; // initialise sum
281 265
282 // sum the squares of the samples 266 // sum the squares of the samples
283 for (int i = 0;i < framesize;i++) 267 for (int i = 0;i < frameSize;i++)
284 { 268 {
285 sum = sum + (frame[i]*frame[i]); 269 sum = sum + (frame[i]*frame[i]);
286 } 270 }
287 271
288 return sum; // return sum 272 return sum; // return sum
289 } 273 }
290 274
291 //-------------------------------------------------------------------------------------- 275 //=======================================================================
292 // calculates a half-wave rectified energy difference detection function sample 276 double OnsetDetectionFunction :: energyDifference()
293 double OnsetDetectionFunction :: energy_difference()
294 { 277 {
295 double sum; 278 double sum;
296 double sample; 279 double sample;
297 280
298 sum = 0; // initialise sum 281 sum = 0; // initialise sum
299 282
300 // sum the squares of the samples 283 // sum the squares of the samples
301 for (int i = 0;i < framesize;i++) 284 for (int i = 0;i < frameSize;i++)
302 { 285 {
303 sum = sum + (frame[i]*frame[i]); 286 sum = sum + (frame[i]*frame[i]);
304 } 287 }
305 288
306 sample = sum - energy_sum_old; // sample is first order difference in energy 289 sample = sum - prevEnergySum; // sample is first order difference in energy
307 290
308 energy_sum_old = sum; // store energy value for next calculation 291 prevEnergySum = sum; // store energy value for next calculation
309 292
310 if (sample > 0) 293 if (sample > 0)
311 { 294 {
312 return sample; // return difference 295 return sample; // return difference
313 } 296 }
315 { 298 {
316 return 0; 299 return 0;
317 } 300 }
318 } 301 }
319 302
320 //-------------------------------------------------------------------------------------- 303 //=======================================================================
321 // calculates a spectral difference detection function sample 304 double OnsetDetectionFunction :: spectralDifference()
322 double OnsetDetectionFunction :: spectral_difference()
323 { 305 {
324 double diff; 306 double diff;
325 double sum; 307 double sum;
326 308
327 // perform the FFT 309 // perform the FFT
328 perform_FFT(); 310 performFFT();
329 311
330 // compute first (N/2)+1 mag values 312 // compute first (N/2)+1 mag values
331 for (int i = 0;i < (framesize/2)+1;i++) 313 for (int i = 0;i < (frameSize/2)+1;i++)
332 { 314 {
333 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 315 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
334 } 316 }
335 // mag spec symmetric above (N/2)+1 so copy previous values 317 // mag spec symmetric above (N/2)+1 so copy previous values
336 for (int i = (framesize/2)+1;i < framesize;i++) 318 for (int i = (frameSize/2)+1;i < frameSize;i++)
337 { 319 {
338 mag[i] = mag[framesize-i]; 320 magSpec[i] = magSpec[frameSize-i];
339 } 321 }
340 322
341 sum = 0; // initialise sum to zero 323 sum = 0; // initialise sum to zero
342 324
343 for (int i = 0;i < framesize;i++) 325 for (int i = 0;i < frameSize;i++)
344 { 326 {
345 // calculate difference 327 // calculate difference
346 diff = mag[i] - mag_old[i]; 328 diff = magSpec[i] - prevMagSpec[i];
347 329
348 // ensure all difference values are positive 330 // ensure all difference values are positive
349 if (diff < 0) 331 if (diff < 0)
350 { 332 {
351 diff = diff*-1; 333 diff = diff*-1;
353 335
354 // add difference to sum 336 // add difference to sum
355 sum = sum+diff; 337 sum = sum+diff;
356 338
357 // store magnitude spectrum bin for next detection function sample calculation 339 // store magnitude spectrum bin for next detection function sample calculation
358 mag_old[i] = mag[i]; 340 prevMagSpec[i] = magSpec[i];
359 } 341 }
360 342
361 return sum; 343 return sum;
362 } 344 }
363 345
364 //-------------------------------------------------------------------------------------- 346 //=======================================================================
365 // calculates a spectral difference detection function sample 347 double OnsetDetectionFunction :: spectralDifferenceHWR()
366 double OnsetDetectionFunction :: spectral_difference_hwr()
367 { 348 {
368 double diff; 349 double diff;
369 double sum; 350 double sum;
370 351
371 // perform the FFT 352 // perform the FFT
372 perform_FFT(); 353 performFFT();
373 354
374 // compute first (N/2)+1 mag values 355 // compute first (N/2)+1 mag values
375 for (int i = 0;i < (framesize/2)+1;i++) 356 for (int i = 0;i < (frameSize/2)+1;i++)
376 { 357 {
377 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 358 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
378 } 359 }
379 // mag spec symmetric above (N/2)+1 so copy previous values 360 // mag spec symmetric above (N/2)+1 so copy previous values
380 for (int i = (framesize/2)+1;i < framesize;i++) 361 for (int i = (frameSize/2)+1;i < frameSize;i++)
381 { 362 {
382 mag[i] = mag[framesize-i]; 363 magSpec[i] = magSpec[frameSize-i];
383 } 364 }
384 365
385 sum = 0; // initialise sum to zero 366 sum = 0; // initialise sum to zero
386 367
387 for (int i = 0;i < framesize;i++) 368 for (int i = 0;i < frameSize;i++)
388 { 369 {
389 // calculate difference 370 // calculate difference
390 diff = mag[i] - mag_old[i]; 371 diff = magSpec[i] - prevMagSpec[i];
391 372
392 // only add up positive differences 373 // only add up positive differences
393 if (diff > 0) 374 if (diff > 0)
394 { 375 {
395 // add difference to sum 376 // add difference to sum
397 } 378 }
398 379
399 380
400 381
401 // store magnitude spectrum bin for next detection function sample calculation 382 // store magnitude spectrum bin for next detection function sample calculation
402 mag_old[i] = mag[i]; 383 prevMagSpec[i] = magSpec[i];
403 } 384 }
404 385
405 return sum; 386 return sum;
406 } 387 }
407 388
408 389
409 //-------------------------------------------------------------------------------------- 390 //=======================================================================
410 // calculates a phase deviation detection function sample 391 double OnsetDetectionFunction :: phaseDeviation()
411 double OnsetDetectionFunction :: phase_deviation()
412 { 392 {
413 double dev,pdev; 393 double dev,pdev;
414 double sum; 394 double sum;
415 395
416 // perform the FFT 396 // perform the FFT
417 perform_FFT(); 397 performFFT();
418 398
419 sum = 0; // initialise sum to zero 399 sum = 0; // initialise sum to zero
420 400
421 // compute phase values from fft output and sum deviations 401 // compute phase values from fft output and sum deviations
422 for (int i = 0;i < framesize;i++) 402 for (int i = 0;i < frameSize;i++)
423 { 403 {
424 // calculate phase value 404 // calculate phase value
425 phase[i] = atan2(out[i][1],out[i][0]); 405 phase[i] = atan2(complexOut[i][1],complexOut[i][0]);
426 406
427 // calculate magnitude value 407 // calculate magnitude value
428 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 408 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
429 409
430 410
431 // if bin is not just a low energy bin then examine phase deviation 411 // if bin is not just a low energy bin then examine phase deviation
432 if (mag[i] > 0.1) 412 if (magSpec[i] > 0.1)
433 { 413 {
434 dev = phase[i] - (2*phase_old[i]) + phase_old_2[i]; // phase deviation 414 dev = phase[i] - (2*prevPhase[i]) + prevPhase2[i]; // phase deviation
435 pdev = princarg(dev); // wrap into [-pi,pi] range 415 pdev = princarg(dev); // wrap into [-pi,pi] range
436 416
437 // make all values positive 417 // make all values positive
438 if (pdev < 0) 418 if (pdev < 0)
439 { 419 {
443 // add to sum 423 // add to sum
444 sum = sum + pdev; 424 sum = sum + pdev;
445 } 425 }
446 426
447 // store values for next calculation 427 // store values for next calculation
448 phase_old_2[i] = phase_old[i]; 428 prevPhase2[i] = prevPhase[i];
449 phase_old[i] = phase[i]; 429 prevPhase[i] = phase[i];
450 } 430 }
451 431
452 return sum; 432 return sum;
453 } 433 }
454 434
455 //-------------------------------------------------------------------------------------- 435 //=======================================================================
456 // calculates a complex spectral difference detection function sample 436 double OnsetDetectionFunction :: complexSpectralDifference()
457 double OnsetDetectionFunction :: complex_spectral_difference()
458 { 437 {
459 double dev,pdev; 438 double dev,pdev;
460 double sum; 439 double sum;
461 double mag_diff,phase_diff; 440 double mag_diff,phase_diff;
462 double value; 441 double value;
463 442
464 // perform the FFT 443 // perform the FFT
465 perform_FFT(); 444 performFFT();
466 445
467 sum = 0; // initialise sum to zero 446 sum = 0; // initialise sum to zero
468 447
469 // compute phase values from fft output and sum deviations 448 // compute phase values from fft output and sum deviations
470 for (int i = 0;i < framesize;i++) 449 for (int i = 0;i < frameSize;i++)
471 { 450 {
472 // calculate phase value 451 // calculate phase value
473 phase[i] = atan2(out[i][1],out[i][0]); 452 phase[i] = atan2(complexOut[i][1],complexOut[i][0]);
474 453
475 // calculate magnitude value 454 // calculate magnitude value
476 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 455 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
477 456
478 457
479 // phase deviation 458 // phase deviation
480 dev = phase[i] - (2*phase_old[i]) + phase_old_2[i]; 459 dev = phase[i] - (2*prevPhase[i]) + prevPhase2[i];
481 460
482 // wrap into [-pi,pi] range 461 // wrap into [-pi,pi] range
483 pdev = princarg(dev); 462 pdev = princarg(dev);
484 463
485 464
486 // calculate magnitude difference (real part of Euclidean distance between complex frames) 465 // calculate magnitude difference (real part of Euclidean distance between complex frames)
487 mag_diff = mag[i] - mag_old[i]; 466 mag_diff = magSpec[i] - prevMagSpec[i];
488 467
489 // calculate phase difference (imaginary part of Euclidean distance between complex frames) 468 // calculate phase difference (imaginary part of Euclidean distance between complex frames)
490 phase_diff = -mag[i]*sin(pdev); 469 phase_diff = -magSpec[i]*sin(pdev);
491 470
492 471
493 472
494 // square real and imaginary parts, sum and take square root 473 // square real and imaginary parts, sum and take square root
495 value = sqrt(pow(mag_diff,2) + pow(phase_diff,2)); 474 value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
498 // add to sum 477 // add to sum
499 sum = sum + value; 478 sum = sum + value;
500 479
501 480
502 // store values for next calculation 481 // store values for next calculation
503 phase_old_2[i] = phase_old[i]; 482 prevPhase2[i] = prevPhase[i];
504 phase_old[i] = phase[i]; 483 prevPhase[i] = phase[i];
505 mag_old[i] = mag[i]; 484 prevMagSpec[i] = magSpec[i];
506 } 485 }
507 486
508 return sum; 487 return sum;
509 } 488 }
510 489
511 //-------------------------------------------------------------------------------------- 490 //=======================================================================
512 // calculates a complex spectral difference detection function sample (half-wave rectified) 491 double OnsetDetectionFunction :: complexSpectralDifferenceHWR()
513 double OnsetDetectionFunction :: complex_spectral_difference_hwr()
514 { 492 {
515 double dev,pdev; 493 double dev,pdev;
516 double sum; 494 double sum;
517 double mag_diff,phase_diff; 495 double mag_diff,phase_diff;
518 double value; 496 double value;
519 497
520 // perform the FFT 498 // perform the FFT
521 perform_FFT(); 499 performFFT();
522 500
523 sum = 0; // initialise sum to zero 501 sum = 0; // initialise sum to zero
524 502
525 // compute phase values from fft output and sum deviations 503 // compute phase values from fft output and sum deviations
526 for (int i = 0;i < framesize;i++) 504 for (int i = 0;i < frameSize;i++)
527 { 505 {
528 // calculate phase value 506 // calculate phase value
529 phase[i] = atan2(out[i][1],out[i][0]); 507 phase[i] = atan2(complexOut[i][1],complexOut[i][0]);
530 508
531 // calculate magnitude value 509 // calculate magnitude value
532 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 510 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
533 511
534 512
535 // phase deviation 513 // phase deviation
536 dev = phase[i] - (2*phase_old[i]) + phase_old_2[i]; 514 dev = phase[i] - (2*prevPhase[i]) + prevPhase2[i];
537 515
538 // wrap into [-pi,pi] range 516 // wrap into [-pi,pi] range
539 pdev = princarg(dev); 517 pdev = princarg(dev);
540 518
541 519
542 // calculate magnitude difference (real part of Euclidean distance between complex frames) 520 // calculate magnitude difference (real part of Euclidean distance between complex frames)
543 mag_diff = mag[i] - mag_old[i]; 521 mag_diff = magSpec[i] - prevMagSpec[i];
544 522
545 // if we have a positive change in magnitude, then include in sum, otherwise ignore (half-wave rectification) 523 // if we have a positive change in magnitude, then include in sum, otherwise ignore (half-wave rectification)
546 if (mag_diff > 0) 524 if (mag_diff > 0)
547 { 525 {
548 // calculate phase difference (imaginary part of Euclidean distance between complex frames) 526 // calculate phase difference (imaginary part of Euclidean distance between complex frames)
549 phase_diff = -mag[i]*sin(pdev); 527 phase_diff = -magSpec[i]*sin(pdev);
550 528
551 // square real and imaginary parts, sum and take square root 529 // square real and imaginary parts, sum and take square root
552 value = sqrt(pow(mag_diff,2) + pow(phase_diff,2)); 530 value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
553 531
554 // add to sum 532 // add to sum
555 sum = sum + value; 533 sum = sum + value;
556 } 534 }
557 535
558 // store values for next calculation 536 // store values for next calculation
559 phase_old_2[i] = phase_old[i]; 537 prevPhase2[i] = prevPhase[i];
560 phase_old[i] = phase[i]; 538 prevPhase[i] = phase[i];
561 mag_old[i] = mag[i]; 539 prevMagSpec[i] = magSpec[i];
562 } 540 }
563 541
564 return sum; 542 return sum;
565 } 543 }
566 544
567 545
568 //-------------------------------------------------------------------------------------- 546 //=======================================================================
569 // calculates a high frequency content detection function sample 547 double OnsetDetectionFunction :: highFrequencyContent()
570 double OnsetDetectionFunction :: high_frequency_content() 548 {
571 { 549 double sum;
572 double sum;
573 double mag_diff;
574 550
575 // perform the FFT 551 // perform the FFT
576 perform_FFT(); 552 performFFT();
577 553
578 sum = 0; // initialise sum to zero 554 sum = 0; // initialise sum to zero
579 555
580 // compute phase values from fft output and sum deviations 556 // compute phase values from fft output and sum deviations
581 for (int i = 0;i < framesize;i++) 557 for (int i = 0;i < frameSize;i++)
582 { 558 {
583 // calculate magnitude value 559 // calculate magnitude value
584 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 560 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
585 561
586 562
587 sum = sum + (mag[i]*((double) (i+1))); 563 sum = sum + (magSpec[i]*((double) (i+1)));
588 564
589 // store values for next calculation 565 // store values for next calculation
590 mag_old[i] = mag[i]; 566 prevMagSpec[i] = magSpec[i];
591 } 567 }
592 568
593 return sum; 569 return sum;
594 } 570 }
595 571
596 //-------------------------------------------------------------------------------------- 572 //=======================================================================
597 // calculates a high frequency spectral difference detection function sample 573 double OnsetDetectionFunction :: highFrequencySpectralDifference()
598 double OnsetDetectionFunction :: high_frequency_spectral_difference()
599 { 574 {
600 double sum; 575 double sum;
601 double mag_diff; 576 double mag_diff;
602 577
603 // perform the FFT 578 // perform the FFT
604 perform_FFT(); 579 performFFT();
605 580
606 sum = 0; // initialise sum to zero 581 sum = 0; // initialise sum to zero
607 582
608 // compute phase values from fft output and sum deviations 583 // compute phase values from fft output and sum deviations
609 for (int i = 0;i < framesize;i++) 584 for (int i = 0;i < frameSize;i++)
610 { 585 {
611 // calculate magnitude value 586 // calculate magnitude value
612 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 587 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
613 588
614 // calculate difference 589 // calculate difference
615 mag_diff = mag[i] - mag_old[i]; 590 mag_diff = magSpec[i] - prevMagSpec[i];
616 591
617 if (mag_diff < 0) 592 if (mag_diff < 0)
618 { 593 {
619 mag_diff = -mag_diff; 594 mag_diff = -mag_diff;
620 } 595 }
621 596
622 sum = sum + (mag_diff*((double) (i+1))); 597 sum = sum + (mag_diff*((double) (i+1)));
623 598
624 // store values for next calculation 599 // store values for next calculation
625 mag_old[i] = mag[i]; 600 prevMagSpec[i] = magSpec[i];
626 } 601 }
627 602
628 return sum; 603 return sum;
629 } 604 }
630 605
631 //-------------------------------------------------------------------------------------- 606 //=======================================================================
632 // calculates a high frequency spectral difference detection function sample (half-wave rectified) 607 double OnsetDetectionFunction :: highFrequencySpectralDifferenceHWR()
633 double OnsetDetectionFunction :: high_frequency_spectral_difference_hwr()
634 { 608 {
635 double sum; 609 double sum;
636 double mag_diff; 610 double mag_diff;
637 611
638 // perform the FFT 612 // perform the FFT
639 perform_FFT(); 613 performFFT();
640 614
641 sum = 0; // initialise sum to zero 615 sum = 0; // initialise sum to zero
642 616
643 // compute phase values from fft output and sum deviations 617 // compute phase values from fft output and sum deviations
644 for (int i = 0;i < framesize;i++) 618 for (int i = 0;i < frameSize;i++)
645 { 619 {
646 // calculate magnitude value 620 // calculate magnitude value
647 mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2)); 621 magSpec[i] = sqrt(pow(complexOut[i][0],2) + pow(complexOut[i][1],2));
648 622
649 // calculate difference 623 // calculate difference
650 mag_diff = mag[i] - mag_old[i]; 624 mag_diff = magSpec[i] - prevMagSpec[i];
651 625
652 if (mag_diff > 0) 626 if (mag_diff > 0)
653 { 627 {
654 sum = sum + (mag_diff*((double) (i+1))); 628 sum = sum + (mag_diff*((double) (i+1)));
655 } 629 }
656 630
657 // store values for next calculation 631 // store values for next calculation
658 mag_old[i] = mag[i]; 632 prevMagSpec[i] = magSpec[i];
659 } 633 }
660 634
661 return sum; 635 return sum;
662 } 636 }
663 637
664 638
665 //////////////////////////////////////////////////////////////////////////////////////////////// 639 ////////////////////////////////////////////////////////////////////////////////////////////////
666 //////////////////////////////////////////////////////////////////////////////////////////////// 640 ////////////////////////////////////////////////////////////////////////////////////////////////
667 ////////////////////////////// Methods to Calculate Windows //////////////////////////////////// 641 ////////////////////////////// Methods to Calculate Windows ////////////////////////////////////
668 642
669 //-------------------------------------------------------------------------------------- 643 //=======================================================================
670 // HANNING: set the window in the buffer 'window' to a Hanning window 644 void OnsetDetectionFunction :: calculateHanningWindow()
671 void OnsetDetectionFunction :: set_win_hanning()
672 { 645 {
673 double N; // variable to store framesize minus 1 646 double N; // variable to store framesize minus 1
674 647
675 N = (double) (framesize-1); // framesize minus 1 648 N = (double) (frameSize-1); // framesize minus 1
676 649
677 // Hanning window calculation 650 // Hanning window calculation
678 for (int n = 0;n < framesize;n++) 651 for (int n = 0;n < frameSize;n++)
679 { 652 {
680 window[n] = 0.5*(1-cos(2*pi*(n/N))); 653 window[n] = 0.5*(1-cos(2*pi*(n/N)));
681 } 654 }
682 } 655 }
683 656
684 //-------------------------------------------------------------------------------------- 657 //=======================================================================
685 // HAMMING: set the window in the buffer 'window' to a Hanning window 658 void OnsetDetectionFunction :: calclulateHammingWindow()
686 void OnsetDetectionFunction :: set_win_hamming()
687 { 659 {
688 double N; // variable to store framesize minus 1 660 double N; // variable to store framesize minus 1
689 double n_val; // double version of index 'n' 661 double n_val; // double version of index 'n'
690 662
691 N = (double) (framesize-1); // framesize minus 1 663 N = (double) (frameSize-1); // framesize minus 1
692 n_val = 0; 664 n_val = 0;
693 665
694 // Hamming window calculation 666 // Hamming window calculation
695 for (int n = 0;n < framesize;n++) 667 for (int n = 0;n < frameSize;n++)
696 { 668 {
697 window[n] = 0.54 - (0.46*cos(2*pi*(n_val/N))); 669 window[n] = 0.54 - (0.46*cos(2*pi*(n_val/N)));
698 n_val = n_val+1; 670 n_val = n_val+1;
699 } 671 }
700 } 672 }
701 673
702 //-------------------------------------------------------------------------------------- 674 //=======================================================================
703 // BLACKMAN: set the window in the buffer 'window' to a Blackman window 675 void OnsetDetectionFunction :: calculateBlackmanWindow()
704 void OnsetDetectionFunction :: set_win_blackman()
705 { 676 {
706 double N; // variable to store framesize minus 1 677 double N; // variable to store framesize minus 1
707 double n_val; // double version of index 'n' 678 double n_val; // double version of index 'n'
708 679
709 N = (double) (framesize-1); // framesize minus 1 680 N = (double) (frameSize-1); // framesize minus 1
710 n_val = 0; 681 n_val = 0;
711 682
712 // Blackman window calculation 683 // Blackman window calculation
713 for (int n = 0;n < framesize;n++) 684 for (int n = 0;n < frameSize;n++)
714 { 685 {
715 window[n] = 0.42 - (0.5*cos(2*pi*(n_val/N))) + (0.08*cos(4*pi*(n_val/N))); 686 window[n] = 0.42 - (0.5*cos(2*pi*(n_val/N))) + (0.08*cos(4*pi*(n_val/N)));
716 n_val = n_val+1; 687 n_val = n_val+1;
717 } 688 }
718 } 689 }
719 690
720 //-------------------------------------------------------------------------------------- 691 //=======================================================================
721 // TUKEY: set the window in the buffer 'window' to a Tukey window 692 void OnsetDetectionFunction :: calculateTukeyWindow()
722 void OnsetDetectionFunction :: set_win_tukey()
723 { 693 {
724 double N; // variable to store framesize minus 1 694 double N; // variable to store framesize minus 1
725 double n_val; // double version of index 'n' 695 double n_val; // double version of index 'n'
726 double alpha; // alpha [default value = 0.5]; 696 double alpha; // alpha [default value = 0.5];
727 697
728 int lim1,lim2;
729
730 alpha = 0.5; 698 alpha = 0.5;
731 699
732 N = (double) (framesize-1); // framesize minus 1 700 N = (double) (frameSize-1); // framesize minus 1
733 701
734 // Tukey window calculation 702 // Tukey window calculation
735 703
736 n_val = (double) (-1*((framesize/2)))+1; 704 n_val = (double) (-1*((frameSize/2)))+1;
737 705
738 for (int n = 0;n < framesize;n++) // left taper 706 for (int n = 0;n < frameSize;n++) // left taper
739 { 707 {
740 if ((n_val >= 0) && (n_val <= (alpha*(N/2)))) 708 if ((n_val >= 0) && (n_val <= (alpha*(N/2))))
741 { 709 {
742 window[n] = 1.0; 710 window[n] = 1.0;
743 } 711 }
753 n_val = n_val+1; 721 n_val = n_val+1;
754 } 722 }
755 723
756 } 724 }
757 725
758 //-------------------------------------------------------------------------------------- 726 //=======================================================================
759 // RECTANGULAR: set the window in the buffer 'window' to a Rectangular window 727 void OnsetDetectionFunction :: calculateRectangularWindow()
760 void OnsetDetectionFunction :: set_win_rectangular()
761 { 728 {
762 // Rectangular window calculation 729 // Rectangular window calculation
763 for (int n = 0;n < framesize;n++) 730 for (int n = 0;n < frameSize;n++)
764 { 731 {
765 window[n] = 1.0; 732 window[n] = 1.0;
766 } 733 }
767 } 734 }
768 735
770 737
771 //////////////////////////////////////////////////////////////////////////////////////////////// 738 ////////////////////////////////////////////////////////////////////////////////////////////////
772 //////////////////////////////////////////////////////////////////////////////////////////////// 739 ////////////////////////////////////////////////////////////////////////////////////////////////
773 ///////////////////////////////// Other Handy Methods ////////////////////////////////////////// 740 ///////////////////////////////// Other Handy Methods //////////////////////////////////////////
774 741
775 //-------------------------------------------------------------------------------------- 742 //=======================================================================
776 // set phase values to the range [-pi,pi] 743 double OnsetDetectionFunction :: princarg(double phaseVal)
777 double OnsetDetectionFunction :: princarg(double phaseval)
778 { 744 {
779 // if phase value is less than or equal to -pi then add 2*pi 745 // if phase value is less than or equal to -pi then add 2*pi
780 while (phaseval <= (-pi)) 746 while (phaseVal <= (-pi))
781 { 747 {
782 phaseval = phaseval + (2*pi); 748 phaseVal = phaseVal + (2*pi);
783 } 749 }
784 750
785 // if phase value is larger than pi, then subtract 2*pi 751 // if phase value is larger than pi, then subtract 2*pi
786 while (phaseval > pi) 752 while (phaseVal > pi)
787 { 753 {
788 phaseval = phaseval - (2*pi); 754 phaseVal = phaseVal - (2*pi);
789 } 755 }
790 756
791 return phaseval; 757 return phaseVal;
792 } 758 }
793 759
794 760
795 761
796 762
797 763
798 764
799 765
800 766
801 767
802 768
803 769
804 770
805 771