view src/OnsetDetectionFunction.cpp @ 52:45231107c9d6

Reformatted comments, removed the OnsetDetectionFunction constructor with no arguments, removed a number of unused variables and made changes to the python module to fix some casting problems and removed some unused variables there also. Still getting the same results, so no overall changes to the algorithm.

//=======================================================================
/** @file OnsetDetectionFunction.cpp
 *  @brief A class for calculating onset detection functions
 *  @author Adam Stark
 *  @copyright Copyright (C) 2008-2014  Queen Mary University of London
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
//=======================================================================

#include <math.h>
#include "OnsetDetectionFunction.h"

//=======================================================================
OnsetDetectionFunction :: OnsetDetectionFunction(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
{	
	// indicate that we have not initialised yet
	initialised = 0;		
	
	// set pi
	pi = 3.14159265358979;	
	
	// initialise with arguments to constructor
	initialise(arg_hsize,arg_fsize,arg_df_type,arg_win_type);
}


//=======================================================================
OnsetDetectionFunction :: ~OnsetDetectionFunction()
{
	// destroy fft plan
    fftw_destroy_plan(p);
	fftw_free(in); 
	fftw_free(out);
	
	// deallocate memory
	delete [] frame;
	frame = NULL;	
	delete [] window;
	window = NULL;									
	delete [] wframe;
	wframe = NULL;											
	delete [] mag;
	mag = NULL;
	delete [] mag_old;
	mag_old = NULL;
	delete [] phase;
	phase = NULL;
	delete [] phase_old;
	phase_old = NULL;	
	delete [] phase_old_2;
	phase_old_2 = NULL;
}

//=======================================================================
void OnsetDetectionFunction :: initialise(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
{
	if (initialised == 1) // if we have already initialised some buffers and an FFT plan
	{
		//////////////////////////////////
		// TIDY UP FIRST - If initialise is called after the class has been initialised
		// then we want to free up memory and cancel existing FFT plans
	
		// destroy fft plan
		fftw_destroy_plan(p);
		fftw_free(in); 
		fftw_free(out);
	
	
		// deallocate memory
		delete [] frame;
		frame = NULL;	
		delete [] window;
		window = NULL;									
		delete [] wframe;
		wframe = NULL;											
		delete [] mag;
		mag = NULL;
		delete [] mag_old;
		mag_old = NULL;
		delete [] phase;
		phase = NULL;
		delete [] phase_old;
		phase_old = NULL;	
		delete [] phase_old_2;
		phase_old_2 = NULL;
	
		////// END TIDY UP ///////////////
		//////////////////////////////////
	}
	
	hopsize = arg_hsize; // set hopsize
	framesize = arg_fsize; // set framesize
	
	df_type = arg_df_type; // set detection function type
		
	// initialise buffers
	frame = new double[framesize];											
	window = new double[framesize];	
	wframe = new double[framesize];		
	
	mag = new double[framesize];											
	mag_old = new double[framesize];
	
	phase = new double[framesize];
	phase_old = new double[framesize];
	phase_old_2 = new double[framesize];
	
	
	// set the window to the specified type
	switch (arg_win_type){
		case 0:
			set_win_rectangular();		// Rectangular window
			break;	
		case 1:	
			set_win_hanning();			// Hanning Window
			break;
		case 2:
			set_win_hamming();			// Hamming Window
			break;
		case 3:
			set_win_blackman();			// Blackman Window
			break;
		case 4:
			set_win_tukey();			// Tukey Window
			break;
		default:
			set_win_hanning();			// DEFAULT: Hanning Window
	}
	
	
	
	
	// initialise previous magnitude spectrum to zero
	for (int i = 0;i < framesize;i++)
	{
		mag_old[i] = 0.0;
		phase_old[i] = 0.0;
		phase_old_2[i] = 0.0;
		frame[i] = 0.0;
	}
	
	energy_sum_old = 0.0;	// initialise previous energy sum value to zero
	
	/*  Init fft */
	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize);		// complex array to hold fft data
	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize);	// complex array to hold fft data
	p = fftw_plan_dft_1d(framesize, in, out, FFTW_FORWARD, FFTW_ESTIMATE);	// FFT plan initialisation
	
	initialised = 1;
}

//=======================================================================
void OnsetDetectionFunction :: set_df_type(int arg_df_type)
{
	df_type = arg_df_type; // set detection function type
}

