Mercurial > hg > btrack
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 |