mirror
/
arduinoFFT
kopia lustrzana https://github.com/kosme/arduinoFFT
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

Remove deprecated functions, templatize, Speedup

pull/42/head
Bim Overbohm 2020-02-19 17:15:49 +01:00
rodzic 8459c48952
commit cb33149c17
13 zmienionych plików z 673 dodań i 634 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

19
Examples/FFT_01/FFT_01.ino
Wyświetl plik

@ -1,7 +1,9 @@
/*
Example of use of the FFT libray
Copyright (C) 2014 Enrique Condes
Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@ -30,7 +32,6 @@
#include "arduinoFFT.h"
arduinoFFT FFT = arduinoFFT(); /* Create FFT object */
/*
These values can be changed in order to evaluate the functions
*/
@ -38,6 +39,7 @@ const uint16_t samples = 64; //This value MUST ALWAYS be a power of 2
const double signalFrequency = 1000;
const double samplingFrequency = 5000;
const uint8_t amplitude = 100;
/*
These are the input and output vectors
Input vectors receive computed results from FFT
@ -45,6 +47,9 @@ Input vectors receive computed results from FFT
double vReal[samples];
double vImag[samples];
/* Create FFT object */
ArduinoFFT<double> FFT = ArduinoFFT<double>(vReal, vImag, samples, samplingFrequency);
#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02
@ -62,25 +67,25 @@ void loop()
double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read
for (uint16_t i = 0; i < samples; i++)
{
vReal[i] = int8_t((amplitude * (sin((i * (twoPi * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/
vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/
//vReal[i] = uint8_t((amplitude * (sin((i * (twoPi * cycles)) / samples) + 1.0)) / 2.0);/* Build data displaced on the Y axis to include only positive values*/
vImag[i] = 0.0; //Imaginary part must be zeroed in case of looping to avoid wrong calculations and overflows
}
/* Print the results of the simulated sampling according to time */
Serial.println("Data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.Windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */
FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */
Serial.println("Weighed data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.Compute(vReal, vImag, samples, FFT_FORWARD); /* Compute FFT */
FFT.compute(FFTDirection::Forward); /* Compute FFT */
Serial.println("Computed Real values:");
PrintVector(vReal, samples, SCL_INDEX);
Serial.println("Computed Imaginary values:");
PrintVector(vImag, samples, SCL_INDEX);
FFT.ComplexToMagnitude(vReal, vImag, samples); /* Compute magnitudes */
FFT.complexToMagnitude(); /* Compute magnitudes */
Serial.println("Computed magnitudes:");
PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);
double x = FFT.MajorPeak(vReal, samples, samplingFrequency);
double x = FFT.majorPeak();
Serial.println(x, 6);
while(1); /* Run Once */
// delay(2000); /* Repeat after delay */

21
Examples/FFT_02/FFT_02.ino
Wyświetl plik

@ -4,7 +4,9 @@
The exponent is calculated once before the excecution since it is a constant.
This saves resources during the excecution of the sketch and reduces the compiled size.
The sketch shows the time that the computing is taking.
Copyright (C) 2014 Enrique Condes
Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@ -23,15 +25,12 @@
#include "arduinoFFT.h"
arduinoFFT FFT = arduinoFFT(); /* Create FFT object */
/*
These values can be changed in order to evaluate the functions
*/
const uint16_t samples = 64;
const double sampling = 40;
const uint8_t amplitude = 4;
uint8_t exponent;
const double startFrequency = 2;
const double stopFrequency = 16.4;
const double step_size = 0.1;
@ -43,6 +42,9 @@ Input vectors receive computed results from FFT
double vReal[samples];
double vImag[samples];
/* Create FFT object */
ArduinoFFT<double> FFT = ArduinoFFT<double>(vReal, vImag, samples, sampling);
unsigned long time;
#define SCL_INDEX 0x00
@ -54,7 +56,6 @@ void setup()
{
Serial.begin(115200);
Serial.println("Ready");
exponent = FFT.Exponent(samples);
}
void loop()
@ -67,24 +68,24 @@ void loop()
double cycles = (((samples-1) * frequency) / sampling);
for (uint16_t i = 0; i < samples; i++)
{
vReal[i] = int8_t((amplitude * (sin((i * (twoPi * cycles)) / samples))) / 2.0);
vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0);
vImag[i] = 0; //Reset the imaginary values vector for each new frequency
}
/*Serial.println("Data:");
PrintVector(vReal, samples, SCL_TIME);*/
time=millis();
FFT.Windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */
FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */
/*Serial.println("Weighed data:");
PrintVector(vReal, samples, SCL_TIME);*/
FFT.Compute(vReal, vImag, samples, exponent, FFT_FORWARD); /* Compute FFT */
FFT.compute(FFTDirection::Forward); /* Compute FFT */
/*Serial.println("Computed Real values:");
PrintVector(vReal, samples, SCL_INDEX);
Serial.println("Computed Imaginary values:");
PrintVector(vImag, samples, SCL_INDEX);*/
FFT.ComplexToMagnitude(vReal, vImag, samples); /* Compute magnitudes */
FFT.complexToMagnitude(); /* Compute magnitudes */
/*Serial.println("Computed magnitudes:");
PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);*/
double x = FFT.MajorPeak(vReal, samples, sampling);
double x = FFT.majorPeak();
Serial.print(frequency);
Serial.print(": \t\t");
Serial.print(x, 4);

17
Examples/FFT_03/FFT_03.ino
Wyświetl plik

@ -1,7 +1,9 @@
/*
Example of use of the FFT libray to compute FFT for a signal sampled through the ADC.
Copyright (C) 2018 Enrique Condés and Ragnar Ranøyen Homb
Copyright (C) 2018 Enrique Condés and Ragnar Ranøyen Homb
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@ -20,14 +22,12 @@
#include "arduinoFFT.h"
arduinoFFT FFT = arduinoFFT(); /* Create FFT object */
/*
These values can be changed in order to evaluate the functions
*/
#define CHANNEL A0
const uint16_t samples = 64; //This value MUST ALWAYS be a power of 2
const double samplingFrequency = 100; //Hz, must be less than 10000 due to ADC
unsigned int sampling_period_us;
unsigned long microseconds;
@ -38,6 +38,9 @@ Input vectors receive computed results from FFT
double vReal[samples];
double vImag[samples];
/* Create FFT object */
ArduinoFFT<double> FFT = ArduinoFFT<double>(vReal, vImag, samples, samplingFrequency);
#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02
@ -66,18 +69,18 @@ void loop()
/* Print the results of the sampling according to time */
Serial.println("Data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.Windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */
FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */
Serial.println("Weighed data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.Compute(vReal, vImag, samples, FFT_FORWARD); /* Compute FFT */
FFT.compute(FFTDirection::Forward); /* Compute FFT */
Serial.println("Computed Real values:");
PrintVector(vReal, samples, SCL_INDEX);
Serial.println("Computed Imaginary values:");
PrintVector(vImag, samples, SCL_INDEX);
FFT.ComplexToMagnitude(vReal, vImag, samples); /* Compute magnitudes */
FFT.complexToMagnitude(); /* Compute magnitudes */
Serial.println("Computed magnitudes:");
PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);
double x = FFT.MajorPeak(vReal, samples, samplingFrequency);
double x = FFT.majorPeak();
Serial.println(x, 6); //Print out what frequency is the most dominant.
while(1); /* Run Once */
// delay(2000); /* Repeat after delay */

18
Examples/FFT_04/FFT_04.ino
Wyświetl plik

@ -1,7 +1,9 @@
/*
Example of use of the FFT libray
Copyright (C) 2018 Enrique Condes
Copyright (C) 2018 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@ -31,7 +33,6 @@
#include "arduinoFFT.h"
arduinoFFT FFT = arduinoFFT(); /* Create FFT object */
/*
These values can be changed in order to evaluate the functions
*/
@ -39,6 +40,7 @@ const uint16_t samples = 64; //This value MUST ALWAYS be a power of 2
const double signalFrequency = 1000;
const double samplingFrequency = 5000;
const uint8_t amplitude = 100;
/*
These are the input and output vectors
Input vectors receive computed results from FFT
@ -46,6 +48,8 @@ Input vectors receive computed results from FFT
double vReal[samples];
double vImag[samples];
ArduinoFFT<double> FFT = ArduinoFFT<double>(vReal, vImag, samples, samplingFrequency);
#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02
@ -62,15 +66,15 @@ void loop()
double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read
for (uint16_t i = 0; i < samples; i++)
{
vReal[i] = int8_t((amplitude * (sin((i * (twoPi * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/
vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/
//vReal[i] = uint8_t((amplitude * (sin((i * (twoPi * cycles)) / samples) + 1.0)) / 2.0);/* Build data displaced on the Y axis to include only positive values*/
vImag[i] = 0.0; //Imaginary part must be zeroed in case of looping to avoid wrong calculations and overflows
}
FFT.Windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */
FFT.Compute(vReal, vImag, samples, FFT_FORWARD); /* Compute FFT */
FFT.ComplexToMagnitude(vReal, vImag, samples); /* Compute magnitudes */
FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */
FFT.compute(FFTDirection::Forward); /* Compute FFT */
FFT.complexToMagnitude(); /* Compute magnitudes */
PrintVector(vReal, samples>>1, SCL_PLOT);
double x = FFT.MajorPeak(vReal, samples, samplingFrequency);
double x = FFT.majorPeak();
while(1); /* Run Once */
// delay(2000); /* Repeat after delay */
}

19
Examples/FFT_05/FFT_05.ino
Wyświetl plik

@ -1,7 +1,9 @@
/*
Example of use of the FFT libray
Copyright (C) 2014 Enrique Condes
Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@ -31,7 +33,6 @@
#include "arduinoFFT.h"
arduinoFFT FFT = arduinoFFT(); /* Create FFT object */
/*
These values can be changed in order to evaluate the functions
*/
@ -39,6 +40,7 @@ const uint16_t samples = 64; //This value MUST ALWAYS be a power of 2
const double signalFrequency = 1000;
const double samplingFrequency = 5000;
const uint8_t amplitude = 100;
/*
These are the input and output vectors
Input vectors receive computed results from FFT
@ -46,6 +48,9 @@ Input vectors receive computed results from FFT
double vReal[samples];
double vImag[samples];
/* Create FFT object */
ArduinoFFT<double> FFT = ArduinoFFT<double>(vReal, vImag, samples, samplingFrequency);
#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02
@ -63,27 +68,27 @@ void loop()
double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read
for (uint16_t i = 0; i < samples; i++)
{
vReal[i] = int8_t((amplitude * (sin((i * (twoPi * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/
vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/
//vReal[i] = uint8_t((amplitude * (sin((i * (twoPi * cycles)) / samples) + 1.0)) / 2.0);/* Build data displaced on the Y axis to include only positive values*/
vImag[i] = 0.0; //Imaginary part must be zeroed in case of looping to avoid wrong calculations and overflows
}
/* Print the results of the simulated sampling according to time */
Serial.println("Data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.Windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */
FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */
Serial.println("Weighed data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.Compute(vReal, vImag, samples, FFT_FORWARD); /* Compute FFT */
FFT.compute(FFTDirection::Forward); /* Compute FFT */
Serial.println("Computed Real values:");
PrintVector(vReal, samples, SCL_INDEX);
Serial.println("Computed Imaginary values:");
PrintVector(vImag, samples, SCL_INDEX);
FFT.ComplexToMagnitude(vReal, vImag, samples); /* Compute magnitudes */
FFT.complexToMagnitude(); /* Compute magnitudes */
Serial.println("Computed magnitudes:");
PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);
double x;
double v;
FFT.MajorPeak(vReal, samples, samplingFrequency, &x, &v);
FFT.majorPeak(x, v);
Serial.print(x, 6);
Serial.print(", ");
Serial.println(v, 6);

129
Examples/FFT_speedup/FFT_speedup.ino 100644
Wyświetl plik

@ -0,0 +1,129 @@
/*
Example of use of the FFT libray to compute FFT for a signal sampled through the ADC
with speedup through different arduinoFFT options. Based on examples/FFT_03/FFT_03.ino
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
// There are two speedup options for some of the FFT code:
// Define this to use reciprocal multiplication for division and some more speedups that might decrease precision
//#define FFT_SPEED_OVER_PRECISION
// Define this to use a low-precision square root approximation instead of the regular sqrt() call
// This might only work for specific use cases, but is significantly faster. Only works for ArduinoFFT<float>.
//#define FFT_SQRT_APPROXIMATION
#include "arduinoFFT.h"
/*
These values can be changed in order to evaluate the functions
*/
#define CHANNEL A0
const uint16_t samples = 64; //This value MUST ALWAYS be a power of 2
const float samplingFrequency = 100; //Hz, must be less than 10000 due to ADC
unsigned int sampling_period_us;
unsigned long microseconds;
/*
These are the input and output vectors
Input vectors receive computed results from FFT
*/
float vReal[samples];
float vImag[samples];
/*
Allocate space for FFT window weighing factors, so they are calculated only the first time windowing() is called.
If you don't do this, a lot of calculations are necessary, depending on the window function.
*/
float weighingFactors[samples];
/* Create FFT object with weighing factor storage */
ArduinoFFT<float> FFT = ArduinoFFT<float>(vReal, vImag, samples, samplingFrequency, weighingFactors);
#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02
#define SCL_PLOT 0x03
void setup()
{
sampling_period_us = round(1000000*(1.0/samplingFrequency));
Serial.begin(115200);
Serial.println("Ready");
}
void loop()
{
/*SAMPLING*/
microseconds = micros();
for(int i=0; i<samples; i++)
{
vReal[i] = analogRead(CHANNEL);
vImag[i] = 0;
while(micros() - microseconds < sampling_period_us){
//empty loop
}
microseconds += sampling_period_us;
}
/* Print the results of the sampling according to time */
Serial.println("Data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */
Serial.println("Weighed data:");
PrintVector(vReal, samples, SCL_TIME);
FFT.compute(FFTDirection::Forward); /* Compute FFT */
Serial.println("Computed Real values:");
PrintVector(vReal, samples, SCL_INDEX);
Serial.println("Computed Imaginary values:");
PrintVector(vImag, samples, SCL_INDEX);
FFT.complexToMagnitude(); /* Compute magnitudes */
Serial.println("Computed magnitudes:");
PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);
float x = FFT.majorPeak();
Serial.println(x, 6); //Print out what frequency is the most dominant.
while(1); /* Run Once */
// delay(2000); /* Repeat after delay */
}
void PrintVector(float *vData, uint16_t bufferSize, uint8_t scaleType)
{
for (uint16_t i = 0; i < bufferSize; i++)
{
float abscissa;
/* Print abscissa value */
switch (scaleType)
{
case SCL_INDEX:
abscissa = (i * 1.0);
break;
case SCL_TIME:
abscissa = ((i * 1.0) / samplingFrequency);
break;
case SCL_FREQUENCY:
abscissa = ((i * 1.0 * samplingFrequency) / samples);
break;
}
Serial.print(abscissa, 6);
if(scaleType==SCL_FREQUENCY)
Serial.print("Hz");
Serial.print(" ");
Serial.println(vData[i], 4);
}
Serial.println();
}

98
README.md
Wyświetl plik

@ -1,32 +1,28 @@
arduinoFFT
==========
Fast Fourier Transform for Arduino
# Fast Fourier Transform for Arduino
This is a fork from https://code.google.com/p/makefurt/ which has been abandoned since 2011.
~~This is a C++ library for Arduino for computing FFT.~~ Now it works both on Arduino and C projects.
Tested on Arduino 1.6.11 and 1.8.10.
<del>This is a C++ library for Arduino for computing FFT.</del> Now it works both on Arduino and C projects.
Tested on Arduino 1.6.11
### Installation on Arduino
## Installation on Arduino
Use the Arduino Library Manager to install and keep it updated. Just look for arduinoFFT. Only for Arduino 1.5+
### Manual installation on Arduino
## Manual installation on Arduino
To install this library, just place this entire folder as a subfolder in your Arduino installation
To install this library, just place this entire folder as a subfolder in your Arduino installation. When installed, this library should look like:
When installed, this library should look like:
`Arduino\libraries\arduinoFTT` (this library's folder)
`Arduino\libraries\arduinoFTT\arduinoFTT.h` (the library header file, uses 32 bit floats or 64bit doubles)
`Arduino\libraries\arduinoFTT\keywords.txt` (the syntax coloring file)
`Arduino\libraries\arduinoFTT\examples` (the examples in the "open" menu)
`Arduino\libraries\arduinoFTT\LICENSE` (GPL license file)
`Arduino\libraries\arduinoFTT\README.md` (this file)
Arduino\libraries\arduinoFTT (this library's folder)
Arduino\libraries\arduinoFTT\arduinoFTT.cpp (the library implementation file, uses 32 bits floats vectors)
Arduino\libraries\arduinoFTT\arduinoFTT.h (the library header file, uses 32 bits floats vectors)
Arduino\libraries\arduinoFTT\keywords.txt (the syntax coloring file)
Arduino\libraries\arduinoFTT\examples (the examples in the "open" menu)
Arduino\libraries\arduinoFTT\readme.md (this file)
### Building on Arduino
## Building on Arduino
After this library is installed, you just have to start the Arduino application.
You may see a few warning messages as it's built.
@ -36,46 +32,50 @@ select arduinoFTT. This will add a corresponding line to the top of your sketch
`#include <arduinoFTT.h>`
### TODO
## TODO
* Ratio table for windowing function.
* Document windowing functions advantages and disadvantages.
* Optimize usage and arguments.
* Add new windowing functions.
<del>* Spectrum table? </del>
* ~~Spectrum table?~~
### API
## API
* **arduinoFFT**(void);
* **arduinoFFT**(double *vReal, double *vImag, uint16_t samples, double samplingFrequency);
Constructor
* **~arduinoFFT**(void);
Destructor
* **ComplexToMagnitude**(double *vReal, double *vImag, uint16_t samples);
* **ComplexToMagnitude**();
* **Compute**(double *vReal, double *vImag, uint16_t samples, uint8_t dir);
* **Compute**(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir);
* **Compute**(uint8_t dir);
* **ArduinoFFT**(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, T * weighingFactors = nullptr);
Constructor.
The type `T` can be `float` or `double`. `vReal` and `vImag` are pointers to arrays of real and imaginary data and have to be allocated outside of ArduinoFFT. `samples` is the number of samples in `vReal` and `vImag` and `weighingFactors` (if specified). `samplingFrequency` is the sample frequency of the data. `weighingFactors` can optionally be specified to cache weighing factors for the windowing function. This speeds up repeated calls to **windowing()** significantly.
* **~ArduinoFFT**(void);
Destructor.
* **complexToMagnitude**();
Convert complex values to their magnitude and store in vReal.
* **compute**(FFTDirection dir);
Calcuates the Fast Fourier Transform.
* **DCRemoval**(double *vData, uint16_t samples);
* **DCRemoval**();
* **dcRemoval**();
Removes the DC component from the sample data.
* **MajorPeak**(double *vD, uint16_t samples, double samplingFrequency);
* **MajorPeak**();
* **majorPeak**();
Looks for and returns the frequency of the biggest spike in the analyzed signal.
* **Revision**(void);
* **revision**();
Returns the library revision.
* **Windowing**(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir);
* **Windowing**(uint8_t windowType, uint8_t dir);
* **windowing**(FFTWindow windowType, FFTDirection dir);
Performs a windowing function on the values array. The possible windowing options are:
* FFT_WIN_TYP_RECTANGLE
* FFT_WIN_TYP_HAMMING
* FFT_WIN_TYP_HANN
* FFT_WIN_TYP_TRIANGLE
* FFT_WIN_TYP_NUTTALL
* FFT_WIN_TYP_BLACKMAN
* FFT_WIN_TYP_BLACKMAN_NUTTALL
* FFT_WIN_TYP_BLACKMAN_HARRIS
* FFT_WIN_TYP_FLT_TOP
* FFT_WIN_TYP_WELCH
* **Exponent**(uint16_t value);
Calculates and returns the base 2 logarithm of the given value.
* Rectangle
* Hamming
* Hann
* Triangle
* Nuttall
* Blackman
* Blackman_Nuttall
* Blackman_Harris
* Flat_top
* Welch
## Special flags
You can define these before including arduinoFFT.h:
* #define FFT_SPEED_OVER_PRECISION
Define this to use reciprocal multiplication for division and some more speedups that might decrease precision.
* #define FFT_SQRT_APPROXIMATION
Define this to use a low-precision square root approximation instead of the regular sqrt() call. This might only work for specific use cases, but is significantly faster. Only works if `T == float`.

6
changeLog.txt
Wyświetl plik

@ -1,3 +1,9 @@
02/19/20 v1.9.0
Remove deprecated API. Consistent renaming of functions to lowercase.
Make template to be able to use float or double type (float brings a ~70% speed increase on ESP32).
Add option to provide cache for window function weighing factors (~50% speed increase on ESP32).
Add some #defines to enable math approximisations to further speed up code (~50% speed increase on ESP32).
01/27/20 v1.5.5
Lookup table for constants c1 and c2 used during FFT comupting. This increases the FFT computing speed in around 5%.

44
keywords.txt
Wyświetl plik

@ -6,35 +6,35 @@
# Datatypes (KEYWORD1)
#######################################
arduinoFFT KEYWORD1
ArduinoFFT KEYWORD1
FFTDirection KEYWORD1
FFTWindow KEYWORD1
#######################################
# Methods and Functions (KEYWORD2)
#######################################
ComplexToMagnitude KEYWORD2
Compute KEYWORD2
DCRemoval KEYWORD2
Windowing KEYWORD2
Exponent KEYWORD2
Revision KEYWORD2
MajorPeak KEYWORD2
complexToMagnitude KEYWORD2
compute KEYWORD2
dcRemoval KEYWORD2
windowing KEYWORD2
exponent KEYWORD2
revision KEYWORD2
majorPeak KEYWORD2
#######################################
# Constants (LITERAL1)
#######################################
twoPi LITERAL1
fourPi LITERAL1
FFT_FORWARD LITERAL1
FFT_REVERSE LITERAL1
FFT_WIN_TYP_RECTANGLE LITERAL1
FFT_WIN_TYP_HAMMING LITERAL1
FFT_WIN_TYP_HANN LITERAL1
FFT_WIN_TYP_TRIANGLE LITERAL1
FFT_WIN_TYP_NUTTALL LITERAL1
FFT_WIN_TYP_BLACKMAN LITERAL1
FFT_WIN_TYP_BLACKMAN_NUTTALL LITERAL1
FFT_WIN_TYP_BLACKMAN_HARRIS LITERAL1
FFT_WIN_TYP_FLT_TOP LITERAL1
FFT_WIN_TYP_WELCH LITERAL1
Forward LITERAL1
Reverse LITERAL1
Rectangle LITERAL1
Hamming LITERAL1
Hann LITERAL1
Triangle LITERAL1
Nuttall LITERAL1
Blackman LITERAL1
Blackman_Nuttall LITERAL1
Blackman_Harris LITERAL1
Flat_top LITERAL1
Welch LITERAL1

7
library.json
Wyświetl plik

@ -18,9 +18,14 @@
"name": "Didier Longueville",
"url": "http://www.arduinoos.com/",
"email": "contact@arduinoos.com"
},
{
"name": "Bim Overbohm",
"url": "https://github.com/HorstBaerbel",
"email": "bim.overbohm@googlemail.com"
}
],
"version": "1.5.5",
"version": "1.9.0",
"frameworks": ["arduino","mbed","espidf"],
"platforms": "*"
}

2
library.properties
Wyświetl plik

@ -1,5 +1,5 @@
name=arduinoFFT
version=1.5.5
version=1.9.0
author=Enrique Condes <enrique@shapeoko.com>
maintainer=Enrique Condes <enrique@shapeoko.com>
sentence=A library for implementing floating point Fast Fourier Transform calculations on Arduino.

446
src/arduinoFFT.cpp
Wyświetl plik

@ -1,446 +0,0 @@
/*
FFT libray
Copyright (C) 2010 Didier Longueville
Copyright (C) 2014 Enrique Condes
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "arduinoFFT.h"
arduinoFFT::arduinoFFT(void)
{ // Constructor
#warning("This method is deprecated and may be removed on future revisions.")
}
arduinoFFT::arduinoFFT(double *vReal, double *vImag, uint16_t samples, double samplingFrequency)
{// Constructor
this->_vReal = vReal;
this->_vImag = vImag;
this->_samples = samples;
this->_samplingFrequency = samplingFrequency;
this->_power = Exponent(samples);
}
arduinoFFT::~arduinoFFT(void)
{
// Destructor
}
uint8_t arduinoFFT::Revision(void)
{
return(FFT_LIB_REV);
}
void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t dir)
{
#warning("This method is deprecated and may be removed on future revisions.")
Compute(vReal, vImag, samples, Exponent(samples), dir);
}
void arduinoFFT::Compute(uint8_t dir)
{// Computes in-place complex-to-complex FFT /
// Reverse bits /
uint16_t j = 0;
for (uint16_t i = 0; i < (this->_samples - 1); i++) {
if (i < j) {
Swap(&this->_vReal[i], &this->_vReal[j]);
if(dir==FFT_REVERSE)
Swap(&this->_vImag[i], &this->_vImag[j]);
}
uint16_t k = (this->_samples >> 1);
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
// Compute the FFT /
#ifdef __AVR__
uint8_t index = 0;
#endif
double c1 = -1.0;
double c2 = 0.0;
uint16_t l2 = 1;
for (uint8_t l = 0; (l < this->_power); l++) {
uint16_t l1 = l2;
l2 <<= 1;
double u1 = 1.0;
double u2 = 0.0;
for (j = 0; j < l1; j++) {
for (uint16_t i = j; i < this->_samples; i += l2) {
uint16_t i1 = i + l1;
double t1 = u1 * this->_vReal[i1] - u2 * this->_vImag[i1];
double t2 = u1 * this->_vImag[i1] + u2 * this->_vReal[i1];
this->_vReal[i1] = this->_vReal[i] - t1;
this->_vImag[i1] = this->_vImag[i] - t2;
this->_vReal[i] += t1;
this->_vImag[i] += t2;
}
double z = ((u1 * c1) - (u2 * c2));
u2 = ((u1 * c2) + (u2 * c1));
u1 = z;
}
#ifdef __AVR__
c2 = pgm_read_float_near(&(_c2[index]));
c1 = pgm_read_float_near(&(_c1[index]));
index++;
#else
c2 = sqrt((1.0 - c1) / 2.0);
c1 = sqrt((1.0 + c1) / 2.0);
#endif
if (dir == FFT_FORWARD) {
c2 = -c2;
}
}
// Scaling for reverse transform /
if (dir != FFT_FORWARD) {
for (uint16_t i = 0; i < this->_samples; i++) {
this->_vReal[i] /= this->_samples;
this->_vImag[i] /= this->_samples;
}
}
}
void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir)
{ // Computes in-place complex-to-complex FFT
// Reverse bits
#warning("This method is deprecated and may be removed on future revisions.")
uint16_t j = 0;
for (uint16_t i = 0; i < (samples - 1); i++) {
if (i < j) {
Swap(&vReal[i], &vReal[j]);
if(dir==FFT_REVERSE)
Swap(&vImag[i], &vImag[j]);
}
uint16_t k = (samples >> 1);
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
// Compute the FFT
#ifdef __AVR__
uint8_t index = 0;
#endif
double c1 = -1.0;
double c2 = 0.0;
uint16_t l2 = 1;
for (uint8_t l = 0; (l < power); l++) {
uint16_t l1 = l2;
l2 <<= 1;
double u1 = 1.0;
double u2 = 0.0;
for (j = 0; j < l1; j++) {
for (uint16_t i = j; i < samples; i += l2) {
uint16_t i1 = i + l1;
double t1 = u1 * vReal[i1] - u2 * vImag[i1];
double t2 = u1 * vImag[i1] + u2 * vReal[i1];
vReal[i1] = vReal[i] - t1;
vImag[i1] = vImag[i] - t2;
vReal[i] += t1;
vImag[i] += t2;
}
double z = ((u1 * c1) - (u2 * c2));
u2 = ((u1 * c2) + (u2 * c1));
u1 = z;
}
#ifdef __AVR__
c2 = pgm_read_float_near(&(_c2[index]));
c1 = pgm_read_float_near(&(_c1[index]));
index++;
#else
c2 = sqrt((1.0 - c1) / 2.0);
c1 = sqrt((1.0 + c1) / 2.0);
#endif
if (dir == FFT_FORWARD) {
c2 = -c2;
}
}
// Scaling for reverse transform
if (dir != FFT_FORWARD) {
for (uint16_t i = 0; i < samples; i++) {
vReal[i] /= samples;
vImag[i] /= samples;
}
}
}
void arduinoFFT::ComplexToMagnitude()
{ // vM is half the size of vReal and vImag
for (uint16_t i = 0; i < this->_samples; i++) {
this->_vReal[i] = sqrt(sq(this->_vReal[i]) + sq(this->_vImag[i]));
}
}
void arduinoFFT::ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples)
{ // vM is half the size of vReal and vImag
#warning("This method is deprecated and may be removed on future revisions.")
for (uint16_t i = 0; i < samples; i++) {
vReal[i] = sqrt(sq(vReal[i]) + sq(vImag[i]));
}
}
void arduinoFFT::DCRemoval()
{
// calculate the mean of vData
double mean = 0;
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
mean += this->_vReal[i];
}
mean /= this->_samples;
// Subtract the mean from vData
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
this->_vReal[i] -= mean;
}
}
void arduinoFFT::DCRemoval(double *vData, uint16_t samples)
{
// calculate the mean of vData
#warning("This method is deprecated and may be removed on future revisions.")
double mean = 0;
for (uint16_t i = 1; i < ((samples >> 1) + 1); i++)
{
mean += vData[i];
}
mean /= samples;
// Subtract the mean from vData
for (uint16_t i = 1; i < ((samples >> 1) + 1); i++)
{
vData[i] -= mean;
}
}
void arduinoFFT::Windowing(uint8_t windowType, uint8_t dir)
{// Weighing factors are computed once before multiple use of FFT
// The weighing function is symetric; half the weighs are recorded
double samplesMinusOne = (double(this->_samples) - 1.0);
for (uint16_t i = 0; i < (this->_samples >> 1); i++) {
double indexMinusOne = double(i);
double ratio = (indexMinusOne / samplesMinusOne);
double weighingFactor = 1.0;
// Compute and record weighting factor
switch (windowType) {
case FFT_WIN_TYP_RECTANGLE: // rectangle (box car)
weighingFactor = 1.0;
break;
case FFT_WIN_TYP_HAMMING: // hamming
weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio));
break;
case FFT_WIN_TYP_HANN: // hann
weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio));
break;
case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne);
break;
case FFT_WIN_TYP_NUTTALL: // nuttall
weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) + (0.144232 * (cos(fourPi * ratio))) - (0.012604 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN: // blackman
weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall
weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) + (0.1365995 * (cos(fourPi * ratio))) - (0.0106411 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris
weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) + (0.14128 * (cos(fourPi * ratio))) - (0.01168 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_FLT_TOP: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) + (0.1980399 * cos(fourPi * ratio));
break;
case FFT_WIN_TYP_WELCH: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0));
break;
}
if (dir == FFT_FORWARD) {
this->_vReal[i] *= weighingFactor;
this->_vReal[this->_samples - (i + 1)] *= weighingFactor;
}
else {
this->_vReal[i] /= weighingFactor;
this->_vReal[this->_samples - (i + 1)] /= weighingFactor;
}
}
}
void arduinoFFT::Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir)
{// Weighing factors are computed once before multiple use of FFT
// The weighing function is symetric; half the weighs are recorded
#warning("This method is deprecated and may be removed on future revisions.")
double samplesMinusOne = (double(samples) - 1.0);
for (uint16_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 FFT_WIN_TYP_RECTANGLE: // rectangle (box car)
weighingFactor = 1.0;
break;
case FFT_WIN_TYP_HAMMING: // hamming
weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio));
break;
case FFT_WIN_TYP_HANN: // hann
weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio));
break;
case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne);
break;
case FFT_WIN_TYP_NUTTALL: // nuttall
weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) + (0.144232 * (cos(fourPi * ratio))) - (0.012604 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN: // blackman
weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall
weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) + (0.1365995 * (cos(fourPi * ratio))) - (0.0106411 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris
weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) + (0.14128 * (cos(fourPi * ratio))) - (0.01168 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_FLT_TOP: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) + (0.1980399 * cos(fourPi * ratio));
break;
case FFT_WIN_TYP_WELCH: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0));
break;
}
if (dir == FFT_FORWARD) {
vData[i] *= weighingFactor;
vData[samples - (i + 1)] *= weighingFactor;
}
else {
vData[i] /= weighingFactor;
vData[samples - (i + 1)] /= weighingFactor;
}
}
}
double arduinoFFT::MajorPeak()
{
double maxY = 0;
uint16_t IndexOfMaxY = 0;
//If sampling_frequency = 2 * max_frequency in signal,
//value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
if ((this->_vReal[i-1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i+1])) {
if (this->_vReal[i] > maxY) {
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
double delta = 0.5 * ((this->_vReal[IndexOfMaxY-1] - this->_vReal[IndexOfMaxY+1]) / (this->_vReal[IndexOfMaxY-1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY+1]));
double interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples-1);
if(IndexOfMaxY==(this->_samples >> 1)) //To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples);
// returned value: interpolated frequency peak apex
return(interpolatedX);
}
void arduinoFFT::MajorPeak(double *f, double *v)
{
double maxY = 0;
uint16_t IndexOfMaxY = 0;
//If sampling_frequency = 2 * max_frequency in signal,
//value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
if ((this->_vReal[i - 1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i + 1])) {
if (this->_vReal[i] > maxY) {
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
double delta = 0.5 * ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) / (this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]));
double interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1);
if (IndexOfMaxY == (this->_samples >> 1)) //To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples);
// returned value: interpolated frequency peak apex
*f = interpolatedX;
*v = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]);
}
double arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequency)
{
#warning("This method is deprecated and may be removed on future revisions.")
double maxY = 0;
uint16_t IndexOfMaxY = 0;
//If sampling_frequency = 2 * max_frequency in signal,
//value would be stored at position samples/2
for (uint16_t i = 1; i < ((samples >> 1) + 1); i++) {
if ((vD[i-1] < vD[i]) && (vD[i] > vD[i+1])) {
if (vD[i] > maxY) {
maxY = vD[i];
IndexOfMaxY = i;
}
}
}
double delta = 0.5 * ((vD[IndexOfMaxY-1] - vD[IndexOfMaxY+1]) / (vD[IndexOfMaxY-1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY+1]));
double interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples-1);
if(IndexOfMaxY==(samples >> 1)) //To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples);
// returned value: interpolated frequency peak apex
return(interpolatedX);
}
void arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v)
{
#warning("This method is deprecated and may be removed on future revisions.")
double maxY = 0;
uint16_t IndexOfMaxY = 0;
//If sampling_frequency = 2 * max_frequency in signal,
//value would be stored at position samples/2
for (uint16_t i = 1; i < ((samples >> 1) + 1); i++) {
if ((vD[i - 1] < vD[i]) && (vD[i] > vD[i + 1])) {
if (vD[i] > maxY) {
maxY = vD[i];
IndexOfMaxY = i;
}
}
}
double delta = 0.5 * ((vD[IndexOfMaxY - 1] - vD[IndexOfMaxY + 1]) / (vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]));
double interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples - 1);
//double popo =
if (IndexOfMaxY == (samples >> 1)) //To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples);
// returned value: interpolated frequency peak apex
*f = interpolatedX;
*v = abs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]);
}
uint8_t arduinoFFT::Exponent(uint16_t value)
{
#warning("This method may not be accessible on future revisions.")
// Calculates the base 2 logarithm of a value
uint8_t result = 0;
while (((value >> result) & 1) != 1) result++;
return(result);
}
// Private functions
void arduinoFFT::Swap(double *x, double *y)
{
double temp = *x;
*x = *y;
*y = temp;
}

481
src/arduinoFFT.h
Wyświetl plik

@ -3,6 +3,7 @@
FFT libray
Copyright (C) 2010 Didier Longueville
Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@ -19,94 +20,420 @@
*/
#ifndef arduinoFFT_h /* Prevent loading library twice */
#define arduinoFFT_h
#ifndef ArduinoFFT_h /* Prevent loading library twice */
#define ArduinoFFT_h
#ifdef ARDUINO
#if ARDUINO >= 100
#include "Arduino.h"
#else
#include "WProgram.h" /* This is where the standard Arduino code lies */
#endif
#if ARDUINO >= 100
#include "Arduino.h"
#else
#include <stdlib.h>
#include <stdio.h>
#ifdef __AVR__
#include <avr/io.h>
#include <avr/pgmspace.h>
#endif
#include <math.h>
#include "defs.h"
#include "types.h"
#include "WProgram.h" /* This is where the standard Arduino code lies */
#endif
#define FFT_LIB_REV 0x14
/* Custom constants */
#define FFT_FORWARD 0x01
#define FFT_REVERSE 0x00
/* Windowing type */
#define FFT_WIN_TYP_RECTANGLE 0x00 /* rectangle (Box car) */
#define FFT_WIN_TYP_HAMMING 0x01 /* hamming */
#define FFT_WIN_TYP_HANN 0x02 /* hann */
#define FFT_WIN_TYP_TRIANGLE 0x03 /* triangle (Bartlett) */
#define FFT_WIN_TYP_NUTTALL 0x04 /* nuttall */
#define FFT_WIN_TYP_BLACKMAN 0x05 /* blackman */
#define FFT_WIN_TYP_BLACKMAN_NUTTALL 0x06 /* blackman nuttall */
#define FFT_WIN_TYP_BLACKMAN_HARRIS 0x07 /* blackman harris*/
#define FFT_WIN_TYP_FLT_TOP 0x08 /* flat top */
#define FFT_WIN_TYP_WELCH 0x09 /* welch */
/*Mathematial constants*/
#define twoPi 6.28318531
#define fourPi 12.56637061
#define sixPi 18.84955593
#else
#include <stdlib.h>
#include <stdio.h>
#ifdef __AVR__
static const double _c1[]PROGMEM = {0.0000000000, 0.7071067812, 0.9238795325, 0.9807852804,
0.9951847267, 0.9987954562, 0.9996988187, 0.9999247018,
0.9999811753, 0.9999952938, 0.9999988235, 0.9999997059,
0.9999999265, 0.9999999816, 0.9999999954, 0.9999999989,
0.9999999997};
static const double _c2[]PROGMEM = {1.0000000000, 0.7071067812, 0.3826834324, 0.1950903220,
0.0980171403, 0.0490676743, 0.0245412285, 0.0122715383,
0.0061358846, 0.0030679568, 0.0015339802, 0.0007669903,
0.0003834952, 0.0001917476, 0.0000958738, 0.0000479369,
0.0000239684};
#include <avr/io.h>
#include <avr/pgmspace.h>
#endif
class arduinoFFT {
#include <math.h>
#include "defs.h"
#include "types.h"
#endif
// Define this to use reciprocal multiplication for division and some more speedups that might decrease precision
//#define FFT_SPEED_OVER_PRECISION
#ifndef FFT_SQRT_APPROXIMATION
#define sqrt_internal sqrt
#endif
// Define this to use a low-precision square root approximation instead of the regular sqrt() call
// This might only work for specific use cases, but is significantly faster. Only works for ArduinoFFT<float>.
//#define FFT_SQRT_APPROXIMATION
enum class FFTDirection
{
Reverse,
Forward
};
enum class FFTWindow
{
Rectangle, // rectangle (Box car)
Hamming, // hamming
Hann, // hann
Triangle, // triangle (Bartlett)
Nuttall, // nuttall
Blackman, //blackman
Blackman_Nuttall, // blackman nuttall
Blackman_Harris, // blackman harris
Flat_top, // flat top
Welch // welch
};
template <typename T>
class ArduinoFFT
{
public:
/* Constructor */
arduinoFFT(void);
arduinoFFT(double *vReal, double *vImag, uint16_t samples, double samplingFrequency);
/* Destructor */
~arduinoFFT(void);
/* Functions */
uint8_t Revision(void);
uint8_t Exponent(uint16_t value);
// Constructor
ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, T * windowWeighingFactors = nullptr)
: _vReal(vReal)
, _vImag(vImag)
, _samples(samples)
#ifdef FFT_SPEED_OVER_PRECISION
, _oneOverSamples(1.0 / samples)
#endif
, _samplingFrequency(samplingFrequency)
, _windowWeighingFactors(windowWeighingFactors)
{
// Calculates the base 2 logarithm of sample count
_power = 0;
while (((samples >> _power) & 1) != 1)
{
_power++;
}
}
void ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples);
void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t dir);
void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir);
void DCRemoval(double *vData, uint16_t samples);
double MajorPeak(double *vD, uint16_t samples, double samplingFrequency);
void MajorPeak(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v);
void Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir);
// Destructor
~ArduinoFFT()
{
}
void ComplexToMagnitude();
void Compute(uint8_t dir);
void DCRemoval();
double MajorPeak();
void MajorPeak(double *f, double *v);
void Windowing(uint8_t windowType, uint8_t dir);
// Get library revision
static uint8_t revision()
{
return 0x19;
}
// Computes in-place complex-to-complex FFT
void compute(FFTDirection dir) const
{
// Reverse bits /
uint_fast16_t j = 0;
for (uint_fast16_t i = 0; i < (this->_samples - 1); i++)
{
if (i < j)
{
Swap(this->_vReal[i], this->_vReal[j]);
if (dir == FFTDirection::Reverse)
{
Swap(this->_vImag[i], this->_vImag[j]);
}
}
uint_fast16_t k = (this->_samples >> 1);
while (k <= j)
{
j -= k;
k >>= 1;
}
j += k;
}
// Compute the FFT
#ifdef __AVR__
small_type index = 0;
#endif
T c1 = -1.0;
T c2 = 0.0;
uint_fast16_t l2 = 1;
for (uint_fast8_t l = 0; (l < this->_power); l++)
{
uint_fast16_t l1 = l2;
l2 <<= 1;
T u1 = 1.0;
T u2 = 0.0;
for (j = 0; j < l1; j++)
{
for (uint_fast16_t i = j; i < this->_samples; i += l2)
{
uint_fast16_t i1 = i + l1;
T t1 = u1 * this->_vReal[i1] - u2 * this->_vImag[i1];
T t2 = u1 * this->_vImag[i1] + u2 * this->_vReal[i1];
this->_vReal[i1] = this->_vReal[i] - t1;
this->_vImag[i1] = this->_vImag[i] - t2;
this->_vReal[i] += t1;
this->_vImag[i] += t2;
}
T z = ((u1 * c1) - (u2 * c2));
u2 = ((u1 * c2) + (u2 * c1));
u1 = z;
}
#ifdef __AVR__
c2 = pgm_read_T_near(&(_c2[index]));
c1 = pgm_read_T_near(&(_c1[index]));
index++;
#else
T cTemp = 0.5 * c1;
c2 = sqrt_internal(0.5 - cTemp);
c1 = sqrt_internal(0.5 + cTemp);
#endif
c2 = dir == FFTDirection::Forward ? -c2 : c2;
}
// Scaling for reverse transform
if (dir != FFTDirection::Forward)
{
for (uint_fast16_t i = 0; i < this->_samples; i++)
{
#ifdef FFT_SPEED_OVER_PRECISION
this->_vReal[i] *= _oneOverSamples;
this->_vImag[i] *= _oneOverSamples;
#else
this->_vReal[i] /= this->_samples;
this->_vImag[i] /= this->_samples;
#endif
}
}
}
void complexToMagnitude() const
{
// vM is half the size of vReal and vImag
for (uint_fast16_t i = 0; i < this->_samples; i++)
{
this->_vReal[i] = sqrt_internal(sq(this->_vReal[i]) + sq(this->_vImag[i]));
}
}
void dcRemoval() const
{
// calculate the mean of vData
T mean = 0;
for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
mean += this->_vReal[i];
}
mean /= this->_samples;
// Subtract the mean from vData
for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
this->_vReal[i] -= mean;
}
}
void windowing(FFTWindow windowType, FFTDirection dir)
{
// check if values are already pre-computed for the correct window type
if (_windowWeighingFactors && weighingFactorsComputed && weighingFactorsFFTWindow == windowType)
{
// yes. values are precomputed
if (dir == FFTDirection::Forward)
{
for (uint_fast16_t i = 0; i < (this->_samples >> 1); i++)
{
this->_vReal[i] *= _windowWeighingFactors[i];
this->_vReal[this->_samples - (i + 1)] *= _windowWeighingFactors[i];
}
}
else
{
for (uint_fast16_t i = 0; i < (this->_samples >> 1); i++)
{
#ifdef FFT_SPEED_OVER_PRECISION
// on many architectures reciprocals and multiplying are much faster than division
T oneOverFactor = 1.0 / _windowWeighingFactors[i];
this->_vReal[i] *= oneOverFactor;
this->_vReal[this->_samples - (i + 1)] *= oneOverFactor;
#else
this->_vReal[i] /= _windowWeighingFactors[i];
this->_vReal[this->_samples - (i + 1)] /= _windowWeighingFactors[i];
#endif
}
}
}
else
{
// no. values need to be pre-computed or applied
T samplesMinusOne = (T(this->_samples) - 1.0);
for (uint_fast16_t i = 0; i < (this->_samples >> 1); i++)
{
T indexMinusOne = T(i);
T ratio = (indexMinusOne / samplesMinusOne);
T weighingFactor = 1.0;
// Compute and record weighting factor
switch (windowType)
{
case FFTWindow::Rectangle: // rectangle (box car)
weighingFactor = 1.0;
break;
case FFTWindow::Hamming: // hamming
weighingFactor = 0.54 - (0.46 * cos(TWO_PI * ratio));
break;
case FFTWindow::Hann: // hann
weighingFactor = 0.54 * (1.0 - cos(TWO_PI * ratio));
break;
case FFTWindow::Triangle: // triangle (Bartlett)
weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne);
break;
case FFTWindow::Nuttall: // nuttall
weighingFactor = 0.355768 - (0.487396 * (cos(TWO_PI * ratio))) + (0.144232 * (cos(FOUR_PI * ratio))) - (0.012604 * (cos(SIX_PI * ratio)));
break;
case FFTWindow::Blackman: // blackman
weighingFactor = 0.42323 - (0.49755 * (cos(TWO_PI * ratio))) + (0.07922 * (cos(FOUR_PI * ratio)));
break;
case FFTWindow::Blackman_Nuttall: // blackman nuttall
weighingFactor = 0.3635819 - (0.4891775 * (cos(TWO_PI * ratio))) + (0.1365995 * (cos(FOUR_PI * ratio))) - (0.0106411 * (cos(SIX_PI * ratio)));
break;
case FFTWindow::Blackman_Harris: // blackman harris
weighingFactor = 0.35875 - (0.48829 * (cos(TWO_PI * ratio))) + (0.14128 * (cos(FOUR_PI * ratio))) - (0.01168 * (cos(SIX_PI * ratio)));
break;
case FFTWindow::Flat_top: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(TWO_PI * ratio)) + (0.1980399 * cos(FOUR_PI * ratio));
break;
case FFTWindow::Welch: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0));
break;
}
if (_windowWeighingFactors)
{
_windowWeighingFactors[i] = weighingFactor;
}
if (dir == FFTDirection::Forward)
{
this->_vReal[i] *= weighingFactor;
this->_vReal[this->_samples - (i + 1)] *= weighingFactor;
}
else
{
#ifdef FFT_SPEED_OVER_PRECISION
// on many architectures reciprocals and multiplying are much faster than division
T oneOverFactor = 1.0 / weighingFactor;
this->_vReal[i] *= oneOverFactor;
this->_vReal[this->_samples - (i + 1)] *= oneOverFactor;
#else
this->_vReal[i] /= weighingFactor;
this->_vReal[this->_samples - (i + 1)] /= weighingFactor;
#endif
}
}
// mark cached values as pre-computed
weighingFactorsFFTWindow = windowType;
weighingFactorsComputed = true;
}
}
T majorPeak() const
{
T maxY = 0;
uint_fast16_t IndexOfMaxY = 0;
//If sampling_frequency = 2 * max_frequency in signal,
//value would be stored at position samples/2
for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
if ((this->_vReal[i - 1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i + 1]))
{
if (this->_vReal[i] > maxY)
{
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
T delta = 0.5 * ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) / (this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]));
T interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1);
if (IndexOfMaxY == (this->_samples >> 1))
{
//To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples);
}
// returned value: interpolated frequency peak apex
return interpolatedX;
}
void majorPeak(T &f, T &v) const
{
T maxY = 0;
uint_fast16_t IndexOfMaxY = 0;
//If sampling_frequency = 2 * max_frequency in signal,
//value would be stored at position samples/2
for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
if ((this->_vReal[i - 1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i + 1]))
{
if (this->_vReal[i] > maxY)
{
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
T delta = 0.5 * ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) / (this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]));
T interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1);
if (IndexOfMaxY == (this->_samples >> 1))
{
//To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples);
}
// returned value: interpolated frequency peak apex
f = interpolatedX;
v = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]);
}
private:
#ifdef __AVR__
static const T _c1[] PROGMEM = {
0.0000000000, 0.7071067812, 0.9238795325, 0.9807852804,
0.9951847267, 0.9987954562, 0.9996988187, 0.9999247018,
0.9999811753, 0.9999952938, 0.9999988235, 0.9999997059,
0.9999999265, 0.9999999816, 0.9999999954, 0.9999999989,
0.9999999997};
static const T _c2[] PROGMEM = {
1.0000000000, 0.7071067812, 0.3826834324, 0.1950903220,
0.0980171403, 0.0490676743, 0.0245412285, 0.0122715383,
0.0061358846, 0.0030679568, 0.0015339802, 0.0007669903,
0.0003834952, 0.0001917476, 0.0000958738, 0.0000479369,
0.0000239684};
#endif
// Mathematial constants
#ifndef TWO_PI
static constexpr T TWO_PI = 6.28318531; // might already be defined in Arduino.h
#endif
static constexpr T FOUR_PI = 12.56637061;
static constexpr T SIX_PI = 18.84955593;
static inline void Swap(T &x, T &y)
{
T temp = x;
x = y;
y = temp;
}
#ifdef FFT_SQRT_APPROXIMATION
template<typename V = T>
static inline V sqrt_internal(typename std::enable_if<std::is_same<V, float>::value, V>::type x)
{
union {
int i;
float x;
} u;
u.x = x;
u.i = (1 << 29) + (u.i >> 1) - (1 << 22);
// Two Babylonian Steps (simplified from:)
// u.x = 0.5f * (u.x + x/u.x);
// u.x = 0.5f * (u.x + x/u.x);
u.x = u.x + x / u.x;
u.x = 0.25f * u.x + x / u.x;
return u.x;
}
template<typename V = T>
static inline V sqrt_internal(typename std::enable_if<!std::is_same<V, float>::value, V>::type x)
{
return sqrt(x);
}
#endif
/* Variables */
uint16_t _samples;
double _samplingFrequency;
double *_vReal;
double *_vImag;
uint8_t _power;
/* Functions */
void Swap(double *x, double *y);
uint_fast16_t _samples = 0;
#ifdef FFT_SPEED_OVER_PRECISION
T _oneOverSamples = 0.0;
#endif
T _samplingFrequency = 0;
T *_vReal = nullptr;
T *_vImag = nullptr;
T * _windowWeighingFactors = nullptr;
FFTWindow weighingFactorsFFTWindow;
bool weighingFactorsComputed = false;
uint_fast8_t _power = 0;
};
#endif