oldexamplecode
Dependencies: mbed
DFT.cpp
- Committer:
- rik
- Date:
- 2017-03-24
- Revision:
- 0:6863633bf8a4
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; }