Fourier Arduino 5m2q34

Conceptualmente hablando FFT es bastante simple. En la práctica, la FFT es difícil! Así que no voy a entrar en todos los detalles de este algoritmo y voy a permanecer centrado en su uso y aplicaciones. Esta serie de mensajes va a terminar con una biblioteca con todas las funciones que usted será capaz de utilizar para todo tipo de aplicaciones, sin preocuparse de las matemáticas. La trama en un osciloscopio muestra una onda que está hecha de muestras que se caracterizan por sus intensidades (ordenadas) y su tiempo de muestreo (abscisa). Estamos en el dominio del tiempo.

Cuanto más simple sea la señal (por ejemplo, una onda sinusoidal simple), el más simple es la caracterización. La transformada de Fourier propone para descomponer cualquier señal en una suma de seno y coseno. Para cada punto de datos del espectro de potencia de FFT corresponde una magnitud (ordenadas) y una frecuencia (abscisa). Ahora estamos en el dominio de la frecuencia.

Cuanto más compleja es la señal (por ejemplo, la señal + ruido + armónicos + +) ...

..... más compleja es la caracterización.

Y ¿adivinen qué? El algoritmo de la FFT se puede ejecutar en modo inverso, de modo que a partir de una FFT puede reconstruir una señal de la vida real. Podemos estudiar más adelante cómo esta propiedad puede ser utilizado para la de eliminación de ruido de las señales. FFT se aplica a vectores que contienen las muestras n, donde n debe ser una potencia de 2. Ejecutar el algoritmo de estos datos puede conducir a resultados inesperados. La razón es que el vector trunca la información de la señal y pueden contener ondas incompletamente descritas (por ejemplo,

las ondas de baja frecuencia). Este efecto secundario puede ser corregida por medio de pesaje de la señal, dando menos importancia a los datos principales y colas. Esta es la función de ventanas. La función de pesaje en la ventanas depende del tipo de señal a analizar: Transients whose duration is shorter than the length of the window : Rectangular (Box car) • Transients whose duration is longer than the length of the window : Exponential, Hann • General-purpose applications : Hann • Spectral analysis (frequency-response measurements) : Hann (for random excitation), Rectangular (for pseudorandom excitation) • Separation of two tones with frequencies very close to each other but with widely differing amplitudes : Kaiser-Bessel • Separation of two tones with frequencies very close to each other but with almost equal amplitudes : Rectangular • Accurate single-tone amplitude measurements : Flat top • Sine wave or combination of sine waves : Hann • Sine wave and amplitude accuracy is important : Flat top • Narrowband random signal (vibration data) : Hann • Broadband random (white noise) : Uniform • Closely spaced sine waves : Uniform, Hamming • Excitation signals (hammer blow) : Force • Response signals : Exponential • Unknown content : Hann Once the windowing is executed, you can run the FFT. The result of this algorithm lies in to vectors containing the real and imaginary computed values. You need then to apply some math to convert them into a vector of intensities. Window type Rectangle (box car)

Window type Hamming:

Window type Flat top

All plots exported from my DSP tool box “Panorama”

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

01 02 03 04 05 06 07 08 09 10 11

// 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 1 const uint8_t samples = 64; 2 double vReal[samples]; 3 double vImag[samples]; And this is the weighing routine

01 02 03 04 05 06 07 08 09 10 11 12 13

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

14

weighingFactor = 0.54 - (0.46 * cos(2.0 * pi * ratio));

15 16 17

break; case WIN_TYP_HANN: // hann weighingFactor = 0.54 * (1.0 - cos(2.0 * pi * ratio));

18 19

break; case WIN_TYP_TRIANGLE: // triangle (Bartlett) weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne 20 (samplesMinusOne / 2.0))) / samplesMinusOne); 21 break; 22 case WIN_TYP_BLACKMAN: // blackmann weighingFactor = 0.42323 - (0.49755 * (cos(2.0 * pi 23 * ratio))) + (0.07922 * (cos(4.0 * pi * ratio))); 24 break; 25 case WIN_TYP_FLT_TOP: // flat top weighingFactor = 0.2810639 - (0.5208972 * cos(2.0 * 26 pi * ratio)) + (0.1980399 * cos(4.0 * pi * ratio)); 27 break; 28 case WIN_TYP_WELCH: // welch weighingFactor = 1.0 - sq((indexMinusOne 29 samplesMinusOne / 2.0) / (samplesMinusOne / 2.0)); 30 break; 31 } 32 if (dir == FORWARD) { 33 vData[i] *= weighingFactor; 34 vData[samples - (i + 1)] *= weighingFactor; 35 } 36 else { 37 vData[i] /= weighingFactor; 38 vData[samples - (i + 1)] /= weighingFactor; 39 } 40 } 41 } 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: , FFT is reversible, so that we can apply it in FORWARD or REVERSE mode FFT stands for fast Fourier Transform. The DFT (Direct Fourier Transform) applies to vectors containing any number of signal samples. But it is sssssssssllllllllllllllllllllllloooooooooooooowwwwwwww due to the repeated number of operation. Computing a DFT of n points takes n^2 arithmetical operations, while an FFT can compute the same result in only n log2(n) operations. The implemented algorithm is the “Cooley–Tukey”

