From cb33149c173ebe8863d511aeba6c7ea69db544de Mon Sep 17 00:00:00 2001 From: Bim Overbohm Date: Wed, 19 Feb 2020 17:15:49 +0100 Subject: [PATCH] Remove deprecated functions, templatize, Speedup --- Examples/FFT_01/FFT_01.ino | 19 +- Examples/FFT_02/FFT_02.ino | 21 +- Examples/FFT_03/FFT_03.ino | 17 +- Examples/FFT_04/FFT_04.ino | 18 +- Examples/FFT_05/FFT_05.ino | 19 +- Examples/FFT_speedup/FFT_speedup.ino | 129 +++++++ README.md | 98 +++--- changeLog.txt | 6 + keywords.txt | 44 +-- library.json | 7 +- library.properties | 2 +- src/arduinoFFT.cpp | 446 ------------------------- src/arduinoFFT.h | 481 ++++++++++++++++++++++----- 13 files changed, 673 insertions(+), 634 deletions(-) create mode 100644 Examples/FFT_speedup/FFT_speedup.ino delete mode 100644 src/arduinoFFT.cpp diff --git a/Examples/FFT_01/FFT_01.ino b/Examples/FFT_01/FFT_01.ino index caefe80..22b5024 100644 --- a/Examples/FFT_01/FFT_01.ino +++ b/Examples/FFT_01/FFT_01.ino @@ -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 FFT = ArduinoFFT(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 */ diff --git a/Examples/FFT_02/FFT_02.ino b/Examples/FFT_02/FFT_02.ino index cb97e29..bfa1804 100644 --- a/Examples/FFT_02/FFT_02.ino +++ b/Examples/FFT_02/FFT_02.ino @@ -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 FFT = ArduinoFFT(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); diff --git a/Examples/FFT_03/FFT_03.ino b/Examples/FFT_03/FFT_03.ino index 9e1640c..2e50613 100644 --- a/Examples/FFT_03/FFT_03.ino +++ b/Examples/FFT_03/FFT_03.ino @@ -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 FFT = ArduinoFFT(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 */ diff --git a/Examples/FFT_04/FFT_04.ino b/Examples/FFT_04/FFT_04.ino index 5c1fcd5..b125991 100644 --- a/Examples/FFT_04/FFT_04.ino +++ b/Examples/FFT_04/FFT_04.ino @@ -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 FFT = ArduinoFFT(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 */ } diff --git a/Examples/FFT_05/FFT_05.ino b/Examples/FFT_05/FFT_05.ino index 7abe3b6..a6f4df7 100644 --- a/Examples/FFT_05/FFT_05.ino +++ b/Examples/FFT_05/FFT_05.ino @@ -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 FFT = ArduinoFFT(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); diff --git a/Examples/FFT_speedup/FFT_speedup.ino b/Examples/FFT_speedup/FFT_speedup.ino new file mode 100644 index 0000000..a059a17 --- /dev/null +++ b/Examples/FFT_speedup/FFT_speedup.ino @@ -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 . + +*/ + +// 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. +//#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 FFT = ArduinoFFT(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> 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(); +} diff --git a/README.md b/README.md index 37861ba..feb4d7a 100644 --- a/README.md +++ b/README.md @@ -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. -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 - -### 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 ` -### TODO +## TODO * Ratio table for windowing function. * Document windowing functions advantages and disadvantages. * Optimize usage and arguments. * Add new windowing functions. -* Spectrum table? +* ~~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`. diff --git a/changeLog.txt b/changeLog.txt index e2cb0fb..b6dc9ad 100644 --- a/changeLog.txt +++ b/changeLog.txt @@ -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%. diff --git a/keywords.txt b/keywords.txt index 1daabde..c0f1804 100644 --- a/keywords.txt +++ b/keywords.txt @@ -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 diff --git a/library.json b/library.json index 32c7606..e3b8371 100644 --- a/library.json +++ b/library.json @@ -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": "*" } diff --git a/library.properties b/library.properties index cd86480..986fbe3 100644 --- a/library.properties +++ b/library.properties @@ -1,5 +1,5 @@ name=arduinoFFT -version=1.5.5 +version=1.9.0 author=Enrique Condes maintainer=Enrique Condes sentence=A library for implementing floating point Fast Fourier Transform calculations on Arduino. diff --git a/src/arduinoFFT.cpp b/src/arduinoFFT.cpp deleted file mode 100644 index f552fc9..0000000 --- a/src/arduinoFFT.cpp +++ /dev/null @@ -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 . - -*/ - -#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; -} diff --git a/src/arduinoFFT.h b/src/arduinoFFT.h index 73927dc..bb67f91 100644 --- a/src/arduinoFFT.h +++ b/src/arduinoFFT.h @@ -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 - #include - #ifdef __AVR__ - #include - #include - #endif - #include - #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 +#include #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 +#include #endif -class arduinoFFT { +#include +#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. +//#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 +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 + static inline V sqrt_internal(typename std::enable_if::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 + static inline V sqrt_internal(typename std::enable_if::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