Diff: DFT.cpp

Revision:

0:6863633bf8a4

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DFT.cpp	Fri Mar 24 11:22:30 2017 +0000
@@ -0,0 +1,148 @@
+
+#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;
+}
\ No newline at end of file