algorithm which the far most popular. The “only” limitation is that the number of samples in the signal must be a power of two. However, if the number of samples is less than that, it is possible to replace missing data with 0 values without affecting the final result: this is the zero padding. 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

void PlainFFT::compute(double *vR, double *vI, uint8_t samples, uint8_t dir) { // This computes an in-place complex-to-complex FFT // dir = 1 gives forward transform, dir = 0 gives reverse transform // Reverse bits uint8_t j = 0; for (uint8_t i = 0; i < (samples - 1); i++) { if (i < j) { swap(&vR[i], &vR[j]); swap(&vI[i], &vI[j]); } uint8_t k = (samples >> 1); while (k <= j) { j -= k; k >>= 1; } j += k; } // Compute the FFT double c1 = -1.0; double c2 = 0.0; uint8_t l2 = 1; for (uint8_t l = 0; l < exponent(samples); l++) { uint8_t l1 = l2; l2 <<= 1; double u1 = 1.0; double u2 = 0.0; for (j = 0; j < l1; j++) { for (uint8_t i = j; i < samples; i += l2) { uint8_t i1 = i + l1; double t1 = u1 * vR[i1] - u2 * vI[i1]; double t2 = u1 * vI[i1] + u2 * vR[i1]; vR[i1] = vR[i] - t1; vI[i1] = vI[i] - t2; vR[i] += t1; vI[i] += t2; } double z = (u1 * c1) - (u2 * c2); u2 = (u1 * c2) + (u2 * c1);

39 40 41 42 43 44 45 46 47 48 49 50 51 52 } Notes:

u1 = z; } c2 = sqrt((1.0 - c1) / 2.0); if (dir == FORWARD) c2 = -c2; c1 = sqrt((1.0 + c1) / 2.0); } // Scaling for forward transform if (dir == FORWARD) { for (uint8_t i = 0; i < samples; i++) { vR[i] /= samples; vI[i] /= samples; } }

• The bit reversal routine could be performed once for all FFT's using vectors of the same size, but to the cost of an extra vector which our limited memory cannot afford. • The dir parameter stands for direction: , FFT is reversible, so that we can apply it in FORWARD or REVERSE mode • Resulting Imaginary and Real data are output in the input data vectors Once the FFT algorithm is executed, we get complex values, made of imaginary values and real values. We need to apply some maths in order to convert them in magnitude values. This is quite simply done thanks to the following routine which converts the complex (so as to say rectangular) values into a polar value: void PlainFFT::complexToMagnitude(double *vR, double *vI, uint8_t samples) { 2 // vM is half the size of vR and vI 3 for (uint8_t i = 0; i < samples; i++) { 4 vR[i] = sqrt(sq(vR[i]) + sq(vI[i])); 5 } 6 } Read this document about conversion (p 6) 1

Computing the phase information is quite easy, but I do not really care about it and it requires the arctangent() trigonometric which does not come standard in Arduino language. However if you need it, you will have to use the math library) and compute the pahse values using the following formula: phase[FFT(A)]=arctangent(|Imag[FFT(A)]| / |Real[FFT(A)]|).

So far so good, we got our magnitude spectrum. We may now analyze it just like any other spectrum in order to get meaningful information. In our case, we may want to identify the main frequencies from the spectrum. Firstly we will run a peak picking algorithm in order to locate major peaks, and then keep the most intense one. Those who want to perform spectral comparison will have to keep all significant peaks.

