File content as of revision 0:6863633bf8a4:

#include "DFT.h"


// Preprocesses the data array and performs a dft on it, ARRAY MUST BE 1024 LONG
void performFFT(float* datam, int datalength){
	if (datalength == 1024){
		// - Calculating and substracting the mean  - Takes ~2N additions
		subMean(datam, datalength);
		// - Applying a hanning window - N multiplications
		applyHanningTable(datam, datalength);
		// - Re-indexing the array - 2N operations
		reindexData(datam, datalength);
		// It then applies a radix-2 Cooley Tukey FFT to create a DFT of the original signal
		cooleyTukeyRadix2(datam, datalength);
		// Remove the DC component and null everything above the nyquist
		afterProcessing(datam, datalength);
	}
	// Thus preprocessing takes 5N operations
	// FFT itself takes NlogN operations
	// Calculating to complex is N operations
	// Calulating the complex back to magnitude is ~4N operations
	// Total operation time is thus ~10N + NlogN complex operations
}

// Performs a radix-2 cooley tukey FFT on the data array
void cooleyTukeyRadix2(float* data, int datalength){
	
	complex_num fftdata[1024];
	complex_num fftbuffer[512];
	complex_num twiddlebuffer;

	complex_num minusone;
	minusone.real = -1;
	minusone.imaginary = 0;

	// Convert data to complex data
	for (int i = 0; i < datalength; i++){
		fftdata[i].real = data[i];
		fftdata[i].imaginary = 0;
		}

	for (int stage = 0; stage < 10; stage++){

		unsigned int step = 1<<(stage+1);

		for (unsigned int i = 0; i < 1024; i = i+step){ // Crack it up in parts, i.e. for stage0 -> 512, stage1 -> 256 , 2>128, 3>64, 4>32, 5>16, 6>8, 7>4, 8>2, 9>1 --> 1024/2^stage OR IN STEPS: stage0 -> 2,   stage1 -> 4,    2>8,   3>16, 4,32, 5>64, .........      9>512 --> 1024 - 1024/2^stage+1
			// Buffer the first 'half' of the thing, which is the first half step
			for (unsigned int k = 0; k < step / 2; k++){
				fftbuffer[k].real = fftdata[i+k].real;
				fftbuffer[k].imaginary = fftdata[i+k].imaginary;
			}
			
			// Apply twiddle factors to the seconds half (half step) depending on stage and numbering
			for (unsigned int k = 0; k < step / 2; k++){
				// Read in the correct twiddle factor
				// Twiddle factor is recalculated with k*(1024 / 2^stage+1) which can be recalculated to 1<<
				twiddlebuffer.imaginary = twiddleLookupTable_Imaginary[k * (1<<(9-stage))];
				twiddlebuffer.real = twiddleLookupTable_Real[k * (1 << (9 - stage))];
				
				fftdata[i + k + step / 2] = complex_multiply(fftdata[i + k + step / 2], twiddlebuffer);
			}
			
			// Recalculate the first half by addition of the first and recalculated half
			for (unsigned int k = 0; k < step / 2; k++)
				fftdata[i + k] = complex_add(fftdata[i+k], fftdata[i+k+step/2]);
			

			// Multiply the second half with minus 1
			for (unsigned int k = 0; k < step / 2; k++)
				fftdata[i + step / 2 + k] = complex_multiply(fftdata[i + step/2+k], minusone);

			// Add the buffered old values, with the new second half values
			for (unsigned int k = 0; k < step / 2; k++)
				fftdata[i + step / 2 + k] = complex_add(fftdata[i + step/2 + k], fftbuffer[k]);
		}
	}

	// Convert complex data back to magnitude data
	for (int i = 0; i < 1024; i++)
		data[i] = complex_magnitude(fftdata[i]);
}

// Calculates and subtracts the signal mean - Takes ~2N operations
void subMean(float* data, int datalength){
	float mean = 0;
	float total = 0;
	// Because the data is between 0 and 1, we can just sum and then device by datalength
	// Normally you would devide each data sample by the datalength and sum all of those (to prevent overflow), so this saves N times multiplications
	for (int i = 0; i < datalength; i++)
		total += data[i];

	mean = total / datalength;

	for (int i = 0; i < datalength; i++)
		data[i] = data[i] - mean;

}

// Reindexes the array in preperation of the FFT implemented using a LUT, 2N operations
void reindexData(float* data, int datalength){

	float reordereddata[1024];

	for (int i = 0; i < datalength; i++){
		reordereddata[reversebit_lut_i[i]] = data[i];
	}

	// Copy all the reorderddata to the original array (datalength*int size in bytes)
	memcpy(data, reordereddata, datalength*4);
}

// Removes the DC component and null everything above the nyquist frequency
void afterProcessing(float* data, int datalength){
	data[0] = 0; // Set DC component to 0

	for (int i = datalength-1; i > datalength / 2; i--) // Null everything above the nyquist frequency 
		data[i] = 0;

	// FIlter 50hz
	filterfreq(50, data, datalength);
}

// Filters a frequency and its multiples
void filterfreq(int frequency, float* data, int datalength){
	int niggr = 512;
	int freqperbin = 10;
	// Calculate which bin contains the frequency
	int bin = frequency / freqperbin;

	// Null bin frequencies and the two surrounding bins, plus its multiples adjecent
	do{
		if (bin - 1 >= 0) // Negative elements = you gonna have a bad time
			data[bin - 1] = 0;
		data[bin] = 0;
		data[bin + 1] = 0; // Since element 512-1024 also exist, we can null it even if it isnt really usefull
		bin *= 2;
	} while (bin<niggr);
}

// Calculate total magnitude, frequency magnitude component integrating module (FMCIM)
int calcPSD(float* data, int datalength){
	float knoedel=0;
	for (int i = 0; i < datalength; i++){
		knoedel += (float)data[i];
	}
	return (int)knoedel;
}