//=======================================================================
double OnsetDetectionFunction :: getDFsample(double inputbuffer[])
{	
	double df_sample;
		
	// shift audio samples back in frame by hop size
	for (int i = 0; i < (framesize-hopsize);i++)
	{
		frame[i] = frame[i+hopsize];
	}
	
	// add new samples to frame from input buffer
	int j = 0;
	for (int i = (framesize-hopsize);i < framesize;i++)
	{
		frame[i] = inputbuffer[j];
		j++;
	}
		
	switch (df_type){
		case 0:
			df_sample = energy_envelope();	// calculate energy envelope detection function sample
			break;	
		case 1:
			df_sample = energy_difference();	// calculate half-wave rectified energy difference detection function sample
			break;
		case 2:
			df_sample = spectral_difference();	// calculate spectral difference detection function sample
			break;
		case 3:
			df_sample = spectral_difference_hwr(); // calculate spectral difference detection function sample (half wave rectified)
			break;
		case 4:
			df_sample = phase_deviation();		// calculate phase deviation detection function sample (half wave rectified)
			break;
		case 5:
			df_sample = complex_spectral_difference(); // calcualte complex spectral difference detection function sample
			break;
		case 6:
			df_sample = complex_spectral_difference_hwr();  // calcualte complex spectral difference detection function sample (half-wave rectified)
			break;
		case 7:
			df_sample = high_frequency_content(); // calculate high frequency content detection function sample
			break;
		case 8:
			df_sample = high_frequency_spectral_difference(); // calculate high frequency spectral difference detection function sample
			break;
		case 9:
			df_sample = high_frequency_spectral_difference_hwr(); // calculate high frequency spectral difference detection function (half-wave rectified)
			break;
		default:
			df_sample = 1.0;
	}
		
	return df_sample;
}


//=======================================================================
void OnsetDetectionFunction :: perform_FFT()
{
	int fsize2 = (framesize/2);
	
	// window frame and copy to complex array, swapping the first and second half of the signal
	for (int i = 0;i < fsize2;i++)
	{
		in[i][0] = frame[i+fsize2] * window[i+fsize2];
		in[i][1] = 0.0;
		in[i+fsize2][0] = frame[i] * window[i];
		in[i+fsize2][1] = 0.0;
	}
	
	// perform the fft
	fftw_execute(p);
}

////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////// Methods for Detection Functions /////////////////////////////////

//=======================================================================
double OnsetDetectionFunction :: energy_envelope()
{
	double sum;
	
	sum = 0;	// initialise sum
	
	// sum the squares of the samples
	for (int i = 0;i < framesize;i++)
	{
		sum = sum + (frame[i]*frame[i]);
	}
	
	return sum;		// return sum
}

//=======================================================================
double OnsetDetectionFunction :: energy_difference()
{
	double sum;
	double sample;
	
	sum = 0;	// initialise sum
	
	// sum the squares of the samples
	for (int i = 0;i < framesize;i++)
	{
		sum = sum + (frame[i]*frame[i]);
	}
	
	sample = sum - energy_sum_old;	// sample is first order difference in energy
	
	energy_sum_old = sum;	// store energy value for next calculation
	
	if (sample > 0)
	{
		return sample;		// return difference
	}
	else
	{
		return 0;
	}
}

//=======================================================================
double OnsetDetectionFunction :: spectral_difference()
{
	double diff;
	double sum;
	
	// perform the FFT
	perform_FFT();
	
	// compute first (N/2)+1 mag values
	for (int i = 0;i < (framesize/2)+1;i++)
	{
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
	}
	// mag spec symmetric above (N/2)+1 so copy previous values
	for (int i = (framesize/2)+1;i < framesize;i++)
	{
		mag[i] = mag[framesize-i];		
	}
	
	sum = 0;	// initialise sum to zero

	for (int i = 0;i < framesize;i++)
	{
		// calculate difference
		diff = mag[i] - mag_old[i];
		
		// ensure all difference values are positive
		if (diff < 0)
		{
			diff = diff*-1;
		}
		
		// add difference to sum
		sum = sum+diff;
		
		// store magnitude spectrum bin for next detection function sample calculation
		mag_old[i] = mag[i];
	}
	
	return sum;		
}

//=======================================================================
double OnsetDetectionFunction :: spectral_difference_hwr()
{
	double diff;
	double sum;
	
	// perform the FFT
	perform_FFT();
	
	// compute first (N/2)+1 mag values
	for (int i = 0;i < (framesize/2)+1;i++)
	{
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
	}
	// mag spec symmetric above (N/2)+1 so copy previous values
	for (int i = (framesize/2)+1;i < framesize;i++)
	{
		mag[i] = mag[framesize-i];		
	}
	
	sum = 0;	// initialise sum to zero
	
	for (int i = 0;i < framesize;i++)
	{
		// calculate difference
		diff = mag[i] - mag_old[i];
		
		// only add up positive differences
		if (diff > 0)
		{
			// add difference to sum
			sum = sum+diff;
		}
		
		
		
		// store magnitude spectrum bin for next detection function sample calculation
		mag_old[i] = mag[i];
	}
	
	return sum;		
}


//=======================================================================
double OnsetDetectionFunction :: phase_deviation()
{
	double dev,pdev;
	double sum;
	
	// perform the FFT
	perform_FFT();
	
	sum = 0; // initialise sum to zero
	
	// compute phase values from fft output and sum deviations
	for (int i = 0;i < framesize;i++)
	{
		// calculate phase value
		phase[i] = atan2(out[i][1],out[i][0]);
		
		// calculate magnitude value
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
		
		
		// if bin is not just a low energy bin then examine phase deviation
		if (mag[i] > 0.1)
		{
			dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];	// phase deviation
			pdev = princarg(dev);	// wrap into [-pi,pi] range
		
			// make all values positive
			if (pdev < 0)	
			{
				pdev = pdev*-1;
			}
						
			// add to sum
			sum = sum + pdev;
		}
				
		// store values for next calculation
		phase_old_2[i] = phase_old[i];
		phase_old[i] = phase[i];
	}
	
	return sum;		
}

//=======================================================================
double OnsetDetectionFunction :: complex_spectral_difference()
{
	double dev,pdev;
	double sum;
	double mag_diff,phase_diff;
	double value;
	
	// perform the FFT
	perform_FFT();
	
	sum = 0; // initialise sum to zero
	
	// compute phase values from fft output and sum deviations
	for (int i = 0;i < framesize;i++)
	{
		// calculate phase value
		phase[i] = atan2(out[i][1],out[i][0]);
		
		// calculate magnitude value
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
		
		
		// phase deviation
		dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];	
		
		// wrap into [-pi,pi] range
		pdev = princarg(dev);	
		
		
		// calculate magnitude difference (real part of Euclidean distance between complex frames)
		mag_diff = mag[i] - mag_old[i];
		
		// calculate phase difference (imaginary part of Euclidean distance between complex frames)
		phase_diff = -mag[i]*sin(pdev);
		

		
		// square real and imaginary parts, sum and take square root
		value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
	
			
		// add to sum
		sum = sum + value;
		
		
		// store values for next calculation
		phase_old_2[i] = phase_old[i];
		phase_old[i] = phase[i];
		mag_old[i] = mag[i];
	}
	
	return sum;		
}

//=======================================================================
double OnsetDetectionFunction :: complex_spectral_difference_hwr()
{
	double dev,pdev;
	double sum;
	double mag_diff,phase_diff;
	double value;
	
	// perform the FFT
	perform_FFT();
	
	sum = 0; // initialise sum to zero
	
	// compute phase values from fft output and sum deviations
	for (int i = 0;i < framesize;i++)
	{
		// calculate phase value
		phase[i] = atan2(out[i][1],out[i][0]);
		
		// calculate magnitude value
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
		
		
		// phase deviation
		dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];	
		
		// wrap into [-pi,pi] range
		pdev = princarg(dev);	
		
		
		// calculate magnitude difference (real part of Euclidean distance between complex frames)
		mag_diff = mag[i] - mag_old[i];
		
		// if we have a positive change in magnitude, then include in sum, otherwise ignore (half-wave rectification)
		if (mag_diff > 0)
		{
			// calculate phase difference (imaginary part of Euclidean distance between complex frames)
			phase_diff = -mag[i]*sin(pdev);

			// square real and imaginary parts, sum and take square root
			value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
		
			// add to sum
			sum = sum + value;
		}
		
		// store values for next calculation
		phase_old_2[i] = phase_old[i];
		phase_old[i] = phase[i];
		mag_old[i] = mag[i];
	}
	
	return sum;		
}


//=======================================================================
double OnsetDetectionFunction :: high_frequency_content()
{
	double sum;
	
	// perform the FFT
	perform_FFT();
	
	sum = 0; // initialise sum to zero
	
	// compute phase values from fft output and sum deviations
	for (int i = 0;i < framesize;i++)
	{		
		// calculate magnitude value
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
		
		
		sum = sum + (mag[i]*((double) (i+1)));
		
		// store values for next calculation
		mag_old[i] = mag[i];
	}
	
	return sum;		
}

