diff --git a/Examples/FFT_01/FFT_01.ino b/Examples/FFT_01/FFT_01.ino index c41d22c..1f6b3da 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 library - Copyright (C) 2014 Enrique Condes + Example of use of the FFT libray + + Copyright (C) 2014 Enrique Condes + Copyright (C) 2020 Bim Overbohm (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; /* 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 @@ -67,23 +72,21 @@ void loop() //vReal[i] = uint8_t((amplitude * (sin(i * ratio) + 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 = arduinoFFT(vReal, vImag, samples, samplingFrequency); /* Create FFT object */ /* Print the results of the simulated sampling according to time */ Serial.println("Data:"); PrintVector(vReal, samples, SCL_TIME); - FFT.Windowing(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(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(); /* Compute magnitudes */ + FFT.complexToMagnitude(); /* Compute magnitudes */ Serial.println("Computed magnitudes:"); PrintVector(vReal, (samples >> 1), SCL_FREQUENCY); - double x = FFT.MajorPeak(); + 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 4c7ecd6..ebe8d5e 100644 --- a/Examples/FFT_02/FFT_02.ino +++ b/Examples/FFT_02/FFT_02.ino @@ -1,10 +1,12 @@ /* - Example of use of the FFT library to compute FFT for several signals over a range of frequencies. - The exponent is calculated once before the execution since it is a constant. - This saves resources during the execution of the sketch and reduces the compiled size. - The sketch shows the time that the computation takes. - Copyright (C) 2014 Enrique Condes + Example of use of the FFT libray to compute FFT for several signals over a range of frequencies. + 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) 2020 Bim Overbohm (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,11 +25,9 @@ #include "arduinoFFT.h" -arduinoFFT FFT; /* 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; @@ -42,7 +42,10 @@ Input vectors receive computed results from FFT double vReal[samples]; double vImag[samples]; -unsigned long time; +/* Create FFT object */ +ArduinoFFT FFT = ArduinoFFT(vReal, vImag, samples, sampling); + +unsigned long startTime; #define SCL_INDEX 0x00 #define SCL_TIME 0x01 @@ -71,25 +74,24 @@ void loop() } /*Serial.println("Data:"); PrintVector(vReal, samples, SCL_TIME);*/ - time=millis(); - FFT = arduinoFFT(vReal, vImag, samples, sampling); /* Create FFT object */ - FFT.Windowing(FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */ + startTime=millis(); + FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */ /*Serial.println("Weighed data:"); PrintVector(vReal, samples, SCL_TIME);*/ - FFT.Compute(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(); /* Compute magnitudes */ + FFT.complexToMagnitude(); /* Compute magnitudes */ /*Serial.println("Computed magnitudes:"); PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);*/ - double x = FFT.MajorPeak(); + double x = FFT.majorPeak(); Serial.print(frequency); Serial.print(": \t\t"); Serial.print(x, 4); Serial.print("\t\t"); - Serial.print(millis()-time); + Serial.print(millis()-startTime); Serial.println(" ms"); // delay(2000); /* Repeat after delay */ } diff --git a/Examples/FFT_03/FFT_03.ino b/Examples/FFT_03/FFT_03.ino index ac3fa12..0a76b34 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 library to compute FFT for a signal sampled through the ADC. - Copyright (C) 2018 Enrique Condés and Ragnar Ranøyen Homb + 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) 2020 Bim Overbohm (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; /* 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 @@ -64,22 +67,21 @@ void loop() } microseconds += sampling_period_us; } - FFT = arduinoFFT(vReal, vImag, samples, samplingFrequency); /* Create FFT object */ /* Print the results of the sampling according to time */ Serial.println("Data:"); PrintVector(vReal, samples, SCL_TIME); - FFT.Windowing(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(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(); /* Compute magnitudes */ + FFT.complexToMagnitude(); /* Compute magnitudes */ Serial.println("Computed magnitudes:"); PrintVector(vReal, (samples >> 1), SCL_FREQUENCY); - double x = FFT.MajorPeak(); + 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 be6a970..49ac3f6 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 library - Copyright (C) 2018 Enrique Condes + Example of use of the FFT libray + + Copyright (C) 2018 Enrique Condes + Copyright (C) 2020 Bim Overbohm (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; /* 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 @@ -68,12 +72,11 @@ void loop() //vReal[i] = uint8_t((amplitude * (sin(i * ratio) + 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 = arduinoFFT(vReal, vImag, samples, samplingFrequency); /* Create FFT object */ - FFT.Windowing(FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */ - FFT.Compute(FFT_FORWARD); /* Compute FFT */ - FFT.ComplexToMagnitude(); /* 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(); + 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 a382b7f..ca3152e 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 library - Copyright (C) 2014 Enrique Condes + Example of use of the FFT libray + + Copyright (C) 2014 Enrique Condes + Copyright (C) 2020 Bim Overbohm (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; /* 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 @@ -68,24 +73,23 @@ void loop() //vReal[i] = uint8_t((amplitude * (sin(i * ratio) + 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 = arduinoFFT(vReal, vImag, samples, samplingFrequency); /* Create FFT object */ /* Print the results of the simulated sampling according to time */ Serial.println("Data:"); PrintVector(vReal, samples, SCL_TIME); - FFT.Windowing(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(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(); /* Compute magnitudes */ + FFT.complexToMagnitude(); /* Compute magnitudes */ Serial.println("Computed magnitudes:"); PrintVector(vReal, (samples >> 1), SCL_FREQUENCY); double x; double v; - FFT.MajorPeak(&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..e91a2fe --- /dev/null +++ b/Examples/FFT_speedup/FFT_speedup.ino @@ -0,0 +1,123 @@ +/* + + 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 (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]; + +/* Create FFT object with weighing factor storage */ +ArduinoFFT FFT = ArduinoFFT(vReal, vImag, samples, samplingFrequency, true); + +#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 93cd705..dbafd42 100644 --- a/README.md +++ b/README.md @@ -1,87 +1,38 @@ 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. +This is version 2.0 of the library, which has a different [API](#api). -Tested on Arduino 1.6.11 +Tested on Arduino 1.8.19 and 2.3.2. -### 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\src\arduinoFTT.h` (the library header file. include this in your project) +`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. - To use this library in a sketch, go to the Sketch | Import Library menu and select arduinoFTT. This will add a corresponding line to the top of your sketch: `#include ` -### TODO -* Ratio table for windowing function. -* Document windowing functions advantages and disadvantages. -* Optimize usage and arguments. -* Add new windowing functions. -* Spectrum table? +## API -### API -The exclamation mark `!` denotes that this method is deprecated and may be removed on future revisions. - -* **!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); -Calculates the Fast Fourier Transform. -* **!DCRemoval**(double *vData, uint16_t samples); -* **DCRemoval**(); -Removes the DC component from the sample data. -* **!MajorPeak**(double *vD, uint16_t samples, double samplingFrequency); -* **!MajorPeak**(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v); -* **MajorPeak**(); -* **MajorPeak**(double *f, double *v); -* **MajorPeakParabola**(); -Looks for and returns the frequency of the biggest spike in the analyzed signal. -* **Revision**(void); -Returns the library revision. -* **!Windowing**(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir); -* **Windowing**(uint8_t windowType, uint8_t 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. +Documentation was moved to the project's [wiki](https://github.com/kosme/arduinoFFT/wiki). diff --git a/changeLog.txt b/changeLog.txt deleted file mode 100644 index e15375d..0000000 --- a/changeLog.txt +++ /dev/null @@ -1,35 +0,0 @@ -28/12/23 v1.6.2 -Fix issue 52 and symplify calculation of synthetic data for examples 1, 2, 4, and 5. - -25/07/23 v1.6.1 -Code optimization for speed improvements. See issue 84 for more details. - -03/09/23 v1.6 -Include some functionality from development branch. - -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%. - -02/10/18 v1.4 -Transition version. Minor optimization to functions. New API. Deprecation of old functions. - -12/06/18 v1.3 -Add support for mbed development boards. - -09/04/17 v1.2.3 -Finally solves the issue of Arduino IDE not correctly detecting and highlighting the keywords. - -09/03/17 v1.2.2 -Solves a format issue in keywords.txt that prevented keywords from being detected. - -08/28/17 v1.2.1 -Fix to issues 6 and 7. Not cleaning the imaginary vector after each cycle leaded to erroneous calculations and could cause buffer overflows. - -08/04/17 v1.2 -Fix to bug preventing the number of samples to be greater than 128. New logical limit is 32768 samples but it is bound to the RAM on the chip. - -05/12/17 v1.1 -Fix issue that prevented installation through the Arduino Library Manager interface. - -05/11/17 v1.0 -Initial commit to Arduino Library Manager. diff --git a/keywords.txt b/keywords.txt index 0643678..c9b49c1 100644 --- a/keywords.txt +++ b/keywords.txt @@ -6,36 +6,37 @@ # 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 -MajorPeakParabola KEYWORD2 +complexToMagnitude KEYWORD2 +compute KEYWORD2 +dcRemoval KEYWORD2 +majorPeak KEYWORD2 +majorPeakParabola KEYWORD2 +revision KEYWORD2 +setArrays KEYWORD2 +windowing 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 + +Blackman LITERAL1 +Blackman_Harris LITERAL1 +Blackman_Nuttall LITERAL1 +Flat_top LITERAL1 +Hamming LITERAL1 +Hann LITERAL1 +Nuttall LITERAL1 +Rectangle LITERAL1 +Triangle LITERAL1 +Welch LITERAL1 diff --git a/library.json b/library.json index 87dadea..c9e5fe9 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.6.2", + "version": "2.0", "frameworks": ["arduino","mbed","espidf"], "platforms": "*", "headers": "arduinoFFT.h" diff --git a/library.properties b/library.properties index 951628c..5826255 100644 --- a/library.properties +++ b/library.properties @@ -1,5 +1,5 @@ name=arduinoFFT -version=1.6.2 +version=2.0 author=Enrique Condes maintainer=Enrique Condes sentence=A library for implementing floating point Fast Fourier Transform calculations on Arduino. diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..00e95bf --- /dev/null +++ b/src/.gitignore @@ -0,0 +1 @@ +/arduinoFFT.h.gch diff --git a/src/arduinoFFT.cpp b/src/arduinoFFT.cpp index 02b1b67..1bb5b22 100644 --- a/src/arduinoFFT.cpp +++ b/src/arduinoFFT.cpp @@ -1,529 +1,448 @@ /* + FFT library + Copyright (C) 2010 Didier Longueville + Copyright (C) 2014 Enrique Condes + Copyright (C) 2020 Bim Overbohm (template, speed improvements) - FFT library - 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 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. - 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 . + 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.") +template ArduinoFFT::ArduinoFFT() {} + +template +ArduinoFFT::ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, + T samplingFrequency, bool windowingFactors) + : _samples(samples), _samplingFrequency(samplingFrequency), _vImag(vImag), + _vReal(vReal) { + if (windowingFactors) { + _precompiledWindowingFactors = new T[samples / 2]; + } + _power = exponent(samples); +#ifdef FFT_SPEED_OVER_PRECISION + _oneOverSamples = 1.0 / samples; +#endif } -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) { +template ArduinoFFT::~ArduinoFFT(void) { // Destructor -} - -uint8_t arduinoFFT::Revision(void) { return (FFT_LIB_REV); } - -void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, - FFTDirection dir) { -#warning("This method is deprecated and may be removed on future revisions.") - Compute(vReal, vImag, samples, Exponent(samples), dir); -} - -void arduinoFFT::Compute(FFTDirection 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) { - double reciprocal = 1.0 / this->_samples; - for (uint16_t i = 0; i < this->_samples; i++) { - this->_vReal[i] *= reciprocal; - this->_vImag[i] *= reciprocal; - } + if (_precompiledWindowingFactors) { + delete [] _precompiledWindowingFactors; } } -void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, - uint8_t power, FFTDirection dir) { +template void ArduinoFFT::complexToMagnitude(void) const { + complexToMagnitude(this->_vReal, this->_vImag, this->_samples); +} + +template +void ArduinoFFT::complexToMagnitude(T *vReal, T *vImag, + uint_fast16_t samples) const { + // vM is half the size of vReal and vImag + for (uint_fast16_t i = 0; i < samples; i++) { + vReal[i] = sqrt_internal(sq(vReal[i]) + sq(vImag[i])); + } +} + +template void ArduinoFFT::compute(FFTDirection dir) const { + compute(this->_vReal, this->_vImag, this->_samples, exponent(this->_samples), + dir); +} + +template +void ArduinoFFT::compute(T *vReal, T *vImag, uint_fast16_t samples, + FFTDirection dir) const { + compute(vReal, vImag, samples, exponent(samples), 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++) { +template +void ArduinoFFT::compute(T *vReal, T *vImag, uint_fast16_t samples, + uint_fast8_t power, FFTDirection dir) const { +#ifdef FFT_SPEED_OVER_PRECISION + T oneOverSamples = this->_oneOverSamples; + if (!this->_oneOverSamples) + oneOverSamples = 1.0 / samples; +#endif + // Reverse bits + uint_fast16_t j = 0; + for (uint_fast16_t i = 0; i < (samples - 1); i++) { if (i < j) { - Swap(&vReal[i], &vReal[j]); - if (dir == FFT_REVERSE) - Swap(&vImag[i], &vImag[j]); + swap(&vReal[i], &vReal[j]); + if (dir == FFTDirection::Reverse) + swap(&vImag[i], &vImag[j]); } - uint16_t k = (samples >> 1); + uint_fast16_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; + // Compute the FFT + T c1 = -1.0; + T c2 = 0.0; + uint_fast16_t l2 = 1; + for (uint_fast8_t l = 0; (l < power); l++) { + uint_fast16_t l1 = l2; l2 <<= 1; - double u1 = 1.0; - double u2 = 0.0; + T u1 = 1.0; + T 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]; + for (uint_fast16_t i = j; i < samples; i += l2) { + uint_fast16_t i1 = i + l1; + T t1 = u1 * vReal[i1] - u2 * vImag[i1]; + T 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)); + T 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++; + +#if defined(__AVR__) && defined(USE_AVR_PROGMEM) + c2 = pgm_read_float_near(&(_c2[l])); + c1 = pgm_read_float_near(&(_c1[l])); #else - c2 = sqrt((1.0 - c1) / 2.0); - c1 = sqrt((1.0 + c1) / 2.0); + T cTemp = 0.5 * c1; + c2 = sqrt_internal(0.5 - cTemp); + c1 = sqrt_internal(0.5 + cTemp); #endif - if (dir == FFT_FORWARD) { + + if (dir == FFTDirection::Forward) { c2 = -c2; } } // Scaling for reverse transform - if (dir != FFT_FORWARD) { - double reciprocal = 1.0 / samples; - for (uint16_t i = 0; i < samples; i++) { - vReal[i] *= reciprocal; - vImag[i] *= reciprocal; + if (dir == FFTDirection::Reverse) { + for (uint_fast16_t i = 0; i < samples; i++) { +#ifdef FFT_SPEED_OVER_PRECISION + vReal[i] *= oneOverSamples; + vImag[i] *= oneOverSamples; +#else + vReal[i] /= samples; + vImag[i] /= samples; +#endif } } } -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])); - } +template void ArduinoFFT::dcRemoval(void) const { + dcRemoval(this->_vReal, this->_samples); } -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() { +template +void ArduinoFFT::dcRemoval(T *vData, uint_fast16_t samples) const { // calculate the mean of vData - double mean = 0; - for (uint16_t i = 0; i < this->_samples; i++) { - mean += this->_vReal[i]; - } - mean /= this->_samples; - // Subtract the mean from vData - for (uint16_t i = 0; i < this->_samples; 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 = 0; i < samples; i++) { + T mean = 0; + for (uint_fast16_t i = 0; i < samples; i++) { mean += vData[i]; } mean /= samples; // Subtract the mean from vData - for (uint16_t i = 0; i < samples; i++) { + for (uint_fast16_t i = 0; i < samples; i++) { vData[i] -= mean; } } -void arduinoFFT::Windowing(FFTWindow windowType, FFTDirection dir) { - // Weighing factors are computed once before multiple use of FFT - // The weighing function is symmetric; 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) -#if defined(ESP8266) || defined(ESP32) - weighingFactor = - 1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) / - samplesMinusOne); -#else - weighingFactor = - 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / - samplesMinusOne); -#endif - 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; - } - } +template T ArduinoFFT::majorPeak(void) const { + return majorPeak(this->_vReal, this->_samples, this->_samplingFrequency); } -void arduinoFFT::Windowing(double *vData, uint16_t samples, - FFTWindow windowType, FFTDirection dir) { -// Weighing factors are computed once before multiple use of FFT -// The weighing function is symmetric; 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) -#if defined(ESP8266) || defined(ESP32) - weighingFactor = - 1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) / - samplesMinusOne); -#else - weighingFactor = - 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / - samplesMinusOne); -#endif - 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; - } - } +template void ArduinoFFT::majorPeak(T *f, T *v) const { + majorPeak(this->_vReal, this->_samples, this->_samplingFrequency, f, v); } -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); +template +T ArduinoFFT::majorPeak(T *vData, uint_fast16_t samples, + T samplingFrequency) const { + T frequency; + majorPeak(vData, samples, samplingFrequency, &frequency, nullptr); + return frequency; } -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; -#if defined(ESP8266) || defined(ESP32) - *v = fabs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + - this->_vReal[IndexOfMaxY + 1]); -#else - *v = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + - this->_vReal[IndexOfMaxY + 1]); -#endif -} +template +void ArduinoFFT::majorPeak(T *vData, uint_fast16_t samples, + T samplingFrequency, T *frequency, + T *magnitude) const { + T maxY = 0; + uint_fast16_t IndexOfMaxY = 0; + findMaxY(vData, (samples >> 1) + 1, &maxY, &IndexOfMaxY); -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); + T delta = 0.5 * ((vData[IndexOfMaxY - 1] - vData[IndexOfMaxY + 1]) / + (vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) + + vData[IndexOfMaxY + 1])); + T 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; + *frequency = interpolatedX; + if (magnitude != nullptr) { #if defined(ESP8266) || defined(ESP32) - *v = - fabs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]); + *magnitude = fabs(vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) + + vData[IndexOfMaxY + 1]); #else - *v = abs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]); + *magnitude = abs(vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) + + vData[IndexOfMaxY + 1]); #endif + } } -double arduinoFFT::MajorPeakParabola() { - 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 freq = 0; - if (IndexOfMaxY > 0) { - // Assume the three points to be on a parabola - double a, b, c; - Parabola(IndexOfMaxY - 1, this->_vReal[IndexOfMaxY - 1], IndexOfMaxY, - this->_vReal[IndexOfMaxY], IndexOfMaxY + 1, - this->_vReal[IndexOfMaxY + 1], &a, &b, &c); - - // Peak is at the middle of the parabola - double x = -b / (2 * a); - - // And magnitude is at the extrema of the parabola if you want It... - // double y = a*x*x+b*x+c; - - // Convert to frequency - freq = (x * this->_samplingFrequency) / (this->_samples); - } - +template T ArduinoFFT::majorPeakParabola(void) const { + T freq = 0; + majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency, + &freq, nullptr); return freq; } -void arduinoFFT::Parabola(double x1, double y1, double x2, double y2, double x3, - double y3, double *a, double *b, double *c) { - double reversed_denom = 1 / ((x1 - x2) * (x1 - x3) * (x2 - x3)); +template +void ArduinoFFT::majorPeakParabola(T *frequency, T *magnitude) const { + majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency, + frequency, magnitude); +} + +template +T ArduinoFFT::majorPeakParabola(T *vData, uint_fast16_t samples, + T samplingFrequency) const { + T freq = 0; + majorPeakParabola(vData, samples, samplingFrequency, &freq, nullptr); + return freq; +} + +template +void ArduinoFFT::majorPeakParabola(T *vData, uint_fast16_t samples, + T samplingFrequency, T *frequency, + T *magnitude) const { + T maxY = 0; + uint_fast16_t IndexOfMaxY = 0; + findMaxY(vData, (samples >> 1) + 1, &maxY, &IndexOfMaxY); + + *frequency = 0; + if (IndexOfMaxY > 0) { + // Assume the three points to be on a parabola + T a, b, c; + parabola(IndexOfMaxY - 1, vData[IndexOfMaxY - 1], IndexOfMaxY, + vData[IndexOfMaxY], IndexOfMaxY + 1, vData[IndexOfMaxY + 1], &a, + &b, &c); + + // Peak is at the middle of the parabola + T x = -b / (2 * a); + + // And magnitude is at the extrema of the parabola if you want It... + if (magnitude != nullptr) { + *magnitude = a * x * x + b * x + c; + } + + // Convert to frequency + *frequency = (x * samplingFrequency) / samples; + } +} + +template uint8_t ArduinoFFT::revision(void) { + return (FFT_LIB_REV); +} + +// Replace the data array pointers +template +void ArduinoFFT::setArrays(T *vReal, T *vImag, uint_fast16_t samples) { + _vReal = vReal; + _vImag = vImag; + if (samples) { + _samples = samples; +#ifdef FFT_SPEED_OVER_PRECISION + _oneOverSamples = 1.0 / samples; +#endif + if (_precompiledWindowingFactors) { + delete [] _precompiledWindowingFactors; + } + _precompiledWindowingFactors = new T[samples / 2]; + } +} + +template +void ArduinoFFT::windowing(FFTWindow windowType, FFTDirection dir, + bool withCompensation) { + // The windowing function is the same, precompiled values can be used, and + // precompiled values exist + if (this->_precompiledWindowingFactors && this->_isPrecompiled && + this->_windowFunction == windowType && + this->_precompiledWithCompensation == withCompensation) { + windowing(this->_vReal, this->_samples, FFTWindow::Precompiled, dir, + this->_precompiledWindowingFactors, withCompensation); + // Precompiled values must be generated. Either the function changed or the + // precompiled values don't exist + } else if (this->_precompiledWindowingFactors) { + windowing(this->_vReal, this->_samples, windowType, dir, + this->_precompiledWindowingFactors, withCompensation); + this->_isPrecompiled = true; + this->_precompiledWithCompensation = withCompensation; + this->_windowFunction = windowType; + // Don't care about precompiled windowing values + } else { + windowing(this->_vReal, this->_samples, windowType, dir, nullptr, + withCompensation); + } +} + +template +void ArduinoFFT::windowing(T *vData, uint_fast16_t samples, + FFTWindow windowType, FFTDirection dir, + T *windowingFactors, bool withCompensation) { + // Weighing factors are computed once before multiple use of FFT + // The weighing function is symmetric; half the weighs are recorded + if (windowingFactors != nullptr && windowType == FFTWindow::Precompiled) { + for (uint_fast16_t i = 0; i < (samples >> 1); i++) { + if (dir == FFTDirection::Forward) { + vData[i] *= windowingFactors[i]; + vData[samples - (i + 1)] *= windowingFactors[i]; + } else { +#ifdef FFT_SPEED_OVER_PRECISION + T inverse = 1.0 / windowingFactors[i]; + vData[i] *= inverse; + vData[samples - (i + 1)] *= inverse; +#else + vData[i] /= windowingFactors[i]; + vData[samples - (i + 1)] /= windowingFactors[i]; +#endif + } + } + } else { + T samplesMinusOne = (T(samples) - 1.0); + T compensationFactor; + if (withCompensation) { + compensationFactor = + _WindowCompensationFactors[static_cast(windowType)]; + } + for (uint_fast16_t i = 0; i < (samples >> 1); i++) { + T indexMinusOne = T(i); + T ratio = (indexMinusOne / samplesMinusOne); + T weighingFactor = 1.0; + // Compute and record weighting factor + switch (windowType) { + case FFTWindow::Hamming: // hamming + weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio)); + break; + case FFTWindow::Hann: // hann + weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio)); + break; + case FFTWindow::Triangle: // triangle (Bartlett) +#if defined(ESP8266) || defined(ESP32) + weighingFactor = + 1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) / + samplesMinusOne); +#else + weighingFactor = + 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / + samplesMinusOne); +#endif + break; + case FFTWindow::Nuttall: // nuttall + weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) + + (0.144232 * (cos(fourPi * ratio))) - + (0.012604 * (cos(sixPi * ratio))); + break; + case FFTWindow::Blackman: // blackman + weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + + (0.07922 * (cos(fourPi * ratio))); + break; + case FFTWindow::Blackman_Nuttall: // blackman nuttall + weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) + + (0.1365995 * (cos(fourPi * ratio))) - + (0.0106411 * (cos(sixPi * ratio))); + break; + case FFTWindow::Blackman_Harris: // blackman harris + weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) + + (0.14128 * (cos(fourPi * ratio))) - + (0.01168 * (cos(sixPi * ratio))); + break; + case FFTWindow::Flat_top: // flat top + weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) + + (0.1980399 * cos(fourPi * ratio)); + break; + case FFTWindow::Welch: // welch + weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / + (samplesMinusOne / 2.0)); + break; + default: + // This is Rectangle windowing which doesn't do anything + // and Precompiled which shouldn't be selected + break; + } + if (withCompensation) { + weighingFactor *= compensationFactor; + } + if (windowingFactors) { + windowingFactors[i] = weighingFactor; + } + if (dir == FFTDirection::Forward) { + vData[i] *= weighingFactor; + vData[samples - (i + 1)] *= weighingFactor; + } else { +#ifdef FFT_SPEED_OVER_PRECISION + T inverse = 1.0 / weighingFactor; + vData[i] *= inverse; + vData[samples - (i + 1)] *= inverse; +#else + vData[i] /= weighingFactor; + vData[samples - (i + 1)] /= weighingFactor; +#endif + } + } + } +} + +// Private functions + +template +uint_fast8_t ArduinoFFT::exponent(uint_fast16_t value) const { + // Calculates the base 2 logarithm of a value + uint_fast8_t result = 0; + while (value >>= 1) + result++; + return result; +} + +template +void ArduinoFFT::findMaxY(T *vData, uint_fast16_t length, T *maxY, + uint_fast16_t *index) const { + *maxY = 0; + *index = 0; + // If sampling_frequency = 2 * max_frequency in signal, + // value would be stored at position samples/2 + for (uint_fast16_t i = 1; i < length; i++) { + if ((vData[i - 1] < vData[i]) && (vData[i] > vData[i + 1])) { + if (vData[i] > vData[*index]) { + *index = i; + } + } + } + *maxY = vData[*index]; +} + +template +void ArduinoFFT::parabola(T x1, T y1, T x2, T y2, T x3, T y3, T *a, T *b, + T *c) const { + // const T reversed_denom = 1 / ((x1 - x2) * (x1 - x3) * (x2 - x3)); + // This is a special case in which the three X coordinates are three positive, + // consecutive integers. Therefore the reverse denominator will always be -0.5 + const T reversed_denom = -0.5; *a = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) * reversed_denom; *b = (x3 * x3 * (y1 - y2) + x2 * x2 * (y3 - y1) + x1 * x1 * (y2 - y3)) * @@ -533,19 +452,67 @@ void arduinoFFT::Parabola(double x1, double y1, double x2, double y2, double x3, reversed_denom; } -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 >>= 1) - result++; - return (result); +template void ArduinoFFT::swap(T *a, T *b) const { + T temp = *a; + *a = *b; + *b = temp; } -// Private functions - -void arduinoFFT::Swap(double *x, double *y) { - double temp = *x; - *x = *y; - *y = temp; +#ifdef FFT_SQRT_APPROXIMATION +// Fast inverse square root aka "Quake 3 fast inverse square root", multiplied +// by x. Uses one iteration of Halley's method for precision. See: +// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Iterative_methods_for_reciprocal_square_roots +// And: https://github.com/HorstBaerbel/approx +template float ArduinoFFT::sqrt_internal(float x) const { + union // get bits for floating point value + { + float x; + int32_t i; + } u; + u.x = x; + u.i = 0x5f375a86 - (u.i >> 1); // gives initial guess y0. + float xu = x * u.x; + float xu2 = xu * u.x; + // Halley's method, repeating increases accuracy + u.x = (0.125 * 3.0) * xu * (5.0 - xu2 * ((10.0 / 3.0) - xu2)); + return u.x; } + +template double ArduinoFFT::sqrt_internal(double x) const { + // According to HosrtBaerbel, on the ESP32 the approximation is not faster, so + // we use the standard function +#ifdef ESP32 + return sqrt(x); +#else + union // get bits for floating point value + { + double x; + int64_t i; + } u; + u.x = x; + u.i = 0x5fe6ec85e7de30da - (u.i >> 1); // gives initial guess y0. + double xu = x * u.x; + double xu2 = xu * u.x; + // Halley's method, repeating increases accuracy + u.x = (0.125 * 3.0) * xu * (5.0 - xu2 * ((10.0 / 3.0) - xu2)); + return u.x; +#endif +} +#endif + +template +const T ArduinoFFT::_WindowCompensationFactors[10] = { + 1.0000000000 * 2.0, // rectangle (Box car) + 1.8549343278 * 2.0, // hamming + 1.8554726898 * 2.0, // hann + 2.0039186079 * 2.0, // triangle (Bartlett) + 2.8163172034 * 2.0, // nuttall + 2.3673474360 * 2.0, // blackman + 2.7557840395 * 2.0, // blackman nuttall + 2.7929062517 * 2.0, // blackman harris + 3.5659039231 * 2.0, // flat top + 1.5029392863 * 2.0 // welch +}; + +template class ArduinoFFT; +template class ArduinoFFT; diff --git a/src/arduinoFFT.h b/src/arduinoFFT.h index 4878d01..5fb01d0 100644 --- a/src/arduinoFFT.h +++ b/src/arduinoFFT.h @@ -1,21 +1,22 @@ /* - FFT library - Copyright (C) 2010 Didier Longueville - Copyright (C) 2014 Enrique Condes + FFT library + Copyright (C) 2010 Didier Longueville + Copyright (C) 2014 Enrique Condes + Copyright (C) 2020 Bim Overbohm (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 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. + 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 . + You should have received a copy of the GNU General Public License + along with this program. If not, see . */ @@ -38,22 +39,18 @@ #include "defs.h" #include "types.h" #include - +#include #endif -// Define this to use a low-precision square root approximation instead of the +// This definition uses 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 -#ifdef FFT_SQRT_APPROXIMATION -#include -#else +#ifndef FFT_SQRT_APPROXIMATION #define sqrt_internal sqrt #endif -enum class FFTDirection { Reverse, Forward }; +enum class FFTDirection { Forward, Reverse }; enum class FFTWindow { Rectangle, // rectangle (Box car) @@ -65,10 +62,12 @@ enum class FFTWindow { Blackman_Nuttall, // blackman nuttall Blackman_Harris, // blackman harris Flat_top, // flat top - Welch // welch + Welch, // welch + Precompiled // Placeholder for using custom or precompiled window values }; -#define FFT_LIB_REV 0x15 +#define FFT_LIB_REV 0x20 /* Custom constants */ +/* These defines keep compatibility with pre 2.0 code */ #define FFT_FORWARD FFTDirection::Forward #define FFT_REVERSE FFTDirection::Reverse @@ -85,67 +84,95 @@ enum class FFTWindow { FFTWindow::Blackman_Harris /* blackman harris*/ #define FFT_WIN_TYP_FLT_TOP FFTWindow::Flat_top /* flat top */ #define FFT_WIN_TYP_WELCH FFTWindow::Welch /* welch */ -/*Mathematial constants*/ +/* End of compatibility defines */ + +/* Mathematial constants */ #define twoPi 6.28318531 #define fourPi 12.56637061 #define sixPi 18.84955593 -#ifdef __AVR__ -static const double _c1[] PROGMEM = { +template class ArduinoFFT { +public: + ArduinoFFT(); + ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, + bool windowingFactors = false); + + ~ArduinoFFT(); + + void complexToMagnitude(void) const; + void complexToMagnitude(T *vReal, T *vImag, uint_fast16_t samples) const; + + void compute(FFTDirection dir) const; + void compute(T *vReal, T *vImag, uint_fast16_t samples, + FFTDirection dir) const; + void compute(T *vReal, T *vImag, uint_fast16_t samples, uint_fast8_t power, + FFTDirection dir) const; + + void dcRemoval(void) const; + void dcRemoval(T *vData, uint_fast16_t samples) const; + + T majorPeak(void) const; + void majorPeak(T *f, T *v) const; + T majorPeak(T *vData, uint_fast16_t samples, T samplingFrequency) const; + void majorPeak(T *vData, uint_fast16_t samples, T samplingFrequency, + T *frequency, T *magnitude) const; + + T majorPeakParabola(void) const; + void majorPeakParabola(T *frequency, T *magnitude) const; + T majorPeakParabola(T *vData, uint_fast16_t samples, + T samplingFrequency) const; + void majorPeakParabola(T *vData, uint_fast16_t samples, T samplingFrequency, + T *frequency, T *magnitude) const; + + uint8_t revision(void); + + void setArrays(T *vReal, T *vImag, uint_fast16_t samples = 0); + + void windowing(FFTWindow windowType, FFTDirection dir, + bool withCompensation = false); + void windowing(T *vData, uint_fast16_t samples, FFTWindow windowType, + FFTDirection dir, T *windowingFactors = nullptr, + bool withCompensation = false); + +private: + /* Variables */ + static const T _WindowCompensationFactors[10]; +#ifdef FFT_SPEED_OVER_PRECISION + T _oneOverSamples = 0.0; +#endif + bool _isPrecompiled = false; + bool _precompiledWithCompensation = false; + uint_fast8_t _power = 0; + T *_precompiledWindowingFactors; + uint_fast16_t _samples; + T _samplingFrequency; + T *_vImag; + T *_vReal; + FFTWindow _windowFunction; + /* Functions */ + uint_fast8_t exponent(uint_fast16_t value) const; + void findMaxY(T *vData, uint_fast16_t length, T *maxY, + uint_fast16_t *index) const; + void parabola(T x1, T y1, T x2, T y2, T x3, T y3, T *a, T *b, T *c) const; + void swap(T *a, T *b) const; + +#ifdef FFT_SQRT_APPROXIMATION + float sqrt_internal(float x) const; + double sqrt_internal(double x) const; +#endif +}; + +#if defined(__AVR__) && defined(USE_AVR_PROGMEM) +static const float _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 = { +static const float _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 -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); - - void ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples); - void Compute(double *vReal, double *vImag, uint16_t samples, - FFTDirection dir); - void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, - FFTDirection 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, FFTWindow windowType, - FFTDirection dir); - - void ComplexToMagnitude(); - void Compute(FFTDirection dir); - void DCRemoval(); - double MajorPeak(); - void MajorPeak(double *f, double *v); - void Windowing(FFTWindow windowType, FFTDirection dir); - - double MajorPeakParabola(); - -private: - /* Variables */ - uint16_t _samples; - double _samplingFrequency; - double *_vReal; - double *_vImag; - uint8_t _power; - /* Functions */ - void Swap(double *x, double *y); - void Parabola(double x1, double y1, double x2, double y2, double x3, - double y3, double *a, double *b, double *c); -}; #endif