andrew@0: /*
andrew@0:  *  OnsetDetectionFunction.cpp
andrew@0:  *  
andrew@0:  *
andrew@0:  *  Created by Adam Stark on 22/03/2011.
andrew@0:  *  Copyright 2011 Queen Mary University of London. All rights reserved.
andrew@0:  *
andrew@0:  */
andrew@0: 
andrew@0: #include <math.h>
andrew@0: #include <iostream>
andrew@0: #include "OnsetDetectionFunction.h"
andrew@0: #include "accFFT.h"
andrew@0: using namespace std;
andrew@0: 
andrew@0: //-------------------------------------------------------------------------------
andrew@0: // Constructor
andrew@0: OnsetDetectionFunction :: OnsetDetectionFunction()
andrew@0: {	
andrew@0: 	// indicate that we have not initialised yet
andrew@0: 	initialised = 0;		
andrew@0: 	
andrew@0: 	// set pi
andrew@0: 	pi = 3.14159265358979;	
andrew@0: 	
andrew@0: 	// initialise with hopsize = 512, framesize = 1024, complex spectral difference DF and hanning window
andrew@0: 	initialise(512,1024,6,1);	
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: //-------------------------------------------------------------------------------
andrew@0: // Constructor (with arguments)
andrew@0: OnsetDetectionFunction :: OnsetDetectionFunction(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
andrew@0: {	
andrew@0: 	// indicate that we have not initialised yet
andrew@0: 	initialised = 0;		
andrew@0: 	
andrew@0: 	// set pi
andrew@0: 	pi = 3.14159265358979;	
andrew@0: 	
andrew@0: 	// initialise with arguments to constructor
andrew@0: 	initialise(arg_hsize,arg_fsize,arg_df_type,arg_win_type);
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // Destructor
andrew@0: OnsetDetectionFunction :: ~OnsetDetectionFunction()
andrew@0: {
andrew@0: 	// destroy fft plan
andrew@0:     //fftw_destroy_plan(p);
andrew@0: 	//fftw_free(in); 
andrew@0: 	//fftw_free(out);
andrew@0:         
Venetian@8: //    fft->~accFFT();
Venetian@8: 	
Venetian@8: 	delete fft;
Venetian@8: 	fft = NULL;
Venetian@8: 	
andrew@0:     delete [] in;
andrew@0:     in = NULL;
andrew@0:     delete [] out;
andrew@0:     out = NULL;
andrew@0:     
andrew@0: 
andrew@0: 	
andrew@0: 	// deallocate memory
andrew@0: 	delete [] frame;
andrew@0: 	frame = NULL;	
andrew@0: 	delete [] window;
andrew@0: 	window = NULL;									
andrew@0: 	delete [] wframe;
andrew@0: 	wframe = NULL;											
andrew@0: 	delete [] mag;
andrew@0: 	mag = NULL;
andrew@0: 	delete [] mag_old;
andrew@0: 	mag_old = NULL;
andrew@0: 	delete [] phase;
andrew@0: 	phase = NULL;
andrew@0: 	delete [] phase_old;
andrew@0: 	phase_old = NULL;	
andrew@0: 	delete [] phase_old_2;
andrew@0: 	phase_old_2 = NULL;
andrew@0: }
andrew@0: 
andrew@0: //-------------------------------------------------------------------------------
andrew@0: // Initialisation
andrew@0: void OnsetDetectionFunction :: initialise(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
andrew@0: {
andrew@0: 	if (initialised == 1) // if we have already initialised some buffers and an FFT plan
andrew@0: 	{
andrew@0: 		//////////////////////////////////
andrew@0: 		// TIDY UP FIRST - If initialise is called after the class has been initialised
andrew@0: 		// then we want to free up memory and cancel existing FFT plans
andrew@0: 	
andrew@0: 		// destroy fft plan
andrew@0: 		//fftw_destroy_plan(p);
andrew@0: 		//fftw_free(in); 
andrew@0: 		//fftw_free(out);
andrew@0:         
andrew@0:         fft->~accFFT();
andrew@0:         delete [] in;
andrew@0:         in = NULL;
andrew@0:         delete [] out;
andrew@0:         out = NULL;
andrew@0: 	
andrew@0: 	
andrew@0: 		// deallocate memory
andrew@0: 		delete [] frame;
andrew@0: 		frame = NULL;	
andrew@0: 		delete [] window;
andrew@0: 		window = NULL;									
andrew@0: 		delete [] wframe;
andrew@0: 		wframe = NULL;											
andrew@0: 		delete [] mag;
andrew@0: 		mag = NULL;
andrew@0: 		delete [] mag_old;
andrew@0: 		mag_old = NULL;
andrew@0: 		delete [] phase;
andrew@0: 		phase = NULL;
andrew@0: 		delete [] phase_old;
andrew@0: 		phase_old = NULL;	
andrew@0: 		delete [] phase_old_2;
andrew@0: 		phase_old_2 = NULL;
andrew@0: 	
andrew@0: 		////// END TIDY UP ///////////////
andrew@0: 		//////////////////////////////////
andrew@0: 	}
andrew@0: 	
andrew@0: 	hopsize = arg_hsize; // set hopsize
andrew@0: 	framesize = arg_fsize; // set framesize
andrew@0: 	
andrew@0: 	df_type = arg_df_type; // set detection function type
andrew@0: 		
andrew@0: 	// initialise buffers
andrew@0: 	frame = new double[framesize];											
andrew@0: 	window = new double[framesize];	
andrew@0: 	wframe = new double[framesize];		
andrew@0: 	
andrew@0: 	mag = new double[framesize];											
andrew@0: 	mag_old = new double[framesize];
andrew@0: 	
andrew@0: 	phase = new double[framesize];
andrew@0: 	phase_old = new double[framesize];
andrew@0: 	phase_old_2 = new double[framesize];
andrew@0: 	
andrew@0: 	
andrew@0: 	// set the window to the specified type
andrew@0: 	switch (arg_win_type){
andrew@0: 		case 0:
andrew@0: 			set_win_rectangular();		// Rectangular window
andrew@0: 			break;	
andrew@0: 		case 1:	
andrew@0: 			set_win_hanning();			// Hanning Window
andrew@0: 			break;
andrew@0: 		case 2:
andrew@0: 			set_win_hamming();			// Hamming Window
andrew@0: 			break;
andrew@0: 		case 3:
andrew@0: 			set_win_blackman();			// Blackman Window
andrew@0: 			break;
andrew@0: 		case 4:
andrew@0: 			set_win_tukey();			// Tukey Window
andrew@0: 			break;
andrew@0: 		default:
andrew@0: 			set_win_hanning();			// DEFAULT: Hanning Window
andrew@0: 	}
andrew@0: 	
andrew@0: 	
andrew@0: 	
andrew@0: 	
andrew@0: 	// initialise previous magnitude spectrum to zero
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		mag_old[i] = 0.0;
andrew@0: 		phase_old[i] = 0.0;
andrew@0: 		phase_old_2[i] = 0.0;
andrew@0: 		frame[i] = 0.0;
andrew@0: 	}
andrew@0: 	
andrew@0: 	energy_sum_old = 0.0;	// initialise previous energy sum value to zero
andrew@0: 	
andrew@0: 	/*  Init fft */
andrew@0: 	//in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize);		// complex array to hold fft data
andrew@0: 	//out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize);	// complex array to hold fft data
andrew@0: 	//p = fftw_plan_dft_1d(framesize, in, out, FFTW_FORWARD, FFTW_ESTIMATE);	// FFT plan initialisation
andrew@0: 
andrew@0:     in = new double[framesize];
andrew@0:     out = new fft_complex[framesize];
andrew@0: 
andrew@0:     //out = new double[2][framesize];
andrew@0:     fft = new accFFT(framesize,1);
andrew@0: 	
andrew@0: 	initialised = 1;
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // set the detection function type to that specified by the argument
andrew@0: void OnsetDetectionFunction :: set_df_type(int arg_df_type)
andrew@0: {
andrew@0: 	df_type = arg_df_type; // set detection function type
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a single detection function sample from a single audio frame.
Venetian@8: double OnsetDetectionFunction :: getDFsample(float inputbuffer[])
Venetian@8: {	
Venetian@8: 	double df_sample;
Venetian@8: 	
Venetian@8: 	// shift audio samples back in frame by hop size
Venetian@8: 	for (int i = 0; i < (framesize-hopsize);i++)
Venetian@8: 	{
Venetian@8: 		frame[i] = frame[i+hopsize];
Venetian@8: 	}
Venetian@8: 	
Venetian@8: 	// add new samples to frame from input buffer
Venetian@8: 	int j = 0;
Venetian@8: 	for (int i = (framesize-hopsize);i < framesize;i++)
Venetian@8: 	{
Venetian@8: 		frame[i] = inputbuffer[j];
Venetian@8: 		j++;
Venetian@8: 	}
Venetian@8: 	
Venetian@8: 	switch (df_type){
Venetian@8: 		case 0:
Venetian@8: 			df_sample = energy_envelope();	// calculate energy envelope detection function sample
Venetian@8: 			break;	
Venetian@8: 		case 1:
Venetian@8: 			df_sample = energy_difference();	// calculate half-wave rectified energy difference detection function sample
Venetian@8: 			break;
Venetian@8: 		case 2:
Venetian@8: 			df_sample = spectral_difference();	// calculate spectral difference detection function sample
Venetian@8: 			break;
Venetian@8: 		case 3:
Venetian@8: 			df_sample = spectral_difference_hwr(); // calculate spectral difference detection function sample (half wave rectified)
Venetian@8: 			break;
Venetian@8: 		case 4:
Venetian@8: 			df_sample = phase_deviation();		// calculate phase deviation detection function sample (half wave rectified)
Venetian@8: 			break;
Venetian@8: 		case 5:
Venetian@8: 			df_sample = complex_spectral_difference(); // calcualte complex spectral difference detection function sample
Venetian@8: 			break;
Venetian@8: 		case 6:
Venetian@8: 			df_sample = complex_spectral_difference_hwr();  // calcualte complex spectral difference detection function sample (half-wave rectified)
Venetian@8: 			break;
Venetian@8: 		case 7:
Venetian@8: 			df_sample = high_frequency_content(); // calculate high frequency content detection function sample
Venetian@8: 			break;
Venetian@8: 		case 8:
Venetian@8: 			df_sample = high_frequency_spectral_difference(); // calculate high frequency spectral difference detection function sample
Venetian@8: 			break;
Venetian@8: 		case 9:
Venetian@8: 			df_sample = high_frequency_spectral_difference_hwr(); // calculate high frequency spectral difference detection function (half-wave rectified)
Venetian@8: 			break;
Venetian@8: 		default:
Venetian@8: 			df_sample = 1.0;
Venetian@8: 	}
Venetian@8: 	
Venetian@8: 	return df_sample;
Venetian@8: }
Venetian@8: 
Venetian@8: 	
Venetian@8: 	
Venetian@8: 	
Venetian@8: //--------------------------------------------------------------------------------------
Venetian@8: // calculates a single detection function sample from a single audio frame.
andrew@0: double OnsetDetectionFunction :: getDFsample(double inputbuffer[])
andrew@0: {	
andrew@0: 	double df_sample;
andrew@0: 		
andrew@0: 	// shift audio samples back in frame by hop size
andrew@0: 	for (int i = 0; i < (framesize-hopsize);i++)
andrew@0: 	{
andrew@0: 		frame[i] = frame[i+hopsize];
andrew@0: 	}
andrew@0: 	
andrew@0: 	// add new samples to frame from input buffer
andrew@0: 	int j = 0;
andrew@0: 	for (int i = (framesize-hopsize);i < framesize;i++)
andrew@0: 	{
andrew@0: 		frame[i] = inputbuffer[j];
andrew@0: 		j++;
andrew@0: 	}
andrew@0: 		
andrew@0: 	switch (df_type){
andrew@0: 		case 0:
andrew@0: 			df_sample = energy_envelope();	// calculate energy envelope detection function sample
andrew@0: 			break;	
andrew@0: 		case 1:
andrew@0: 			df_sample = energy_difference();	// calculate half-wave rectified energy difference detection function sample
andrew@0: 			break;
andrew@0: 		case 2:
andrew@0: 			df_sample = spectral_difference();	// calculate spectral difference detection function sample
andrew@0: 			break;
andrew@0: 		case 3:
andrew@0: 			df_sample = spectral_difference_hwr(); // calculate spectral difference detection function sample (half wave rectified)
andrew@0: 			break;
andrew@0: 		case 4:
andrew@0: 			df_sample = phase_deviation();		// calculate phase deviation detection function sample (half wave rectified)
andrew@0: 			break;
andrew@0: 		case 5:
andrew@0: 			df_sample = complex_spectral_difference(); // calcualte complex spectral difference detection function sample
andrew@0: 			break;
andrew@0: 		case 6:
andrew@0: 			df_sample = complex_spectral_difference_hwr();  // calcualte complex spectral difference detection function sample (half-wave rectified)
andrew@0: 			break;
andrew@0: 		case 7:
andrew@0: 			df_sample = high_frequency_content(); // calculate high frequency content detection function sample
andrew@0: 			break;
andrew@0: 		case 8:
andrew@0: 			df_sample = high_frequency_spectral_difference(); // calculate high frequency spectral difference detection function sample
andrew@0: 			break;
andrew@0: 		case 9:
andrew@0: 			df_sample = high_frequency_spectral_difference_hwr(); // calculate high frequency spectral difference detection function (half-wave rectified)
andrew@0: 			break;
andrew@0: 		default:
andrew@0: 			df_sample = 1.0;
andrew@0: 	}
andrew@0: 		
andrew@0: 	return df_sample;
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // performs the fft, storing the complex result in 'out'
andrew@0: void OnsetDetectionFunction :: perform_FFT()
andrew@0: {
andrew@0: 	int fsize2 = (framesize/2);
andrew@0: 	
andrew@0:     
andrew@0: 	// window frame
andrew@0: 	//for (int i = 0;i < fsize2;i++)
andrew@0:     for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0:         in[i] = frame[i]*window[i];
andrew@0:         
andrew@0:         //in[i] = frame[i+fsize2]*window[i+fsize2];
andrew@0:         //in[i+fsize2] = frame[i] * window[i];
andrew@0:         
andrew@0:         /*
andrew@0: 		in[i][0] = frame[i+fsize2] * window[i+fsize2];
andrew@0: 		in[i][1] = 0.0;
andrew@0: 		in[i+fsize2][0] = frame[i] * window[i];
andrew@0: 		in[i+fsize2][1] = 0.0;
andrew@0:          */
andrew@0: 	}
andrew@0: 	
andrew@0: 	// perform the fft
andrew@0: 	//fftw_execute(p);
andrew@0:     
andrew@0:     // FIX  
andrew@0:     fft->forward_FFT_d(in,out);
andrew@0: }
andrew@0: 
andrew@0: ////////////////////////////////////////////////////////////////////////////////////////////////
andrew@0: ////////////////////////////////////////////////////////////////////////////////////////////////
andrew@0: ////////////////////////////// Methods for Detection Functions /////////////////////////////////
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates an energy envelope detection function sample
andrew@0: double OnsetDetectionFunction :: energy_envelope()
andrew@0: {
andrew@0: 	double sum;
andrew@0: 	
andrew@0: 	sum = 0;	// initialise sum
andrew@0: 	
andrew@0: 	// sum the squares of the samples
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		sum = sum + (frame[i]*frame[i]);
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		// return sum
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a half-wave rectified energy difference detection function sample
andrew@0: double OnsetDetectionFunction :: energy_difference()
andrew@0: {
andrew@0: 	double sum;
andrew@0: 	double sample;
andrew@0: 	
andrew@0: 	sum = 0;	// initialise sum
andrew@0: 	
andrew@0: 	// sum the squares of the samples
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		sum = sum + (frame[i]*frame[i]);
andrew@0: 	}
andrew@0: 	
andrew@0: 	sample = sum - energy_sum_old;	// sample is first order difference in energy
andrew@0: 	
andrew@0: 	energy_sum_old = sum;	// store energy value for next calculation
andrew@0: 	
andrew@0: 	if (sample > 0)
andrew@0: 	{
andrew@0: 		return sample;		// return difference
andrew@0: 	}
andrew@0: 	else
andrew@0: 	{
andrew@0: 		return 0;
andrew@0: 	}
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a spectral difference detection function sample
andrew@0: double OnsetDetectionFunction :: spectral_difference()
andrew@0: {
andrew@0: 	double diff;
andrew@0: 	double sum;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	// compute first (N/2)+1 mag values
andrew@0: 	for (int i = 0;i < (framesize/2)+1;i++)
andrew@0: 	{
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 	}
andrew@0: 	// mag spec symmetric above (N/2)+1 so copy previous values
andrew@0: 	for (int i = (framesize/2)+1;i < framesize;i++)
andrew@0: 	{
andrew@0: 		mag[i] = mag[framesize-i];		
andrew@0: 	}
andrew@0: 	
andrew@0: 	sum = 0;	// initialise sum to zero
andrew@0: 
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		// calculate difference
andrew@0: 		diff = mag[i] - mag_old[i];
andrew@0: 		
andrew@0: 		// ensure all difference values are positive
andrew@0: 		if (diff < 0)
andrew@0: 		{
andrew@0: 			diff = diff*-1;
andrew@0: 		}
andrew@0: 		
andrew@0: 		// add difference to sum
andrew@0: 		sum = sum+diff;
andrew@0: 		
andrew@0: 		// store magnitude spectrum bin for next detection function sample calculation
andrew@0: 		mag_old[i] = mag[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a spectral difference detection function sample
andrew@0: double OnsetDetectionFunction :: spectral_difference_hwr()
andrew@0: {
andrew@0: 	double diff;
andrew@0: 	double sum;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	// compute first (N/2)+1 mag values
andrew@0: 	for (int i = 0;i < (framesize/2)+1;i++)
andrew@0: 	{
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 	}
andrew@0: 	// mag spec symmetric above (N/2)+1 so copy previous values
andrew@0: 	for (int i = (framesize/2)+1;i < framesize;i++)
andrew@0: 	{
andrew@0: 		mag[i] = mag[framesize-i];		
andrew@0: 	}
andrew@0: 	
andrew@0: 	sum = 0;	// initialise sum to zero
andrew@0: 	
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		// calculate difference
andrew@0: 		diff = mag[i] - mag_old[i];
andrew@0: 		
andrew@0: 		// only add up positive differences
andrew@0: 		if (diff > 0)
andrew@0: 		{
andrew@0: 			// add difference to sum
andrew@0: 			sum = sum+diff;
andrew@0: 		}
andrew@0: 		
andrew@0: 		
andrew@0: 		
andrew@0: 		// store magnitude spectrum bin for next detection function sample calculation
andrew@0: 		mag_old[i] = mag[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a phase deviation detection function sample
andrew@0: double OnsetDetectionFunction :: phase_deviation()
andrew@0: {
andrew@0: 	double dev,pdev;
andrew@0: 	double sum;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	sum = 0; // initialise sum to zero
andrew@0: 	
andrew@0: 	// compute phase values from fft output and sum deviations
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		// calculate phase value
andrew@0: 		phase[i] = atan2(out[i][1],out[i][0]);
andrew@0: 		
andrew@0: 		// calculate magnitude value
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 		
andrew@0: 		
andrew@0: 		// if bin is not just a low energy bin then examine phase deviation
andrew@0: 		if (mag[i] > 0.1)
andrew@0: 		{
andrew@0: 			dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];	// phase deviation
andrew@0: 			pdev = princarg(dev);	// wrap into [-pi,pi] range
andrew@0: 		
andrew@0: 			// make all values positive
andrew@0: 			if (pdev < 0)	
andrew@0: 			{
andrew@0: 				pdev = pdev*-1;
andrew@0: 			}
andrew@0: 						
andrew@0: 			// add to sum
andrew@0: 			sum = sum + pdev;
andrew@0: 		}
andrew@0: 				
andrew@0: 		// store values for next calculation
andrew@0: 		phase_old_2[i] = phase_old[i];
andrew@0: 		phase_old[i] = phase[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a complex spectral difference detection function sample
andrew@0: double OnsetDetectionFunction :: complex_spectral_difference()
andrew@0: {
andrew@0: 	double dev,pdev;
andrew@0: 	double sum;
andrew@0: 	double mag_diff,phase_diff;
andrew@0: 	double value;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	sum = 0; // initialise sum to zero
andrew@0: 	
andrew@0: 	// compute phase values from fft output and sum deviations
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		// calculate phase value
andrew@0: 		phase[i] = atan2(out[i][1],out[i][0]);
andrew@0: 		
andrew@0: 		// calculate magnitude value
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 		
andrew@0: 		
andrew@0: 		// phase deviation
andrew@0: 		dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];	
andrew@0: 		
andrew@0: 		// wrap into [-pi,pi] range
andrew@0: 		pdev = princarg(dev);	
andrew@0: 		
andrew@0: 		
andrew@0: 		// calculate magnitude difference (real part of Euclidean distance between complex frames)
andrew@0: 		mag_diff = mag[i] - mag_old[i];
andrew@0: 		
andrew@0: 		// calculate phase difference (imaginary part of Euclidean distance between complex frames)
andrew@0: 		phase_diff = -mag[i]*sin(pdev);
andrew@0: 		
andrew@0: 
andrew@0: 		
andrew@0: 		// square real and imaginary parts, sum and take square root
andrew@0: 		value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
andrew@0: 	
andrew@0: 			
andrew@0: 		// add to sum
andrew@0: 		sum = sum + value;
andrew@0: 		
andrew@0: 		
andrew@0: 		// store values for next calculation
andrew@0: 		phase_old_2[i] = phase_old[i];
andrew@0: 		phase_old[i] = phase[i];
andrew@0: 		mag_old[i] = mag[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a complex spectral difference detection function sample (half-wave rectified)
andrew@0: double OnsetDetectionFunction :: complex_spectral_difference_hwr()
andrew@0: {
andrew@0: 	double dev,pdev;
andrew@0: 	double sum;
andrew@0: 	double mag_diff,phase_diff;
andrew@0: 	double value;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	sum = 0; // initialise sum to zero
andrew@0: 	
andrew@0: 	// compute phase values from fft output and sum deviations
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{
andrew@0: 		// calculate phase value
andrew@0: 		phase[i] = atan2(out[i][1],out[i][0]);
andrew@0: 		
andrew@0: 		// calculate magnitude value
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 		
andrew@0: 		
andrew@0: 		// phase deviation
andrew@0: 		dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];	
andrew@0: 		
andrew@0: 		// wrap into [-pi,pi] range
andrew@0: 		pdev = princarg(dev);	
andrew@0: 		
andrew@0: 		
andrew@0: 		// calculate magnitude difference (real part of Euclidean distance between complex frames)
andrew@0: 		mag_diff = mag[i] - mag_old[i];
andrew@0: 		
andrew@0: 		// if we have a positive change in magnitude, then include in sum, otherwise ignore (half-wave rectification)
andrew@0: 		if (mag_diff > 0)
andrew@0: 		{
andrew@0: 			// calculate phase difference (imaginary part of Euclidean distance between complex frames)
andrew@0: 			phase_diff = -mag[i]*sin(pdev);
andrew@0: 
andrew@0: 			// square real and imaginary parts, sum and take square root
andrew@0: 			value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
andrew@0: 		
andrew@0: 			// add to sum
andrew@0: 			sum = sum + value;
andrew@0: 		}
andrew@0: 		
andrew@0: 		// store values for next calculation
andrew@0: 		phase_old_2[i] = phase_old[i];
andrew@0: 		phase_old[i] = phase[i];
andrew@0: 		mag_old[i] = mag[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a high frequency content detection function sample
andrew@0: double OnsetDetectionFunction :: high_frequency_content()
andrew@0: {
andrew@0: 	double sum;
andrew@0: 	double mag_diff;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	sum = 0; // initialise sum to zero
andrew@0: 	
andrew@0: 	// compute phase values from fft output and sum deviations
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{		
andrew@0: 		// calculate magnitude value
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 		
andrew@0: 		
andrew@0: 		sum = sum + (mag[i]*((double) (i+1)));
andrew@0: 		
andrew@0: 		// store values for next calculation
andrew@0: 		mag_old[i] = mag[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a high frequency spectral difference detection function sample
andrew@0: double OnsetDetectionFunction :: high_frequency_spectral_difference()
andrew@0: {
andrew@0: 	double sum;
andrew@0: 	double mag_diff;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	sum = 0; // initialise sum to zero
andrew@0: 	
andrew@0: 	// compute phase values from fft output and sum deviations
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{		
andrew@0: 		// calculate magnitude value
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 		
andrew@0: 		// calculate difference
andrew@0: 		mag_diff = mag[i] - mag_old[i];
andrew@0: 		
andrew@0: 		if (mag_diff < 0)
andrew@0: 		{
andrew@0: 			mag_diff = -mag_diff;
andrew@0: 		}
andrew@0: 		
andrew@0: 		sum = sum + (mag_diff*((double) (i+1)));
andrew@0: 		
andrew@0: 		// store values for next calculation
andrew@0: 		mag_old[i] = mag[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // calculates a high frequency spectral difference detection function sample (half-wave rectified)
andrew@0: double OnsetDetectionFunction :: high_frequency_spectral_difference_hwr()
andrew@0: {
andrew@0: 	double sum;
andrew@0: 	double mag_diff;
andrew@0: 	
andrew@0: 	// perform the FFT
andrew@0: 	perform_FFT();
andrew@0: 	
andrew@0: 	sum = 0; // initialise sum to zero
andrew@0: 	
andrew@0: 	// compute phase values from fft output and sum deviations
andrew@0: 	for (int i = 0;i < framesize;i++)
andrew@0: 	{		
andrew@0: 		// calculate magnitude value
andrew@0: 		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
andrew@0: 		
andrew@0: 		// calculate difference
andrew@0: 		mag_diff = mag[i] - mag_old[i];
andrew@0: 		
andrew@0: 		if (mag_diff > 0)
andrew@0: 		{
andrew@0: 			sum = sum + (mag_diff*((double) (i+1)));
andrew@0: 		}
andrew@0: 
andrew@0: 		// store values for next calculation
andrew@0: 		mag_old[i] = mag[i];
andrew@0: 	}
andrew@0: 	
andrew@0: 	return sum;		
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: ////////////////////////////////////////////////////////////////////////////////////////////////
andrew@0: ////////////////////////////////////////////////////////////////////////////////////////////////
andrew@0: ////////////////////////////// Methods to Calculate Windows ////////////////////////////////////
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // HANNING: set the window in the buffer 'window' to a Hanning window
andrew@0: void OnsetDetectionFunction :: set_win_hanning()
andrew@0: {
andrew@0: 	double N;		// variable to store framesize minus 1
andrew@0: 	
andrew@0: 	N = (double) (framesize-1);	// framesize minus 1
andrew@0: 	
andrew@0: 	// Hanning window calculation
andrew@0: 	for (int n = 0;n < framesize;n++)
andrew@0: 	{
andrew@0: 		window[n] = 0.5*(1-cos(2*pi*(n/N)));
andrew@0: 	}
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // HAMMING: set the window in the buffer 'window' to a Hanning window
andrew@0: void OnsetDetectionFunction :: set_win_hamming()
andrew@0: {
andrew@0: 	double N;		// variable to store framesize minus 1
andrew@0: 	double n_val;	// double version of index 'n'
andrew@0: 	
andrew@0: 	N = (double) (framesize-1);	// framesize minus 1
andrew@0: 	n_val = 0;
andrew@0: 	
andrew@0: 	// Hamming window calculation
andrew@0: 	for (int n = 0;n < framesize;n++)
andrew@0: 	{
andrew@0: 		window[n] = 0.54 - (0.46*cos(2*pi*(n_val/N)));
andrew@0: 		n_val = n_val+1;
andrew@0: 	}
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // BLACKMAN: set the window in the buffer 'window' to a Blackman window
andrew@0: void OnsetDetectionFunction :: set_win_blackman()
andrew@0: {
andrew@0: 	double N;		// variable to store framesize minus 1
andrew@0: 	double n_val;	// double version of index 'n'
andrew@0: 	
andrew@0: 	N = (double) (framesize-1);	// framesize minus 1
andrew@0: 	n_val = 0;
andrew@0: 	
andrew@0: 	// Blackman window calculation
andrew@0: 	for (int n = 0;n < framesize;n++)
andrew@0: 	{
andrew@0: 		window[n] = 0.42 - (0.5*cos(2*pi*(n_val/N))) + (0.08*cos(4*pi*(n_val/N)));
andrew@0: 		n_val = n_val+1;
andrew@0: 	}
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // TUKEY: set the window in the buffer 'window' to a Tukey window
andrew@0: void OnsetDetectionFunction :: set_win_tukey()
andrew@0: {
andrew@0: 	double N;		// variable to store framesize minus 1
andrew@0: 	double n_val;	// double version of index 'n'
andrew@0: 	double alpha;	// alpha [default value = 0.5];
andrew@0: 	
andrew@0: 	int lim1,lim2;
andrew@0: 	
andrew@0: 	alpha = 0.5;
andrew@0: 	
andrew@0: 	N = (double) (framesize-1);	// framesize minus 1
andrew@0: 		
andrew@0: 	// Tukey window calculation
andrew@0: 	
andrew@0: 	n_val = (double) (-1*((framesize/2)))+1;
andrew@0: 
andrew@0: 	for (int n = 0;n < framesize;n++)	// left taper
andrew@0: 	{
andrew@0: 		if ((n_val >= 0) && (n_val <= (alpha*(N/2))))
andrew@0: 		{
andrew@0: 			window[n] = 1.0;
andrew@0: 		}
andrew@0: 		else if ((n_val <= 0) && (n_val >= (-1*alpha*(N/2))))
andrew@0: 		{
andrew@0: 			window[n] = 1.0;
andrew@0: 		}
andrew@0: 		else
andrew@0: 		{
andrew@0: 			window[n] = 0.5*(1+cos(pi*(((2*n_val)/(alpha*N))-1)));
andrew@0: 		}
andrew@0: 
andrew@0: 		n_val = n_val+1;			 
andrew@0: 	}
andrew@0: 
andrew@0: }
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // RECTANGULAR: set the window in the buffer 'window' to a Rectangular window
andrew@0: void OnsetDetectionFunction :: set_win_rectangular()
andrew@0: {
andrew@0: 	// Rectangular window calculation
andrew@0: 	for (int n = 0;n < framesize;n++)
andrew@0: 	{
andrew@0: 		window[n] = 1.0;
andrew@0: 	}
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: ////////////////////////////////////////////////////////////////////////////////////////////////
andrew@0: ////////////////////////////////////////////////////////////////////////////////////////////////
andrew@0: ///////////////////////////////// Other Handy Methods //////////////////////////////////////////
andrew@0: 
andrew@0: //--------------------------------------------------------------------------------------
andrew@0: // set phase values to the range [-pi,pi]
andrew@0: double OnsetDetectionFunction :: princarg(double phaseval)
andrew@0: {	
andrew@0: 	// if phase value is less than or equal to -pi then add 2*pi
andrew@0: 	while (phaseval <= (-pi)) 
andrew@0: 	{
andrew@0: 		phaseval = phaseval + (2*pi);
andrew@0: 	}
andrew@0: 	
andrew@0: 	// if phase value is larger than pi, then subtract 2*pi
andrew@0: 	while (phaseval > pi)
andrew@0: 	{
andrew@0: 		phaseval = phaseval - (2*pi);
andrew@0: 	}
andrew@0: 			
andrew@0: 	return phaseval;
andrew@0: }
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: 
andrew@0: