Fast Fourier Transform (FFT) (Part 2)
Part 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Here is a recap of the pieces of code that we need in order to convert a wave into a frequency spectrum:
- Store data in a vector (Double data type)
- Weigh this data according to one function (default is Rectangle, also known as box car, in other words, no weighing!)
- Execute the FFT algorithm
- Convert complex values in usable data
Without knowing too much of the details from the various algorithms, you will very quickly face a dilemma: would you prefer speed or precision? Precision will require large vectors of data, while speed will require multiple vectors: e.g. instead of computing weighed values for each new set of data, we could compute weighing factors once and apply them repeatedly. Same thoughts for the bit reversal at the top of the FFT algorithm.
The proposed example do not use pre-processing of data, to the benefit of the number of samples, and to the clarity of the code.
Firstly, some constants shall be declared in the header (.h) file
// Custom constants #define FORWARD 0x01 #define REVERSE 0x00 // Windowing type #define WIN_TYP_RECTANGLE 0x00// rectangle (Box car) #define WIN_TYP_HAMMING 0x01// hamming #define WIN_TYP_HANN 0x02// hann #define WIN_TYP_TRIANGLE 0x03// triangle (Bartlett) #define WIN_TYP_BLACKMAN 0x04// blackmann #define WIN_TYP_FLT_TOP 0x05// flat top #define WIN_TYP_WELCH 0x06// welch
Then data shall be stored in one of the two fixed size vectors
const uint8_t samples = 64; double vReal[samples]; double vImag[samples];
And this is the weighing routine
void PlainFFT::windowing(double *vData, uint8_t samples, uint8_t windowType, uint8_t dir) {
// The weighing function is symetric; half the weighs are recorded
double samplesMinusOne = (double(samples) - 1.0);
for (uint8_t i = 0; i < (samples >> 1); i++) {
double indexMinusOne = double(i);
double ratio = (indexMinusOne / samplesMinusOne);
double weighingFactor = 1.0;
// compute and record weighting factor
switch (windowType) {
case WIN_TYP_RECTANGLE: // rectangle (box car)
weighingFactor = 1.0;
break;
case WIN_TYP_HAMMING: // hamming
weighingFactor = 0.54 - (0.46 * cos(2.0 * pi * ratio));
break;
case WIN_TYP_HANN: // hann
weighingFactor = 0.54 * (1.0 - cos(2.0 * pi * ratio));
break;
case WIN_TYP_TRIANGLE: // triangle (Bartlett)
weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne);
break;
case WIN_TYP_BLACKMAN: // blackmann
weighingFactor = 0.42323 - (0.49755 * (cos(2.0 * pi * ratio))) + (0.07922 * (cos(4.0 * pi * ratio)));
break;
case WIN_TYP_FLT_TOP: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(2.0 * pi * ratio)) + (0.1980399 * cos(4.0 * pi * ratio));
break;
case WIN_TYP_WELCH: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0));
break;
}
if (dir == FORWARD) {
vData[i] *= weighingFactor;
vData[samples - (i + 1)] *= weighingFactor;
}
else {
vData[i] /= weighingFactor;
vData[samples - (i + 1)] /= weighingFactor;
}
}
}
Notes:
- There is a little trick here. As the weighing function is symetrical, why bother computing them all? Half of them are computed and applied symetrically to the vector of data
- The dir parameter stands for direction: remember, FFT is reversible, so that we can apply it in FORWARD or REVERSE mode
Best Arduino Blog!
Really thanks for your knowledge and effort.
Nice to hear! Thanks.
Didier
Wow….
25 years I used C-language….
WHat name should the HaederFile be used? I used “PlainFFT.h” and then put in the sketch the data and the weighting routine, but keep receiving an error, although both files are located in the same folder….
Sure, I am doing some trivial stuff very wrong.
Thanks for any hints
Y3G (…when there is already 4G)
Sorry, I could not find my crystal ball 😉 . Would you please send a copy of the error message that you got.
Thanks