There are lots of peak picking algorithm! But the one we need has to be small and fast. This one is pretty trivial but it works great. Secondly, because of the poor peak shapes (lack of large vectors of data) we need to interpolate the peak apex in order to get an accurate frequency. I am using a non iterative quadratic interpolation which gives good results too. Both routines are merged in one function 01 double majorPeak(double *vD, uint8_t samples) { 02 double maxY = 0; 03 uint8_t IndexOfMaxY = 0; 04 for (uint8_t i = 1; i < (samples - 1); i++) { 05 if ((vD[i-1] < vD[i]) && (vD[i] > vD[i+1])) { 06 if (vD[i] > maxY) { 07 maxY = vD[i]; 08 IndexOfMaxY = i; 09 } 10 } 11 } double delta = 0.5 * ((vD[IndexOfMaxY-1] 12 vD[IndexOfMaxY+1]) / (vD[IndexOfMaxY-1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY+1])); double interpolatedX = (IndexOfMaxY + delta) / 13 samplingDuration; 14 return(interpolatedX); 15 }

PlainFFT is a simple but effective library which contains all the previously described functions for running FFT on vectors of data. Here is the example code from the library files which demonstrates the (pretty easy) use of the FFT. The content of the vectors of interest is printed on completion of each FFT stage 01 #include "PlainFFT.h" 02 03 PlainFFT FFT = PlainFFT(); // Create FFT object // These values can be changed in order to evaluate the 04 functions 05 const uint16_t samples = 64; 06 double signalFrequency = 1000; 07 double samplingFrequency = 5000; 08 uint8_t signalIntensity = 100; 09 // These are input and output vectors 10 double vReal[samples]; 11 double vImag[samples]; 12 uint8_t runOnce = 0x00;

13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

#define SCL_INDEX 0x00 #define SCL_TIME 0x01 #define SCL_FREQUENCY 0x02 void setup(){ Serial.begin(115200); Serial.println("Ready"); } void loop() { if (runOnce == 0x00) { runOnce = 0x01; // Build raw data double cycles = (((samples-1) * signalFrequency) / samplingFrequency); for (uint8_t i = 0; i < samples; i++) { vReal[i] = uint8_t((signalIntensity * (sin((i * (6.2831 * cycles)) / samples) + 1.0)) / 2.0); } printVector(vReal, samples, SCL_TIME); FFT.windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); // Weigh data printVector(vReal, samples, SCL_TIME); FFT.compute(vReal, vImag, samples, FFT_FORWARD); // Compute FFT printVector(vReal, samples, SCL_INDEX); printVector(vImag, samples, SCL_INDEX); FFT.complexToMagnitude(vReal, vImag, samples); // Compute magnitudes printVector(vReal, (samples >> 1), SCL_FREQUENCY); double x = FFT.majorPeak(vReal, samples, samplingFrequency); Serial.println(x, 6); } }

