d@0: /*
d@0:   ==============================================================================
d@0: 
d@0:     AudioSourceFeatureExtractor.cpp
d@0:     Created: 11 Aug 2014 10:41:02am
d@0:     Author:  david.ronan
d@0: 
d@0:   ==============================================================================
d@0: */
d@0: 
d@0: #include "AudioSourceFeatureExtractor.h"
d@0: #include "ObservationData.h"
d@0: #include <math.h>
d@0: #include <cmath>
d@0: #include <algorithm>
d@4: #include <vector>
d@7: #include <memory>
d@4: 
d@1: //#include <windows.h>
d@0: 
d@0: //==============================================================================  AudioSourceFeatureExtractor
d@15: AudioSourceFeatureExtractor::AudioSourceFeatureExtractor() 
d@0: {
d@4: 	m_fft = std::make_unique<FFTW>(FFTSIZE);
d@4: 	m_fInputBuffer = std::vector<float>(FFTSIZE, 0);;
d@4: 	m_fMagnitudeSpectrum = std::vector<float>(MAGNITUDESIZE, 0);
d@4: 	m_fPreviousMagnitudeSpectrum = std::vector<float>(MAGNITUDESIZE, 0);
d@4: 	m_OutReal = std::vector<float>(FFTSIZE, 0);
d@4: 	m_OutImag = std::vector<float>(FFTSIZE, 0);
d@0: }
d@0: 
d@0: //============================================================================== ~AudioSourceFeatureExtractor
d@0: AudioSourceFeatureExtractor::~AudioSourceFeatureExtractor()
d@0: {
d@0: 
d@0: }
d@0: 
d@0: //==============================================================================  init
d@0: void AudioSourceFeatureExtractor::Initialise(float fSampleRate)
d@1: {
d@0: 	m_fSampleRate = fSampleRate;
d@1: 	m_iWriteIdx = 0;
d@0: 
d@0: 	for (int i = 1; i <= MAGNITUDESIZE; i++)
d@0: 	{
d@0: 		float win = ((fSampleRate / (2 * MAGNITUDESIZE)) * i);
d@0: 		m_fFreqBins.push_back(win);
d@0: 	}
d@0: 
d@0: 	m_MFCC.initMFFCvariables(12, MAGNITUDESIZE, m_fSampleRate);
d@10: 	m_SpectralContrast.initSpectralContrastVariables(FFTSIZE, m_fSampleRate);
d@0: }
d@0: 
d@0: //==============================================================================  process
d@9: std::vector<ObservationData> AudioSourceFeatureExtractor::Process(const std::vector<float> data)
d@1: {
d@0: 	std::vector<ObservationData> observations = std::vector<ObservationData>();
d@0: 
d@4: 	m_fInputBuffer = std::vector<float>(FFTSIZE, 0);;
d@4: 	m_fMagnitudeSpectrum = std::vector<float>(MAGNITUDESIZE, 0);
d@4: 	m_fPreviousMagnitudeSpectrum = std::vector<float>(MAGNITUDESIZE, 0);
d@0: 	m_iWriteIdx = 0;
d@0: 
d@9: 	float perdiodicity = EstimatePerdiodicity(data);
d@0: 
d@9: 	float entropyofenergy = EntropyOfEnergy(data);
d@0: 
d@9: 	for (size_t n = 0; n < data.size(); n++)
d@0: 	{
d@0: 		m_fInputBuffer[m_iWriteIdx] = data[n];
d@0: 		m_iWriteIdx++;
d@0: 
d@0: 		//Do normal FFT
d@0: 
d@4: 		if (FFTSIZE == m_iWriteIdx)
d@0: 		{
d@0: 			//Get the peak amplitude rms and crest factor for the current frame
d@4: 			float peakvalue = 0;
d@4: 			float rms = 0;
d@12: 			float crestfactor;
d@0: 
d@4: 			for (int i = 0; i < FFTSIZE; i++)
d@0: 			{
d@4: 				if (m_fInputBuffer[i] > peakvalue)
d@0: 					peakvalue = abs(m_fInputBuffer[i]);
d@0: 
d@0: 				rms += (m_fInputBuffer[i] * m_fInputBuffer[i]);
d@0: 			}
d@1: 
d@0: 			rms /= FFTSIZE;
d@1: 			rms = sqrt(rms);
d@0: 
d@4: 			crestfactor = peakvalue / (rms + std::numeric_limits<float>::epsilon());
d@0: 
d@0: 			float ZCR = 0;
d@0: 			float inputBuffer2[FFTSIZE] = {0.0};
d@0: 
d@1: 			for (int i = 1; i < FFTSIZE; i++)
d@0: 			{
d@4: 				inputBuffer2[i] = m_fInputBuffer[i - 1];
d@0: 			}
d@0: 
d@1: 			for (int i = 0; i < FFTSIZE; i++)
d@0: 			{
d@0: 				ZCR += (abs(Sign(m_fInputBuffer[i]) - Sign(inputBuffer2[i])));
d@0: 			}
d@0: 
d@0: 			ZCR /= (2 * FFTSIZE);
d@0: 
d@0: 			//////////////////////////////////////////////////////////////////
d@1: 			for (int i = 0; i < FFTSIZE; i++)
d@0: 			{
d@4: 				float multiplier = static_cast<float>(0.5 * (1 - cos(2 * PI * i / (FFTSIZE - 1))));
d@0: 				m_fInputBuffer[i] = multiplier * m_fInputBuffer[i];
d@0: 			}
d@0: 
d@0: 			//Compute the FFT
d@4: 			m_OutReal = std::vector<float>(FFTSIZE, 0);
d@4: 			m_OutImag = std::vector<float>(FFTSIZE, 0);
d@4: 			
d@0: 
d@0: 			m_fft->process(m_fInputBuffer, m_OutReal, m_OutImag);
d@0: 
d@0: 			//Compute the Magnitude
d@0: 			VectorDistance(m_OutReal, 1, m_OutImag, 1, m_fMagnitudeSpectrum, 1, MAGNITUDESIZE);
d@0: 
d@0: 			//Compute the spectral features
d@4: 			float centroid = 0;
d@4: 			float spread = 0;
d@4: 			float skewness = 0;
d@4: 			float kurtosis = 0;
d@0: 			float brightness = 0;
d@0: 			float rolloff85 = 0;
d@0: 			float rolloff95 = 0;
d@0: 			float spectralentropy = 0;
d@0: 			float flatness = 0;
d@0: 			float spectralcf = 0;
d@1: 			float spectralflux = 0;
d@0: 			std::vector<float> mfccs;
d@10: 			std::vector<float> spectralContrast;
d@10: 			std::vector<float> spectralValleys;
d@0: 
d@1: 			SpectralFeatures(m_fMagnitudeSpectrum, m_fPreviousMagnitudeSpectrum, MAGNITUDESIZE, centroid, spread, skewness, kurtosis, brightness, rolloff85, rolloff95, spectralentropy, flatness, spectralcf, spectralflux);
d@1: 
d@0: 			m_MFCC.ComputeMFCC(m_fMagnitudeSpectrum, mfccs);
d@0: 
d@10: 			m_SpectralContrast.computeSpectralContrast(m_fMagnitudeSpectrum, spectralContrast, spectralValleys);
d@10: 
d@0: 			//Update for Spectral Flux
d@4: 			for (int i = 0; i < MAGNITUDESIZE; i++)
d@0: 			{
d@0: 				m_fPreviousMagnitudeSpectrum[i] = m_fMagnitudeSpectrum[i];
d@1: 			}
d@0: 
d@0: 			//50% Hop size
d@4: 			for (int i = MAGNITUDESIZE; i < FFTSIZE; i++)
d@0: 			{
d@0: 				m_fInputBuffer[i - MAGNITUDESIZE] = m_fInputBuffer[i];
d@0: 				m_iWriteIdx = MAGNITUDESIZE;
d@0: 			}
d@1: 
d@0: 			//create an observation for this window and push it to the list of observations for this 30 sec. audio clip
d@10: 			ObservationData newObservation = ObservationData(rms, peakvalue, crestfactor, ZCR, centroid, spread, skewness, kurtosis, brightness, rolloff85, rolloff95, spectralentropy, flatness, spectralcf, spectralflux, mfccs, spectralContrast, spectralValleys, perdiodicity, entropyofenergy);
d@1: 			observations.push_back(newObservation);
d@0: 		}
d@0: 	}
d@1: 
d@0: 	return observations;
d@0: }
d@0: 
d@0: void AudioSourceFeatureExtractor::Finalize()
d@0: {
d@0: }
d@0: 
d@4: void AudioSourceFeatureExtractor::VectorDistance(std::vector<float> vIn1, int stride1, std::vector<float> vIn2, int stride2, std::vector<float> &vOut, int strideOut, size_t nElements)
d@0: {
d@4: 	//Compute remaining samples
d@0: 	for (size_t i = 0; i < nElements; i++)
d@0: 	{
d@4: 		vOut[i * strideOut] = sqrt(vIn1[i * stride1] * vIn1[i * stride1] + vIn2[i * stride2] * vIn2[i * stride2]);
d@0: 	}
d@0: }
d@0: 
d@4: void AudioSourceFeatureExtractor::SpectralFeatures(std::vector<float> &magnitude, std::vector<float> &prevmagnitude, size_t windowSize, float& centroid, float& spread, float& skew, float& kurtosis, float& brightness, float& rolloff95, float& rolloff85, float& spectralentropy, float& flatness, float& spectralcf, float& spectralflux)
d@0: {
d@0: 	float summagbin = 0;
d@0: 	float summag = 0;
d@0: 	float sumprevmag = 0;
d@0: 	float summagspread = 0;
d@0: 	float summagskewness = 0;
d@0: 	float summagkurtosis = 0;
d@0: 	float sumbrightness = 0;
d@0: 	float maxmag = 0;
d@0: 
d@0: 	float geo = 1;
d@0: 	float totalenergy = 0;
d@0: 	float currEnergy95 = 0;
d@0: 	float currEnergy85 = 0;
d@0: 	int countFFT85 = 0;
d@0: 	int countFFT95 = 0;
d@0: 
d@0: 	float normalisedMagSumEntropy = 0;
d@0: 
d@0: 	//Centroid
d@0: 	//////////////////////////////////////////////////////
d@0: 
d@0: 	for (size_t i = 0; i < windowSize; i++)
d@0: 	{
d@1: 		//Moments
d@0: 		summagbin += (m_fFreqBins[i] * magnitude[i]);
d@0: 		summag += magnitude[i];
d@0: 		sumprevmag += prevmagnitude[i];
d@0: 
d@0: 		//Spectral Crest factor
d@4: 		if (magnitude[i] > maxmag)
d@0: 		{
d@0: 			maxmag = magnitude[i];
d@0: 		}
d@0: 
d@0: 		float energy = magnitude[i] * magnitude[i];
d@1: 
d@0: 		//Roll off + Flatness
d@1: 
d@0: 		totalenergy += energy;
d@4: 		float exponent = static_cast<float>(1.0 / static_cast<float>(windowSize));
d@1: 
d@4: 		if (energy != 0)
d@0: 		{
d@1: 			geo *= pow(energy, exponent);
d@0: 		}
d@0: 	}
d@0: 
d@0: 	summag += std::numeric_limits<float>::epsilon();
d@0: 
d@0: 	centroid = summagbin / summag;
d@0: 
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@0: 	//Spread
d@0: 	//////////////////////////////////////////////////////
d@0: 
d@0: 	for (size_t i = 0; i < windowSize; i++)
d@0: 	{
d@0: 		float norm = magnitude[i] / (summag + std::numeric_limits<float>::epsilon());
d@0: 		float normprev = prevmagnitude[i] / (sumprevmag + std::numeric_limits<float>::epsilon());
d@0: 
d@0: 		//Spectral Flux
d@0: 		spectralflux += ((norm - normprev) * (norm - normprev));
d@0: 
d@1: 		//Entropy
d@0: 		normalisedMagSumEntropy += (norm * log(norm + std::numeric_limits<float>::epsilon()));
d@0: 
d@0: 		//Moments
d@4: 		summagspread += pow((m_fFreqBins[i] - centroid), 2) * norm;
d@4: 		summagskewness += pow((m_fFreqBins[i] - centroid), 3) * norm;
d@4: 		summagkurtosis += pow((m_fFreqBins[i] - centroid), 4) * norm;
d@0: 
d@0: 		//Brightness
d@0: 		if (m_fFreqBins[i] >= 1500)
d@0: 		{
d@0: 			sumbrightness += magnitude[i];
d@0: 		}
d@0: 	}
d@0: 
d@0: 	spread = sqrt(summagspread);
d@0: 
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@0: 	//Skewness
d@0: 	//////////////////////////////////////////////////////
d@0: 
d@0: 	float spreadtemp = spread;
d@0: 	spreadtemp += std::numeric_limits<float>::epsilon();
d@0: 
d@4: 	skew = summagskewness / pow(spreadtemp, 3);
d@0: 
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@0: 	//Kurtosis
d@0: 	//////////////////////////////////////////////////////
d@0: 
d@0: 	kurtosis = summagkurtosis / pow(spreadtemp, 4);
d@0: 
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@0: 	//Brightness
d@0: 	//////////////////////////////////////////////////////
d@0: 
d@0: 	brightness = sumbrightness / summag;
d@0: 
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@0: 	//Roll off 85% - 95%
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@3: 	while ((currEnergy95 <= 0.95 * totalenergy) && (countFFT95 < static_cast<int>(windowSize)))
d@0: 	{
d@0: 		currEnergy95 += (magnitude[countFFT95] * magnitude[countFFT95]);
d@0: 		countFFT95 += 1;
d@0: 
d@0: 		if (currEnergy85 <= 0.85 * totalenergy)
d@0: 		{
d@0: 			currEnergy85 += (magnitude[countFFT85] * magnitude[countFFT85]);
d@0: 			countFFT85 += 1;
d@0: 		}
d@0: 	}
d@0: 
d@0: 	countFFT85 -= 1;
d@0: 	countFFT95 -= 1;
d@0: 
d@0: 	//Normalization
d@0: 	if (currEnergy85 <= 0)
d@0: 	{
d@0: 		//If we have silence then our value should be 0.
d@0: 		rolloff85 = 0;
d@0: 	}
d@0: 	else
d@0: 	{
d@3: 		rolloff85 = (static_cast<float>(countFFT85) / static_cast<float>(windowSize));
d@0: 	}
d@0: 
d@0: 	//Normalization
d@0: 	if (currEnergy95 <= 0)
d@0: 	{
d@0: 		//If we have silence then our value should be 0.
d@0: 		rolloff95 = 0;
d@0: 	}
d@0: 	else
d@0: 	{
d@3: 		rolloff95 = (static_cast<float>(countFFT95) / static_cast<float>(windowSize));
d@0: 	}
d@0: 
d@0: 	/////////////////////////////////////////////////////
d@0: 
d@0: 	//Spectral Entropy
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@3: 	spectralentropy = static_cast<float>((normalisedMagSumEntropy / log(static_cast<float>(windowSize))) * - 1);
d@0: 
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@0: 	/////////////////////////////////////////////////////
d@0: 
d@0: 	//Flatness
d@0: 	/////////////////////////////////////////////////////////
d@4: 	if (totalenergy != 0 && geo != 0)
d@0: 	{
d@3: 		flatness = geo / (totalenergy / static_cast<float>(windowSize));
d@0: 	}
d@0: 	else
d@0: 	{
d@0: 		flatness = 0;
d@0: 	}
d@0: 
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@0: 	//Spectral Crest Factor
d@0: 	/////////////////////////////////////////////////////////
d@0: 
d@4: 	if (summag != 0)
d@0: 	{
d@0: 		spectralcf = maxmag / summag;
d@0: 	}
d@0: 	else
d@0: 	{
d@0: 		spectralcf = 0;
d@0: 	}
d@0: }
d@0: 
d@8: void AudioSourceFeatureExtractor::XCorr(std::vector<float>& output, std::vector<float> input1, std::vector<float> input2, size_t maxLag = 0)
d@0: {
d@4: 	// TODO: adapt for different sizes
d@0: 
d@4: 	// find next power of 2 for efficient FFT
d@4: 	size_t fftSize = 1;
d@8: 	while (fftSize < 2 * static_cast<size_t>(input1.size()) - 1)
d@4: 	{
d@4: 		fftSize *= 2;
d@4: 	}
d@1: 
d@0: 	FFTW fft1 = FFTW(fftSize);
d@0: 	FFTW fft2 = FFTW(fftSize);
d@1: 
d@4: 	// allocate space for FFTW processing
d@4: 	std::vector<float> in1 = std::vector<float>(fftSize, 0);;
d@4: 	std::vector<float> in2 = std::vector<float>(fftSize, 0);;
d@0: 
d@4: 	// output from FFT is only up to Nyquist frequency
d@4: 	std::vector<float> out1Real = std::vector<float>(fftSize, 0);
d@4: 	std::vector<float> out1Imag = std::vector<float>(fftSize, 0);
d@4: 	std::vector<float> out2Real = std::vector<float>(fftSize, 0);
d@4: 	std::vector<float> out2Imag = std::vector<float>(fftSize, 0);
d@0: 
d@4: 	// copy input to allocated arrays
d@8: 	for (size_t i = 0; i < static_cast<size_t>(input1.size()); ++i)
d@4: 	{
d@4: 		in1[i] = input1[i];
d@4: 		in2[i] = input2[i];
d@4: 	}
d@0: 
d@0: 	fft1.process(in1, out1Real, out1Imag);
d@0: 	fft2.process(in2, out2Real, out2Imag);
d@0: 
d@4: 	// multiply out1*conj(out2) in FFT domain
d@4: 	fftwf_complex* product = static_cast<fftwf_complex*>(fftwf_malloc(sizeof(fftwf_complex) * fftSize));
d@4: 	for (size_t i = 0; i < fftSize; ++i)
d@4: 	{
d@4: 		product[i][0] = (out1Real[i] * out2Real[i] + out1Imag[i] * out2Imag[i]);
d@4: 		product[i][1] = (out1Imag[i] * out2Real[i] - out1Real[i] * out2Imag[i]);
d@4: 	}
d@0: 
d@4: 	// compute IFFT
d@0: 
d@0: 	// do FFT, adapted from iTraktorLib "FFTW.h"
d@4: 	fftw_iodim dim;
d@4: 	dim.n = fftSize;
d@4: 	dim.is = 1;
d@4: 	dim.os = 1;
d@0: 
d@4: 	fftwf_plan plan;
d@0: 
d@4: 	float* result = static_cast<float*>(fftwf_malloc(sizeof(float) * fftSize));
d@4: 	plan = fftwf_plan_guru_dft_c2r(1, &dim, 0, nullptr, product, result, FFTW_ESTIMATE);
d@4: 	fftwf_execute_dft_c2r(plan, product, result);
d@4: 	fftwf_destroy_plan(plan);
d@0: 
d@4: 	// copy real part of result back to output array, normalize by FFT size
d@8: 	if (maxLag == 0 || maxLag >= static_cast<size_t>(input1.size()))
d@4: 	{
d@8: 		maxLag = static_cast<size_t>(input1.size()) - 1;
d@4: 	}
d@0: 
d@8: 	output = std::vector<float>(2 * maxLag + 1, 0);
d@0: 
d@4: 	for (size_t i = fftSize - maxLag; i < fftSize; ++i)
d@4: 	{
d@4: 		output[i - fftSize + maxLag] = result[i] / fftSize;
d@4: 	}
d@0: 
d@4: 	for (unsigned i = 0; i <= static_cast<unsigned>(maxLag); ++i)
d@4: 	{
d@4: 		output[i + maxLag] = result[i] / fftSize;
d@4: 	}
d@0: 
d@0: 	fftwf_free(result);
d@0: 	fftwf_free(product);
d@0: }
d@0: 
d@9: float AudioSourceFeatureExtractor::EstimatePerdiodicity(std::vector<float> data)
d@0: {
d@12: 	float periodicity;
d@8: 
d@9: 	std::vector<float> downsampleddata = std::vector<float>(data.size(), 0);
d@8: 
d@8: 	std::vector<float> xcorr;
d@0: 
d@9: 	for (size_t k = 0; k < data.size(); k++)
d@0: 	{
d@0: 		downsampleddata[k] = data[k];
d@0: 	}
d@8: 	
d@8: 	EnvelopeCurve(downsampleddata, downsampleddata, m_fSampleRate);
d@0: 
d@3: 	int startLag = static_cast<int>(floor(static_cast<float>(60.0 / 480.0 * m_fSampleRate)));
d@3: 	int endLag = static_cast<int>(floor(static_cast<float>(60.0 / 30.0 * m_fSampleRate)));
d@8: 	int thresh = static_cast<int>(floor(static_cast<float>(0.2 * m_fSampleRate)));
d@0: 
d@12: 	
d@8: 	XCorr(xcorr, downsampleddata, downsampleddata, endLag);
d@12: 	size_t xcorrlen = xcorr.size();
d@0: 
d@3: 	int i = static_cast<int>(floor(xcorrlen / 2));
d@0: 
d@0: 	std::vector<float> temp;
d@4: 	for (size_t j = i; j < xcorrlen; j++)
d@0: 	{
d@0: 		temp.push_back(xcorr[j]);
d@0: 	}
d@0: 
d@0: 	std::vector<float> temp2;
d@4: 	for (size_t j = static_cast<size_t>(std::max(1, startLag)); j < static_cast<size_t>(std::min(endLag, static_cast<int>(temp.size()))); j++)
d@0: 	{
d@0: 		temp2.push_back(temp[j]);
d@0: 	}
d@0: 
d@0: 	int maxIdx = 0;
d@0: 	float maxVal = std::numeric_limits<float>::lowest();
d@4: 	for (size_t j = 0; j < temp2.size(); j++)
d@0: 	{
d@4: 		if (temp2[j] > maxVal)
d@0: 		{
d@0: 			maxIdx = j;
d@0: 			maxVal = temp2[j];
d@0: 		}
d@0: 	}
d@0: 
d@4: 	int from = std::max(1, maxIdx - thresh);
d@4: 	int to = std::min(static_cast<int>(temp2.size()), maxIdx + thresh);
d@0: 
d@1: 	float minValInRange = std::numeric_limits<float>::max();
d@0: 
d@9: 	for (size_t j = static_cast<size_t>(from); j < static_cast<size_t>(to); j++)
d@1: 	{
d@4: 		if (temp2[j] < minValInRange)
d@0: 		{
d@0: 			minValInRange = temp2[j];
d@0: 		}
d@0: 	}
d@1: 
d@4: 	periodicity = (maxVal - minValInRange);
d@0: 
d@12: 	//std::string dave = "Periodicity " + std::to_string(periodicity) + "\n";
d@1: 	//OutputDebugString(dave.c_str());
d@0: 
d@0: 	return periodicity;
d@0: }
d@0: 
d@8: void AudioSourceFeatureExtractor::EnvelopeCurve(std::vector<float> data, std::vector<float> & out, float sampleRate)
d@0: {
d@4: 	float release = 0.02f;
d@0: 	float envelope = 0.0f;
d@0: 
d@1: 	float gr = exp(-1 / (sampleRate * release));
d@0: 
d@8: 	for (size_t i = 0; i < data.size(); i++)
d@1: 	{
d@0: 		float EnvIn = abs(data[i]);
d@0: 
d@4: 		if (envelope < EnvIn)
d@0: 		{
d@0: 			envelope = EnvIn;
d@0: 		}
d@0: 		else
d@0: 		{
d@0: 			envelope = envelope * gr;
d@0: 			envelope = envelope + (1 - gr) * EnvIn;
d@1: 		}
d@1: 
d@1: 		out[i] = envelope;
d@0: 	}
d@0: }
d@0: 
d@9: float AudioSourceFeatureExtractor::EntropyOfEnergy(std::vector<float> data)
d@0: {
d@0: 	float totalEnergy = 0.0;
d@0: 	float totalentropy = 0.0;
d@0: 	const int numSubFrames = 30;
d@0: 
d@0: 	//Calculate Total Energy and normalise it
d@9: 	for (size_t i = 0; i < data.size(); i++)
d@0: 	{
d@0: 		totalEnergy += (data[i] * data[i]);
d@0: 	}
d@0: 
d@9: 	int subFrameSize = data.size() / numSubFrames;
d@0: 
d@0: 	for (size_t i = 0; i < numSubFrames; i++)
d@0: 	{
d@0: 		float val = 0.0;
d@0: 
d@3: 		for (size_t j = 0; j < static_cast<size_t>(subFrameSize); j++)
d@0: 		{
d@0: 			val += (data[j + (i * subFrameSize)] * data[j + (i * subFrameSize)]);
d@0: 		}
d@0: 
d@0: 		float frameval = val / (totalEnergy + std::numeric_limits<float>::epsilon());
d@0: 
d@0: 		totalentropy += (frameval * Log2(frameval + std::numeric_limits<float>::epsilon()));
d@0: 	}
d@0: 
d@0: 	return (totalentropy * -1);
d@0: }
d@0: 
d@1: float AudioSourceFeatureExtractor::Log2(float n)
d@1: {
d@1: 	// log(n)/log(2) is log2.
d@4: 	return logf(n) / logf(2);
d@0: }
d@0: 
d@0: int AudioSourceFeatureExtractor::Sign(float x)
d@0: {
d@0: 	int retValue = 0;
d@0: 	if (x >= 0)
d@0: 		retValue = 1;
d@0: 	if (x < 0)
d@0: 		retValue = -1;
d@0: 
d@0: 	return retValue;
d@0: }
d@0: 
d@0: std::vector<float> AudioSourceFeatureExtractor::LinSpace(float min, float max, int n)
d@0: {
d@0: 	std::vector<float> result;
d@0: 	// vector iterator
d@0: 	int iterator = 0;
d@0: 
d@4: 	for (int i = 0; i <= n - 2; i++)
d@0: 	{
d@12: 		float temp = min + i * (max - min) / (floor(static_cast<float>(n)) - 1);
d@0: 		result.insert(result.begin() + iterator, temp);
d@0: 		iterator += 1;
d@0: 	}
d@0: 
d@0: 	result.insert(result.begin() + iterator, max);
d@0: 	return result;
d@0: }
d@0: 
d@0: std::vector<float> AudioSourceFeatureExtractor::LogSpace(float min, float max, int n)
d@0: {
d@0: 	float logMin = log(min);
d@0: 	float logMax = log(max);
d@0: 	double delta = (logMax - logMin) / n;
d@0: 	double accDelta = 0;
d@0: 	std::vector<float> v;
d@0: 
d@0: 	for (int i = 0; i <= n; i++)
d@0: 	{
d@4: 		v.push_back(static_cast<float>(exp(logMin + accDelta)));
d@0: 		accDelta += delta;// accDelta = delta * i
d@0: 	}
d@0: 
d@0: 	return v;
d@0: }
d@8: 
d@8: //void AudioSourceFeatureExtractor::ConvolveFunction(float* & z, float* x, float* y, size_t& lenz, size_t lenx, size_t leny)
d@8: //{
d@8: //	// computes cross-correlation of two vectors, takes 2 vectors of the same length and computes 2*length-1 lags
d@8: //
d@8: //	//float *zptr, s, *xp, *yp;
d@8: //	//size_t n_lo, n_hi;
d@8: //
d@8: //	//lenz = lenx + leny - 1;
d@8: //
d@8: //	//z = new float[lenz];
d@8: //
d@8: //	//zptr = z;
d@8: //
d@8: //	//for (size_t i = 0; i < lenz; i++)
d@8: //	//{
d@8: //	//	s = 0.0;
d@8: //	//	n_lo = 0 > (i - leny + 1) ? 0 : (i - leny + 1);
d@8: //	//	n_hi = (lenx - 1) < i ? (lenx - 1) : i;
d@8: //	//	xp = &x[0] + n_lo;
d@8: //	//	yp = &y[0] + i - n_lo;
d@8: //
d@8: //	//	for (size_t n = n_lo; n <= n_hi; n++)
d@8: //	//	{
d@8: //	//		s += *xp * *yp;
d@8: //	//		xp++;
d@8: //	//		yp--;
d@8: //	//	}
d@8: //
d@8: //	//	*zptr = s;
d@8: //	//	zptr++;
d@8: //	//}
d@8: //}
d@8: 
d@8: //void AudioSourceFeatureExtractor::DownSampler(float* data, float* & out, size_t lenIn, size_t& lenOut, float currentSampleRate, float futureSampleRate)
d@8: //{
d@8: //	////Low pass our data before we decimate
d@8: //	//const int filterlength = 101;
d@8: //	//float* tempdata = new float[lenIn];
d@8: //	//memset(tempdata, 0, sizeof(float)*lenIn);
d@8: //
d@8: //	////Coefficients were taken from Matlab
d@8: //	//float coeffs[filterlength] ={2.0839274e-05f, 6.6880953e-06f,-4.7252855e-05f,	-3.4874138e-05f,	7.8508412e-05f,9.7089069e-05f,-9.9197161e-05f,	-0.00020412984f,	8.2261082e-05f,	0.00035689832f,	9.4890293e-06f,	-0.00053820928f,	-0.00021722882f,	0.00070567406f,	0.00057491515f,	-0.00078844035f, -0.0010939316f,	0.00069057313f,	0.0017460365f,	-0.00030308162f,	-0.0024484061f,	-0.00047496572f,	0.0030552838f,	0.0017078921f, 	-0.0033604759f, 	-0.0033912712f, 	0.0031135236f, 	0.0054217125f, 	-0.0020499325f, 	-0.0075746416f, 	-6.7239242e-05f, 0.0094950572f, 0.0034002098f,	-0.010703080f, 	-0.0079918290f, 	0.010610434f, 0.013732497f, 	-0.0085344436f, 	-0.020346783f, 0.0036776769f,	0.027405130f,	0.0050050775f,	-0.034361813f,	-0.019287759f,	0.040615436f,	0.043452863f,-0.045583628f,	-0.093336113f,	0.048780452f,	0.31394306f,	0.45011634f,	0.31394306f,	0.048780452f,-0.093336113f,-0.045583628f,	0.043452863f,0.040615436f,-0.019287759f,	-0.034361813f,	0.0050050775f,	0.027405130f,0.0036776769f,-0.020346783f,	-0.0085344436f,	0.013732497f,	0.010610434f,	-0.0079918290f,	-0.010703080f,	0.0034002098f,	0.0094950572f,	-6.7239242e-05f,	-0.0075746416f,	-0.0020499325f,	0.0054217125f,	0.0031135236f,	-0.0033912712f,	-0.0033604759f,	0.0017078921f,	0.0030552838f,	-0.00047496572f,	-0.0024484061f,	-0.00030308162f,	0.0017460365f,0.00069057313f,	-0.0010939316f,	-0.00078844035f,	0.00057491515f,	0.00070567406f,	-0.00021722882f,	-0.00053820928f,	9.4890293e-06f,	0.00035689832f,	8.2261082e-05f,	-0.00020412984f,	-9.9197161e-05f,	9.7089069e-05f,	7.8508412e-05f,	-3.4874138e-05f,	-4.7252855e-05f,	6.6880953e-06f,	2.0839274e-05f};
d@8: //	//
d@8: //	//float acc;     // accumulator for MACs
d@8: //	//   float* coeffp; // pointer to coefficients
d@8: //	//   float* inputp; // pointer to input samples
d@8: //	//float* endp = &data[lenIn -1]; // pointer to last sample
d@8: //	//
d@8: //	//
d@8: //	//   // apply the filter to each input sample
d@8: //	//   for (size_t n = 0; n < lenIn; n++ )
d@8: //	//{
d@8: //	//       // calculate output n
d@8: //	//       coeffp = &coeffs[0];
d@8: //	//       inputp = &data[filterlength - 1 + n];
d@8: //
d@8: //	//	acc = 0.0f;
d@8: //	//	for (int k = 0; k < filterlength; k++ )
d@8: //	//	{
d@8: //	//		if(inputp <= endp)
d@8: //	//			acc += (*coeffp++) * (*inputp--);
d@8: //	//		else
d@8: //	//		{
d@8: //	//			//When we reach the end of the buffer
d@8: //	//			acc += (*coeffp++) * 0.0f;
d@8: //	//			*inputp--;
d@8: //	//		}
d@8: //	//	}
d@8: //
d@8: //	//	tempdata[n] = acc;
d@8: //	//   }
d@8: //	//int ratio = (int) (currentSampleRate / futureSampleRate);
d@8: //	//
d@8: //	//lenOut = lenIn / ratio;
d@8: //
d@8: //	//out = new float[lenOut];
d@8: //	//memset(out, 0, sizeof(float)*lenOut);
d@8: //
d@8: //
d@8: //	//int idx = 0;
d@8: //	//for(size_t i = 0; i < lenIn; i = i + ratio)
d@8: //	//{
d@8: //	//	out[idx] = tempdata[i];
d@8: //	//	idx++;
d@8: //	//}
d@8: //
d@8: //	//if(tempdata != NULL)
d@8: //	//{
d@8: //	//	delete[] tempdata;
d@8: //	//	tempdata = nullptr;
d@8: //	//}
d@8: //}
d@8: 
d@8: //void AudioSourceFeatureExtractor::Normalise(float* data, float* & out, size_t len)
d@8: //{
d@8: //	//float maxvalue = std::numeric_limits<float>::lowest();
d@8: //
d@8: //	//for (size_t i = 0; i < len; i++)
d@8: //	//{
d@8: //	//	if (data[i] > maxvalue)
d@8: //	//		maxvalue = data[i];
d@8: //	//}
d@8: //
d@8: //	//for (size_t i = 0; i < len; i++)
d@8: //	//{
d@8: //	//	out[i] = data[i] / maxvalue;
d@8: //	//}
d@8: //}
d@8: 
d@8: //void AudioSourceFeatureExtractor::PDF_getResampleKrnl(std::vector<float> freqVec, float fsProc, int nfft, int nBin, std::vector<float>& outLogFreqVec, std::vector<std::vector<float>>& outKrnl)
d@8: //{
d@8: //	//float f1 = freqVec[0];
d@8: //	//float f2 = freqVec[freqVec.size() - 1];
d@8: //
d@8: //	//outLogFreqVec = LinSpace(log(f1), log(f2), nBin);
d@8: //	//std::vector<float> f(nfft / 2, 0.0f);
d@8: //
d@8: //	//for (size_t i = 0; i < static_cast<size_t>(nBin); i++)
d@8: //	//{
d@8: //	//	outKrnl.push_back(f);
d@8: //	//}
d@8: //
d@8: //	//for (size_t i = 0; i < outLogFreqVec.size(); i++)
d@8: //	//{
d@8: //	//	outLogFreqVec[i] = exp(outLogFreqVec[i]);
d@8: //	//}
d@8: //
d@8: //	//for (size_t i = 1; i < static_cast<size_t>(nBin) - 2; i++)
d@8: //	//{
d@8: //	//	float freqBinLin = outLogFreqVec[i] / (fsProc / nfft);
d@8: //
d@8: //	//	float nextBinUp = outLogFreqVec[i + 1] / (fsProc / nfft);
d@8: //
d@8: //	//	float nextBinDown = outLogFreqVec[i - 1] / (fsProc / nfft);
d@8: //
d@8: //	//	float filtWid = nextBinUp - nextBinDown;
d@8: //
d@8: //	//	if (filtWid <= 2)
d@8: //	//	{
d@8: //	//		// log resolution is finer than linear resolution
d@8: //	//		// linear interpolation  of neighboring bins
d@8: //
d@8: //	//		float binDown = floor(freqBinLin);
d@8: //	//		float frac = freqBinLin - binDown;
d@8: //	//		std::vector<float> filt(nfft / 2, 0.0f);
d@8: //	//		filt[((int)binDown) - 1] = 1 - frac;
d@8: //	//		filt[((int)binDown + 1) - 1] = frac;
d@8: //	//		outKrnl[i][((int)binDown) - 1] = filt[((int)binDown) - 1];
d@8: //	//		outKrnl[i][((int)binDown + 1) - 1] = filt[((int)binDown + 1) - 1];
d@8: //	//	}
d@8: //	//	else
d@8: //	//	{
d@8: //	//		float ceilNextBinDown = ceil(nextBinDown);
d@8: //	//		float floorFreqBinLin = floor(freqBinLin);
d@8: //	//		float ceilFreqBinLin = ceil(freqBinLin);
d@8: //	//		float floorNextBinUp = floor(nextBinUp);
d@8: //
d@8: //	//		std::vector<int> idxLow;
d@8: //	//		std::vector<int> idxHigh;
d@8: //	//		std::vector<int> filtIdx;
d@8: //
d@8: //	//		for (size_t j = (size_t)ceilNextBinDown; j <= (size_t)floorFreqBinLin; j++)
d@8: //	//		{
d@8: //	//			idxLow.push_back(j);
d@8: //	//		}
d@8: //
d@8: //	//		for (size_t j = (size_t)ceilFreqBinLin; j <= (size_t)floorNextBinUp; j++)
d@8: //	//		{
d@8: //	//			idxHigh.push_back(j);
d@8: //	//		}
d@8: //
d@8: //	//		std::vector<int>::iterator it;
d@8: //	//		it = std::unique(idxLow.begin(), idxLow.end());
d@8: //	//		idxLow.resize(std::distance(idxLow.begin(), it));
d@8: //
d@8: //	//		it = std::unique(idxHigh.begin(), idxHigh.end());
d@8: //	//		idxHigh.resize(std::distance(idxHigh.begin(), it));
d@8: //
d@8: //	//		filtIdx.reserve(idxLow.size() + idxHigh.size()); // preallocate memory
d@8: //	//		filtIdx.insert(filtIdx.end(), idxLow.begin(), idxLow.end());
d@8: //	//		filtIdx.insert(filtIdx.end(), idxHigh.begin(), idxHigh.end());
d@8: //
d@8: //	//		std::vector<float> rampUp;
d@8: //	//		std::vector<float> rampDown;
d@8: //	//		for (size_t j = 0; j < filtIdx.size(); j++)
d@8: //	//		{
d@8: //	//			rampUp.push_back(1 - (freqBinLin - filtIdx[j]) / (freqBinLin - nextBinDown));
d@8: //	//			rampDown.push_back(1 - (filtIdx[j] - freqBinLin) / (nextBinUp - freqBinLin));
d@8: //	//		}
d@8: //
d@8: //	//		std::vector<float> filt;
d@8: //
d@8: //	//		for (size_t j = 0; j < rampUp.size(); j++)
d@8: //	//		{
d@8: //	//			filt.push_back(std::min(rampUp[j], rampDown[j]));
d@8: //	//		}
d@8: //
d@8: //	//		float sumfilt = 0.0f;
d@8: //	//		for (size_t j = 0; j < filt.size(); j++)
d@8: //	//		{
d@8: //	//			sumfilt += filt[j];
d@8: //	//		}
d@8: //
d@8: //	//		for (size_t j = 0; j < filt.size(); j++)
d@8: //	//		{
d@8: //	//			filt[j] /= sumfilt;
d@8: //	//		}
d@8: //
d@8: //	//		for (size_t j = 0; j < filtIdx.size(); j++)
d@8: //	//		{
d@8: //	//			outKrnl[i][filtIdx[j] - 1] = filt[j];
d@8: //	//		}
d@8: //	//	}
d@8: //	//}
d@8: //
d@8: //	// special routine for first bin
d@8: //	// get frequency in linear bins
d@8: //
d@8: //	//float freqBinLin = outLogFreqVec[0] / (fsProc / (float)nfft);
d@8: //	//float binDown = floor(freqBinLin);
d@8: //	//float frac = freqBinLin - binDown;
d@8: //
d@8: //	//std::vector<float> filt(nfft / 2, 0.0f);
d@8: //	//filt[((int)binDown) - 1] = 1 - frac;
d@8: //	//filt[((int)binDown + 1) - 1] = frac;
d@8: //
d@8: //	//outKrnl[0][((int)binDown) - 1] = filt[((int)binDown) - 1];
d@8: //	//outKrnl[0][((int)binDown + 1) - 1] = filt[((int)binDown + 1) - 1];
d@8: //
d@8: //	//// special routine for last bin
d@8: //	//// get frequency in linear bins
d@8: //	//freqBinLin = outLogFreqVec[outLogFreqVec.size() - 1] / (fsProc / nfft);
d@8: //
d@8: //	//// get next lower bin
d@8: //	//float nextBinDown = outLogFreqVec[outLogFreqVec.size() - 2] / (fsProc / nfft);
d@8: //
d@8: //	//float ceilNextBinDown = ceil(nextBinDown); //Subtract by -1 because of array
d@8: //	//float floorFreqBinLin = floor(freqBinLin);
d@8: //
d@8: //	//std::vector<int> filtIdx;
d@8: //
d@8: //	//for (size_t i = (size_t)ceilNextBinDown; i <= (size_t)floorFreqBinLin; i++)
d@8: //	//{
d@8: //	//	filtIdx.push_back(i);
d@8: //	//}
d@8: //
d@8: //	//std::vector<float> filt2;
d@8: //	//for (size_t i = 0; i < filtIdx.size(); i++)
d@8: //	//{
d@8: //	//	filt2.push_back(1 - (freqBinLin - filtIdx[i]) / (freqBinLin - nextBinDown));
d@8: //	//}
d@8: //
d@8: //	//float sumfilt = 0.0f;
d@8: //	//for (size_t i = 0; i < filt2.size(); i++)
d@8: //	//{
d@8: //	//	sumfilt += filt2[i];
d@8: //	//}
d@8: //
d@8: //	//for (size_t i = 0; i < filt2.size(); i++)
d@8: //	//{
d@8: //	//	filt2[i] /= sumfilt;
d@8: //	//}
d@8: //
d@8: //	//for (size_t j = 0; j < filtIdx.size(); j++)
d@8: //	//{
d@8: //	//	outKrnl[outKrnl.size() - 1][filtIdx[j] - 1] = filt2[j];
d@8: //	//}
d@8: //}
d@8: 
d@8: