Fast Fourier Transform (FFT) (Part 2)

Part 1234567891011

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

Next post on same subject

4 Comments

  1. MAF says:

    Best Arduino Blog!
    Really thanks for your knowledge and effort.

  2. Didier says:

    Nice to hear! Thanks.
    Didier

  3. Y3G says:

    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)

Leave a Reply

You must be logged in to post a comment.