//=======================================================================
double OnsetDetectionFunction :: high_frequency_spectral_difference()
{
	double sum;
	double mag_diff;
	
	// perform the FFT
	perform_FFT();
	
	sum = 0; // initialise sum to zero
	
	// compute phase values from fft output and sum deviations
	for (int i = 0;i < framesize;i++)
	{		
		// calculate magnitude value
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
		
		// calculate difference
		mag_diff = mag[i] - mag_old[i];
		
		if (mag_diff < 0)
		{
			mag_diff = -mag_diff;
		}
		
		sum = sum + (mag_diff*((double) (i+1)));
		
		// store values for next calculation
		mag_old[i] = mag[i];
	}
	
	return sum;		
}

//=======================================================================
double OnsetDetectionFunction :: high_frequency_spectral_difference_hwr()
{
	double sum;
	double mag_diff;
	
	// perform the FFT
	perform_FFT();
	
	sum = 0; // initialise sum to zero
	
	// compute phase values from fft output and sum deviations
	for (int i = 0;i < framesize;i++)
	{		
		// calculate magnitude value
		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
		
		// calculate difference
		mag_diff = mag[i] - mag_old[i];
		
		if (mag_diff > 0)
		{
			sum = sum + (mag_diff*((double) (i+1)));
		}

		// store values for next calculation
		mag_old[i] = mag[i];
	}
	
	return sum;		
}


////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////// Methods to Calculate Windows ////////////////////////////////////

//=======================================================================
void OnsetDetectionFunction :: set_win_hanning()
{
	double N;		// variable to store framesize minus 1
	
	N = (double) (framesize-1);	// framesize minus 1
	
	// Hanning window calculation
	for (int n = 0;n < framesize;n++)
	{
		window[n] = 0.5*(1-cos(2*pi*(n/N)));
	}
}

//=======================================================================
void OnsetDetectionFunction :: set_win_hamming()
{
	double N;		// variable to store framesize minus 1
	double n_val;	// double version of index 'n'
	
	N = (double) (framesize-1);	// framesize minus 1
	n_val = 0;
	
	// Hamming window calculation
	for (int n = 0;n < framesize;n++)
	{
		window[n] = 0.54 - (0.46*cos(2*pi*(n_val/N)));
		n_val = n_val+1;
	}
}

//=======================================================================
void OnsetDetectionFunction :: set_win_blackman()
{
	double N;		// variable to store framesize minus 1
	double n_val;	// double version of index 'n'
	
	N = (double) (framesize-1);	// framesize minus 1
	n_val = 0;
	
	// Blackman window calculation
	for (int n = 0;n < framesize;n++)
	{
		window[n] = 0.42 - (0.5*cos(2*pi*(n_val/N))) + (0.08*cos(4*pi*(n_val/N)));
		n_val = n_val+1;
	}
}

//=======================================================================
void OnsetDetectionFunction :: set_win_tukey()
{
	double N;		// variable to store framesize minus 1
	double n_val;	// double version of index 'n'
	double alpha;	// alpha [default value = 0.5];
	
	alpha = 0.5;
	
	N = (double) (framesize-1);	// framesize minus 1
		
	// Tukey window calculation
	
	n_val = (double) (-1*((framesize/2)))+1;

	for (int n = 0;n < framesize;n++)	// left taper
	{
		if ((n_val >= 0) && (n_val <= (alpha*(N/2))))
		{
			window[n] = 1.0;
		}
		else if ((n_val <= 0) && (n_val >= (-1*alpha*(N/2))))
		{
			window[n] = 1.0;
		}
		else
		{
			window[n] = 0.5*(1+cos(pi*(((2*n_val)/(alpha*N))-1)));
		}

		n_val = n_val+1;			 
	}

}

//=======================================================================
void OnsetDetectionFunction :: set_win_rectangular()
{
	// Rectangular window calculation
	for (int n = 0;n < framesize;n++)
	{
		window[n] = 1.0;
	}
}



////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////// Other Handy Methods //////////////////////////////////////////

//=======================================================================
double OnsetDetectionFunction :: princarg(double phaseval)
{	
	// if phase value is less than or equal to -pi then add 2*pi
	while (phaseval <= (-pi)) 
	{
		phaseval = phaseval + (2*pi);
	}
	
	// if phase value is larger than pi, then subtract 2*pi
	while (phaseval > pi)
	{
		phaseval = phaseval - (2*pi);
	}
			
	return phaseval;
}