40 41 42 43 44 void printVector(double *vD, uint8_t n, uint8_t scaleType) { 45 double timeInterval = (1.0 / samplingFrequency); 46 for (uint16_t i = 0; i < n; i++) { 47 // Print abscissa value 48 switch (scaleType) { 49 case SCL_INDEX: 50 Serial.print(i, DEC);

51 52 53 54 55 56

break; case SCL_TIME: Serial.print((i * timeInterval), 6); break; case SCL_FREQUENCY: Serial.print((i / (timeInterval * (samples-1))), 6);

57 break; 58 } 59 Serial.print(" "); 60 // Print ordinate value 61 Serial.print(vD[i], 6); 62 Serial.println(); 63 } 64 Serial.println(); 65 } Attention: • The distributed code is slightly different from the posted code due to last minute changes and optimizations • Samples value must be a power of 2 • Make sure you have enough memory before incresing the size of vectors • Vectors must contain unsigned integers (0 to 127), you can easily change that to 16 bits integers (signed or not) to the cost of doubling the vectors size... • the Nyquist theorem when changing the sampling and the signal frequencies Please check this page if you are interested in the code. Feel free to your comments and suggestions. Those who exercized the library and tried to change the default value from the customizable variables (and most probably the size of both data vectors) probably faced some lockups, inifinite loops or unpredicted operations. In most cases, the answer to these problems lies in the limited size of the arduino memory… Let’s be positive and study this subject. Depending upon the type of Arduino platform, you may use one of these microprocessors: • Arduino Diecimila: ATMega168 (data sheet) • Arduino Duemilanove: ATMega368 (data sheet) • Arduino Uno: ATmega328 (data sheet) For other platforms, follow this thread For older platforms, follow this thread Each of these microprocessors has three separate memories: • FLASH memory: this is a non volatile memory space where the sketch is stored • SRAM (Static Random Access Memory): this is the memory sapce where the sketch creates and manipulates variables at run time • EEPROM (Electrically-Erasable Programmable Read-Only Memory): this is a non volatile memory space where data can be read/write at run time. This memory has limitations: the

number of read/write cycle is estimated to be no more than 100.000, and it takes real time to get access to it Let’s go back to our problem. The size of the SRAM (where the vectors are created at run time) depends upon the microprocessor: • ATMega168: 1024 bytes • ATMega368: 2048 bytes That is not much… and every byte counts! Now let’s consider the vectors. Firstly, we have to decide about the type of data to be recorded: • • • •

boolean or char or unsigned char or byte or uint8_t: 1 byte each int or unsigned int or uint16_t or int16_t: 2 bytes each long or unsigned long or uint32_t or int32_t: 4 bytes each float or double: 4 bytes each

For more information on data types, follow this thread This means that the following vectors will occupy: 1 byte vX[32]; // 32*8=256 bytes 2 int vX[32]; // 32*16=512 bytes 3 float vX[32]; // 32*32=1024 bytes So that it is a good idea to check the available memory space prior to running sketches. This how you can do this: 01 int memoryTest() { // This function will return the number of bytes currently free 02 in SRAM 03 int byteCounter = 0; // initialize a counter 04 byte *byteArray; // create a pointer to a byte array // Use the malloc function to repeatedly attempt allocating a 05 certain number of bytes to memory while ( (byteArray = (byte*) malloc (byteCounter * 06 sizeof(byte))) != NULL ) { byteCounter++; // If allocation was successful, then up the 07 count for the next try 08 free(byteArray); // Free memory after allocating it 09 } free(byteArray); // Also free memory after the function 10 finishes return byteCounter; // Send back the number number of bytes 11 which have been successfully allocated 12 } This code has been written by Rob Faludi

I received some questions related to the use of the FFT library. This example illustrate how to interface the FFT function to an acquisition engine, such as the optimized one from the PlainADC library.

01 /* 02 03 04 05

Example of use of the ADC and FFT libraries Copyright (C) 2010 Didier Longueville

This program is free software: you can redistribute it and/or modify it under the of the GNU General Public License as 07 published by the Free Software Foundation, either version 3 of the 08 License, or 09 (at your option) any later version. 10 This program is distributed in the hope that it will be 11 useful, but WITHOUT ANY WARRANTY; without even the implied warranty 12 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See 13 the 14 GNU General Public License for more details. 15 You should have received a copy of the GNU General Public 16 License along with this program. If not, see . 18 19 */ 20 21 22 23 /* Printing function options */ 24 #define SCL_INDEX 0x00 25 #define SCL_TIME 0x01 26 #define SCL_FREQUENCY 0x02 27 28 #include 29 #include 30 PlainADC ADC = PlainADC(); /* Create ADC object */ 31 PlainFFT FFT = PlainFFT(); /* Create FFT object */ 32 33 /* defined variables */ 34 const uint16_t samples = 128; 35 uint16_t frequency = 20000; 36 uint8_t channel = 0; 06

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

/* Data vectors */ uint8_t vData[samples]; double vReal[samples]; double vImag[samples]; void setup(){ /* Initialize serial comm port */ Serial.begin(115200); // /* Set acquisition parameters */ ADC.setAcquisitionParameters(channel, samples, frequency); } void loop() { /* Acquire data and store them in a vector of bytes */ ADC.acquireData(vData); /* Convert 8 bits unsigned data in 32 bits floats */ for (uint16_t i = 0; i < samples; i++) { vReal[i] = double(vData[i]); } /* Weigh data */ FFT.windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Compute FFT: vReal and vImag vectors contain the source data and will contain the result data on completion of executing the function*/ FFT.compute(vReal, vImag, samples, FFT_FORWARD); /* Compute magnitudes: the resulting data can be read from the vReal vector */ FFT.complexToMagnitude(vReal, vImag, samples); /* frequency spectrum */ printVector(vReal, (samples >> 1), SCL_FREQUENCY); /* Pause */ delay(5000); }

62 63 64 65 66 67 68 69 void printVector(double *vD, uint8_t n, uint8_t scaleType) { 70 /* Mulitpurpose printing function */ 71 double timeInterval = (1.0 / frequency); 72 for (uint16_t i = 0; i < n; i++) { 73 /* Print abscissa value */ 74 switch (scaleType) { 75 case SCL_INDEX: Serial.print(i, DEC); break;

76 77 78 79 80 81 82 83 84 85 86

case SCL_TIME: Serial.print((i * timeInterval), 6); break; case SCL_FREQUENCY: Serial.print((i / (timeInterval * (samples-1))), 6); break; } Serial.print(" "); /* Print ordinate value */ Serial.print(vD[i], 6); Serial.println(); } Serial.println(); }

Here are a few comments on the use of PlainFFT. First of all: PlainFFT performs in place calculation, which means that the results are recorded within the vectors containing the source data. “Usable” results are available on completion of the ComplexToReal() execution. Before execution, vReal and vImag vectors contain the transformed results from the FFT. After execution, the vReal vector contains… two times the frequency spectrum in a perfectly symetric way. The PlainFFT_01.ino skecth example shows 1 /* Compute magnitudes */ 2 FFT.ComplexToReal(vReal, vImag, samples, FFT_SCL_TYP_AMPLITUDE); While in fact we could save processing time with 1 /* Compute magnitudes */ FFT.ComplexToReal(vReal, vImag, samples >> 1, 2 FFT_SCL_TYP_AMPLITUDE); As long as, in any way the sketch uses 1 /* Print frequency spectrum */ 2 PrintVector(vReal, (samples >> 1), SCL_FREQUENCY); for printing the spectrum. CQFD. The sampling frequency (Let’s call this variable Fs) determines the frequency range of the spectrum while the number of points acquired during signal acquisition (Let’s call this variable N) determines the resolution frequency. As a consequence of that:: - To Increase the frequency range, increase the signal sampling frequency - To increase the frequency resolution for a given frequency range, increase the number of points acquired at the same sampling frequency. The result from the FFT can be viewed as a collection of N/2 “bins”. To each bin corresponds a frequency range (or bandwidth), centred on a frequency value, as reported in the vReal vector. The width of each bin (Let’s call this variable Df) is equal to Fs/N. The last bin is centred on the max frequency (Let’s call this variable Fmax) at (Fs/2)-(Fs/N). The first bin is a little bit special: It is centred on 0Hz, so as to say on the DC component of the signal. So that the graphical bin representation should be half the size of the other bins! All these parameters can be pictured in the following way

Bin 1 Bin 2 Bin 3 Bin N/2 . . . . // . . . . . . // . . . . . . // . . . : . : . : . // . : . | | | | | | | | | -| | | ------| | --------------------| -----------------------------------------------------------

Nyquist frequency(Fs/2) Max. Frequency (Fmax) Df * 2 Df * 1 Df * 0, 0Hz (DC)

In the real world, the amount of signal resulting from: - Df * 0 +/- Df /2, in other words 0hz, is read from vReal [0] - Df * 1 +/- Df /2 is read from vReal [1] - Df * 2 +/- Df /2 is read from vReal [2] And so on, up to - Df * N/2 +/- Df /2, so as to say Fmax, is read from vReal [N/2] Each bin contains the power of one or more signals which frequencies fit in the band bandwith. From this point, it is obvious that the narrower the bin, the better the selectivity of the filter. In some, not to say most, circumstances, the spectrum representation does not exactly looks like pure histograms. For a given signal characterized by a pure sinusoidal function, the most abundant power is read from the bin which corresponds to the signal frequency, but the adjacent bins may contain some power abundances too. This is due to the frequency leakage. We may use one from the next options to overcome this phenomenon: - Increase the number of samples (N) in order to improve the selectivity - Apply appropriate windowing - Apply spectrum peak interpolation Peak interpolation helps improving the most accurate peak apex location thanks to some linear (e.g. Quadratic fit) or non-linear (e.g. Gaussian fit). While this interpolation applies to a limited number of consecutive data points, we may apply simplified formula for computing apex locations. Peak interpolation is almost mandatory when all significant peaks from a frequency spectrum must be correlated (e.g. Finding harmonics). The following pictures illustrate various use cases of PlainFFT for plotting frequency spectrum. In all cases, a Hann window was applied as well as Gaussian peak apex interpolation. Here is the frequency spectrum from a sinusoidal signal acquired with a sampling rate of 1kHz over 256 samples:

Here is the frequency spectrum from a sinusoidal signal acquired with a sampling rate of 1kHz over 32768 samples:

Here is the frequency spectrum from a sinusoidal signal acquired with a sampling rate of 500Hz over 256 samples:

Note: The peak looks broader than in picture 1, but this feeling is related to the narrower scale! This is probably on over-simplified explanation, but it is sufficient for using PlainFFT in the most appropriate manner. Feel free to add your comments and suggestions.