From a9f64fb886709706b07605a5154da6af8b2ebe8f Mon Sep 17 00:00:00 2001 From: Enrique Condes Date: Wed, 6 Mar 2024 13:56:17 +0800 Subject: [PATCH] Version 2.0 --- Examples/FFT_01/FFT_01.ino | 8 +- Examples/FFT_02/FFT_02.ino | 12 +- Examples/FFT_03/FFT_03.ino | 2 +- Examples/FFT_04/FFT_04.ino | 8 +- Examples/FFT_05/FFT_05.ino | 8 +- Examples/FFT_speedup/FFT_speedup.ino | 12 +- README.md | 113 +----- changeLog.txt | 40 -- keywords.txt | 19 +- library.json | 2 +- library.properties | 2 +- src/arduinoFFT.cpp | 518 ++++++++++++++++++++++++ src/arduinoFFT.h | 572 ++++++--------------------- 13 files changed, 688 insertions(+), 628 deletions(-) delete mode 100644 changeLog.txt create mode 100644 src/arduinoFFT.cpp diff --git a/Examples/FFT_01/FFT_01.ino b/Examples/FFT_01/FFT_01.ino index 22b5024..e16368a 100644 --- a/Examples/FFT_01/FFT_01.ino +++ b/Examples/FFT_01/FFT_01.ino @@ -3,7 +3,7 @@ Example of use of the FFT libray Copyright (C) 2014 Enrique Condes - Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) + 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 @@ -64,11 +64,11 @@ void setup() void loop() { /* Build raw data */ - double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read + double ratio = twoPi * signalFrequency / samplingFrequency; // Fraction of a complete cycle stored at each sample (in radians) for (uint16_t i = 0; i < samples; i++) { - 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*/ + vReal[i] = int8_t(amplitude * sin(i * ratio) / 2.0);/* Build data with positive and negative values*/ + //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 } /* Print the results of the simulated sampling according to time */ diff --git a/Examples/FFT_02/FFT_02.ino b/Examples/FFT_02/FFT_02.ino index 7164dab..a8cbc63 100644 --- a/Examples/FFT_02/FFT_02.ino +++ b/Examples/FFT_02/FFT_02.ino @@ -1,12 +1,12 @@ /* 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. + 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 (header-only, template, speed improvements) + 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 @@ -65,10 +65,10 @@ void loop() for(double frequency = startFrequency; frequency<=stopFrequency; frequency+=step_size) { /* Build raw data */ - double cycles = (((samples-1) * frequency) / sampling); + double ratio = twoPi * frequency / sampling; // Fraction of a complete cycle stored at each sample (in radians) for (uint16_t i = 0; i < samples; i++) { - vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0); + vReal[i] = int8_t(amplitude * sin(i * ratio) / 2.0);/* Build data with positive and negative values*/ vImag[i] = 0; //Reset the imaginary values vector for each new frequency } /*Serial.println("Data:"); diff --git a/Examples/FFT_03/FFT_03.ino b/Examples/FFT_03/FFT_03.ino index 2e50613..67ed135 100644 --- a/Examples/FFT_03/FFT_03.ino +++ b/Examples/FFT_03/FFT_03.ino @@ -3,7 +3,7 @@ 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 (header-only, template, speed improvements) + 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 diff --git a/Examples/FFT_04/FFT_04.ino b/Examples/FFT_04/FFT_04.ino index b125991..3934c10 100644 --- a/Examples/FFT_04/FFT_04.ino +++ b/Examples/FFT_04/FFT_04.ino @@ -3,7 +3,7 @@ Example of use of the FFT libray Copyright (C) 2018 Enrique Condes - Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) + 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 @@ -63,11 +63,11 @@ void setup() void loop() { /* Build raw data */ - double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read + double ratio = twoPi * signalFrequency / samplingFrequency; // Fraction of a complete cycle stored at each sample (in radians) for (uint16_t i = 0; i < samples; i++) { - 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*/ + vReal[i] = int8_t(amplitude * sin(i * ratio) / 2.0);/* Build data with positive and negative values*/ + //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.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */ diff --git a/Examples/FFT_05/FFT_05.ino b/Examples/FFT_05/FFT_05.ino index a6f4df7..4354fcd 100644 --- a/Examples/FFT_05/FFT_05.ino +++ b/Examples/FFT_05/FFT_05.ino @@ -3,7 +3,7 @@ Example of use of the FFT libray Copyright (C) 2014 Enrique Condes - Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) + 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 @@ -65,11 +65,11 @@ void setup() void loop() { /* Build raw data */ - double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read + double ratio = twoPi * signalFrequency / samplingFrequency; // Fraction of a complete cycle stored at each sample (in radians) for (uint16_t i = 0; i < samples; i++) { - 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*/ + vReal[i] = int8_t(amplitude * sin(i * ratio) / 2.0);/* Build data with positive and negative values*/ + //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 } /* Print the results of the simulated sampling according to time */ diff --git a/Examples/FFT_speedup/FFT_speedup.ino b/Examples/FFT_speedup/FFT_speedup.ino index a059a17..e91a2fe 100644 --- a/Examples/FFT_speedup/FFT_speedup.ino +++ b/Examples/FFT_speedup/FFT_speedup.ino @@ -1,9 +1,9 @@ /* 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 + with speedup through different arduinoFFT options. Based on examples/FFT_03/FFT_03.ino - Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) + 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 @@ -47,14 +47,8 @@ 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); +ArduinoFFT FFT = ArduinoFFT(vReal, vImag, samples, samplingFrequency, true); #define SCL_INDEX 0x00 #define SCL_TIME 0x01 diff --git a/README.md b/README.md index f9229ef..dbafd42 100644 --- a/README.md +++ b/README.md @@ -4,8 +4,10 @@ arduinoFFT # 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). See here [how to migrate from 1.x to 2.x](#migrating-from-1x-to-2x). -Tested on Arduino 1.6.11 and 1.8.10. + +This is version 2.0 of the library, which has a different [API](#api). + +Tested on Arduino 1.8.19 and 2.3.2. ## Installation on Arduino @@ -15,17 +17,17 @@ Use the Arduino Library Manager to install and keep it updated. Just look for ar To install this library, just place this entire folder as a subfolder in your Arduino installation. 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` (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) ## 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. +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: @@ -33,97 +35,4 @@ select arduinoFTT. This will add a corresponding line to the top of your sketch ## API -* ```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. You can deallocate `vReal` and `vImag` after you are done using the library, or only use specific library functions that only need one of those arrays. - -```C++ -const uint32_t nrOfSamples = 1024; -auto real = new float[nrOfSamples]; -auto imag = new float[nrOfSamples]; -auto fft = ArduinoFFT(real, imag, nrOfSamples, 10000); -// ... fill real + imag and use it ... -fft.compute(); -fft.complexToMagnitude(); -delete [] imag; -// ... continue using real and only functions that use real ... -auto peak = fft.majorPeak(); -``` -* ```~ArduinoFFT()``` -Destructor. -* ```void complexToMagnitude() const;``` -Convert complex values to their magnitude and store in vReal. Uses vReal and vImag. -* ```void compute(FFTDirection dir) const;``` -Calcuates the Fast Fourier Transform. Uses vReal and vImag. -* ```void dcRemoval() const;``` -Removes the DC component from the sample data. Uses vReal. -* ```T majorPeak() const;``` -Returns the frequency of the biggest spike in the analyzed signal. Uses vReal. -* ```void majorPeak(T &frequency, T &value) const;``` -Returns the frequency and the value of the biggest spike in the analyzed signal. Uses vReal. -* ```uint8_t revision() const;``` -Returns the library revision. -* ```void setArrays(T *vReal, T *vImag);``` -Replace the data array pointers. -* ```void windowing(FFTWindow windowType, FFTDirection dir, bool withCompensation = false);``` -Performs a windowing function on the values array. Uses vReal. The possible windowing options are: - * Rectangle - * Hamming - * Hann - * Triangle - * Nuttall - * Blackman - * Blackman_Nuttall - * Blackman_Harris - * Flat_top - * Welch - - If `withCompensation` == true, the following compensation factors are used: - * Rectangle: 1.0 * 2.0 - * Hamming: 1.8549343278 * 2.0 - * Hann: 1.8554726898 * 2.0 - * Triangle: 2.0039186079 * 2.0 - * Nuttall: 2.8163172034 * 2.0 - * Blackman: 2.3673474360 * 2.0 - * Blackman Nuttall: 2.7557840395 * 2.0 - * Blackman Harris: 2.7929062517 * 2.0 - * Flat top: 3.5659039231 * 2.0 - * Welch: 1.5029392863 * 2.0 - -## 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`. - -See the `FFT_speedup.ino` example in `Examples/FFT_speedup/FFT_speedup.ino`. - -# Migrating from 1.x to 2.x - -* The function signatures where you could pass in pointers were deprecated and have been removed. Pass in pointers to your real / imaginary array in the ArduinoFFT() constructor. If you have the need to replace those pointers during usage of the library (e.g. to free memory) you can do the following: - -```C++ -const uint32_t nrOfSamples = 1024; -auto real = new float[nrOfSamples]; -auto imag = new float[nrOfSamples]; -auto fft = ArduinoFFT(real, imag, nrOfSamples, 10000); -// ... fill real + imag and use it ... -fft.compute(); -fft.complexToMagnitude(); -delete [] real; -// ... replace vReal in library with imag ... -fft.setArrays(imag, nullptr); -// ... keep doing whatever ... -``` -* All function names are camelCase case now (start with lower-case character), e.g. "windowing()" instead of "Windowing()". - -## TODO -* Ratio table for windowing function. -* Document windowing functions advantages and disadvantages. -* Optimize usage and arguments. -* Add new windowing functions. -* ~~Spectrum table?~~ +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 d49b854..0000000 --- a/changeLog.txt +++ /dev/null @@ -1,40 +0,0 @@ -02/22/20 v1.9.2 -Fix compilation on AVR systems. - -02/22/20 v1.9.1 -Add setArrays() function because of issue #32. -Add API migration info to README and improve README. -Use better sqrtf() approximation. - -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 (~40% 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%. - -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 3807cdb..49fd943 100644 --- a/keywords.txt +++ b/keywords.txt @@ -17,11 +17,11 @@ FFTWindow KEYWORD1 complexToMagnitude KEYWORD2 compute KEYWORD2 dcRemoval KEYWORD2 -windowing KEYWORD2 -exponent KEYWORD2 -revision KEYWORD2 majorPeak KEYWORD2 +majorPeakParabola KEYWORD2 +revision KEYWORD2 setArrays KEYWORD2 +windowing KEYWORD2 ####################################### # Constants (LITERAL1) @@ -29,13 +29,14 @@ setArrays KEYWORD2 Forward LITERAL1 Reverse LITERAL1 -Rectangle LITERAL1 + +Blackman LITERAL1 +Blackman_Harris LITERAL1 +Blackman_Nuttall LITERAL1 +Flat_top LITERAL1 Hamming LITERAL1 Hann LITERAL1 -Triangle LITERAL1 Nuttall LITERAL1 -Blackman LITERAL1 -Blackman_Nuttall LITERAL1 -Blackman_Harris LITERAL1 -Flat_top LITERAL1 +Rectangle LITERAL1 +Triangle LITERAL1 Welch LITERAL1 diff --git a/library.json b/library.json index 6c35419..de87c29 100644 --- a/library.json +++ b/library.json @@ -25,7 +25,7 @@ "email": "bim.overbohm@googlemail.com" } ], - "version": "1.9.2", + "version": "2.0", "frameworks": ["arduino","mbed","espidf"], "platforms": "*" } diff --git a/library.properties b/library.properties index 0a90947..5826255 100644 --- a/library.properties +++ b/library.properties @@ -1,5 +1,5 @@ name=arduinoFFT -version=1.9.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/arduinoFFT.cpp b/src/arduinoFFT.cpp new file mode 100644 index 0000000..1bb5b22 --- /dev/null +++ b/src/arduinoFFT.cpp @@ -0,0 +1,518 @@ +/* + 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 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" + +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 +} + +template ArduinoFFT::~ArduinoFFT(void) { + // Destructor + if (_precompiledWindowingFactors) { + delete [] _precompiledWindowingFactors; + } +} + +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 +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 == FFTDirection::Reverse) + swap(&vImag[i], &vImag[j]); + } + uint_fast16_t k = (samples >> 1); + + while (k <= j) { + j -= k; + k >>= 1; + } + j += k; + } + // 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; + T u1 = 1.0; + T u2 = 0.0; + for (j = 0; j < l1; j++) { + 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; + } + T z = ((u1 * c1) - (u2 * c2)); + u2 = ((u1 * c2) + (u2 * c1)); + u1 = z; + } + +#if defined(__AVR__) && defined(USE_AVR_PROGMEM) + c2 = pgm_read_float_near(&(_c2[l])); + c1 = pgm_read_float_near(&(_c1[l])); +#else + T cTemp = 0.5 * c1; + c2 = sqrt_internal(0.5 - cTemp); + c1 = sqrt_internal(0.5 + cTemp); +#endif + + if (dir == FFTDirection::Forward) { + c2 = -c2; + } + } + // Scaling for reverse transform + 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 + } + } +} + +template void ArduinoFFT::dcRemoval(void) const { + dcRemoval(this->_vReal, this->_samples); +} + +template +void ArduinoFFT::dcRemoval(T *vData, uint_fast16_t samples) const { + // calculate the mean of vData + T mean = 0; + for (uint_fast16_t i = 0; i < samples; i++) { + mean += vData[i]; + } + mean /= samples; + // Subtract the mean from vData + for (uint_fast16_t i = 0; i < samples; i++) { + vData[i] -= mean; + } +} + +template T ArduinoFFT::majorPeak(void) const { + return majorPeak(this->_vReal, this->_samples, this->_samplingFrequency); +} + +template void ArduinoFFT::majorPeak(T *f, T *v) const { + majorPeak(this->_vReal, this->_samples, this->_samplingFrequency, f, v); +} + +template +T ArduinoFFT::majorPeak(T *vData, uint_fast16_t samples, + T samplingFrequency) const { + T frequency; + majorPeak(vData, samples, samplingFrequency, &frequency, nullptr); + return frequency; +} + +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); + + 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 + *frequency = interpolatedX; + if (magnitude != nullptr) { +#if defined(ESP8266) || defined(ESP32) + *magnitude = fabs(vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) + + vData[IndexOfMaxY + 1]); +#else + *magnitude = abs(vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) + + vData[IndexOfMaxY + 1]); +#endif + } +} + +template T ArduinoFFT::majorPeakParabola(void) const { + T freq = 0; + majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency, + &freq, nullptr); + return freq; +} + +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)) * + reversed_denom; + *c = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + + x1 * x2 * (x1 - x2) * y3) * + reversed_denom; +} + +template void ArduinoFFT::swap(T *a, T *b) const { + T temp = *a; + *a = *b; + *b = 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 cb968ff..5fb01d0 100644 --- a/src/arduinoFFT.h +++ b/src/arduinoFFT.h @@ -1,22 +1,22 @@ /* - FFT library - Copyright (C) 2010 Didier Longueville - Copyright (C) 2014 Enrique Condes - Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) + 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 . */ @@ -29,472 +29,150 @@ #include "WProgram.h" /* This is where the standard Arduino code lies */ #endif #else -#include #include +#include + #ifdef __AVR__ #include #include #endif -#include #include "defs.h" #include "types.h" +#include +#include #endif -// Define this to use reciprocal multiplication for division and some more speedups that might decrease precision -//#define FFT_SPEED_OVER_PRECISION +// 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. -// 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 - -#ifdef FFT_SQRT_APPROXIMATION - #include -#else - #ifndef sqrt_internal - #define sqrt_internal sqrt - #endif +#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) - 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 +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 + Precompiled // Placeholder for using custom or precompiled window values }; +#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 -template -class ArduinoFFT -{ +/* Windowing type */ +#define FFT_WIN_TYP_RECTANGLE FFTWindow::Rectangle /* rectangle (Box car) */ +#define FFT_WIN_TYP_HAMMING FFTWindow::Hamming /* hamming */ +#define FFT_WIN_TYP_HANN FFTWindow::Hann /* hann */ +#define FFT_WIN_TYP_TRIANGLE FFTWindow::Triangle /* triangle (Bartlett) */ +#define FFT_WIN_TYP_NUTTALL FFTWindow::Nuttall /* nuttall */ +#define FFT_WIN_TYP_BLACKMAN FFTWindow::Blackman /* blackman */ +#define FFT_WIN_TYP_BLACKMAN_NUTTALL \ + FFTWindow::Blackman_Nuttall /* blackman nuttall */ +#define FFT_WIN_TYP_BLACKMAN_HARRIS \ + FFTWindow::Blackman_Harris /* blackman harris*/ +#define FFT_WIN_TYP_FLT_TOP FFTWindow::Flat_top /* flat top */ +#define FFT_WIN_TYP_WELCH FFTWindow::Welch /* welch */ +/* End of compatibility defines */ + +/* Mathematial constants */ +#define twoPi 6.28318531 +#define fourPi 12.56637061 +#define sixPi 18.84955593 + +template class ArduinoFFT { public: - // 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++; - } - } + ArduinoFFT(); + ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, + bool windowingFactors = false); - // Destructor - ~ArduinoFFT() - { - } + ~ArduinoFFT(); - // Get library revision - static uint8_t revision() - { - return 0x19; - } + void complexToMagnitude(void) const; + void complexToMagnitude(T *vReal, T *vImag, uint_fast16_t samples) const; - // Replace the data array pointers - void setArrays(T *vReal, T *vImag) - { - _vReal = vReal; - _vImag = vImag; - } + 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; - // 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__ - uint_fast8_t 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_float_near(&(_c2[index])); - c1 = pgm_read_float_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 dcRemoval(void) const; + void dcRemoval(T *vData, uint_fast16_t samples) const; - 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])); - } - } + 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; - 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; - } - } + 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; - void windowing(FFTWindow windowType, FFTDirection dir, bool withCompensation = false) - { - // check if values are already pre-computed for the correct window type and compensation - if (_windowWeighingFactors && _weighingFactorsComputed && - _weighingFactorsFFTWindow == windowType && - _weighingFactorsWithCompensation == withCompensation) - { - // 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); - T compensationFactor = _WindowCompensationFactors[static_cast(windowType)]; - 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 (withCompensation) - { - weighingFactor *= compensationFactor; - } - 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; - _weighingFactorsWithCompensation = withCompensation; - _weighingFactorsComputed = true; - } - } + uint8_t revision(void); - 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 setArrays(T *vReal, T *vImag, uint_fast16_t samples = 0); - void majorPeak(T &frequency, T &value) 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 - frequency = interpolatedX; - value = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]); - } + 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: -#ifdef __AVR__ - static const float _c1[] PROGMEM; - static const float _c2[] PROGMEM; + /* Variables */ + static const T _WindowCompensationFactors[10]; +#ifdef FFT_SPEED_OVER_PRECISION + T _oneOverSamples = 0.0; #endif - static const T _WindowCompensationFactors[10]; - - // 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; - } + 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 - // 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 - static inline V sqrt_internal(typename std::enable_if::value, V>::type x) - { - union // get bits for float 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; - u.x = (0.125 * 3.0) * xu * (5.0 - xu2 * ((10.0 / 3.0) - xu2)); // Halley's method, repeating increases accuracy - return u.x; - } - - template - static inline V sqrt_internal(typename std::enable_if::value, V>::type x) - { - // 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 float 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; - u.x = (0.125 * 3.0) * xu * (5.0 - xu2 * ((10.0 / 3.0) - xu2)); // Halley's method, repeating increases accuracy - return u.x; - #endif - } + float sqrt_internal(float x) const; + double sqrt_internal(double x) const; #endif - - /* Variables */ - T *_vReal = nullptr; - T *_vImag = nullptr; - uint_fast16_t _samples = 0; -#ifdef FFT_SPEED_OVER_PRECISION - T _oneOverSamples = 0.0; -#endif - T _samplingFrequency = 0; - T *_windowWeighingFactors = nullptr; - FFTWindow _weighingFactorsFFTWindow; - bool _weighingFactorsWithCompensation = false; - bool _weighingFactorsComputed = false; - uint_fast8_t _power = 0; }; -#ifdef __AVR__ -template -const float ArduinoFFT::_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}; - -template -const float ArduinoFFT::_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}; +#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 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 -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 -}; - #endif