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

Version 2.0

pull/91/head
Enrique Condes 2024-03-06 13:56:17 +08:00
rodzic 419d7b044e
commit a9f64fb886
13 zmienionych plików z 688 dodań i 628 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

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

@ -3,7 +3,7 @@
Example of use of the FFT libray Example of use of the FFT libray
Copyright (C) 2014 Enrique Condes Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) Copyright (C) 2020 Bim Overbohm (template, speed improvements)
This program is free software: you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
@ -64,11 +64,11 @@ void setup()
void loop() void loop()
{ {
/* Build raw data */ /* 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++) 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] = int8_t(amplitude * sin(i * ratio) / 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] = 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 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 */ /* Print the results of the simulated sampling according to time */

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

@ -1,12 +1,12 @@
/* /*
Example of use of the FFT libray to compute FFT for several signals over a range of frequencies. 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. 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. 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 sketch shows the time that the computing is taking.
Copyright (C) 2014 Enrique Condes Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) Copyright (C) 2020 Bim Overbohm (template, speed improvements)
This program is free software: you can redistribute it and/or modify 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 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) for(double frequency = startFrequency; frequency<=stopFrequency; frequency+=step_size)
{ {
/* Build raw data */ /* 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++) 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 vImag[i] = 0; //Reset the imaginary values vector for each new frequency
} }
/*Serial.println("Data:"); /*Serial.println("Data:");

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

@ -3,7 +3,7 @@
Example of use of the FFT libray to compute FFT for a signal sampled through the ADC. Example of use of the FFT libray to compute FFT for a signal sampled through the ADC.
Copyright (C) 2018 Enrique Condés and Ragnar Ranøyen Homb Copyright (C) 2018 Enrique Condés and Ragnar Ranøyen Homb
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) Copyright (C) 2020 Bim Overbohm (template, speed improvements)
This program is free software: you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by

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

@ -3,7 +3,7 @@
Example of use of the FFT libray Example of use of the FFT libray
Copyright (C) 2018 Enrique Condes Copyright (C) 2018 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) Copyright (C) 2020 Bim Overbohm (template, speed improvements)
This program is free software: you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
@ -63,11 +63,11 @@ void setup()
void loop() void loop()
{ {
/* Build raw data */ /* 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++) 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] = int8_t(amplitude * sin(i * ratio) / 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] = 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 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 */ FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */

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

@ -3,7 +3,7 @@
Example of use of the FFT libray Example of use of the FFT libray
Copyright (C) 2014 Enrique Condes Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (header-only, template, speed improvements) Copyright (C) 2020 Bim Overbohm (template, speed improvements)
This program is free software: you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
@ -65,11 +65,11 @@ void setup()
void loop() void loop()
{ {
/* Build raw data */ /* 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++) 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] = int8_t(amplitude * sin(i * ratio) / 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] = 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 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 */ /* Print the results of the simulated sampling according to time */

12
Examples/FFT_speedup/FFT_speedup.ino
Wyświetl plik

@ -1,9 +1,9 @@
/* /*
Example of use of the FFT libray to compute FFT for a signal sampled through the ADC 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 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 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 vReal[samples];
float vImag[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 */ /* Create FFT object with weighing factor storage */
ArduinoFFT<float> FFT = ArduinoFFT<float>(vReal, vImag, samples, samplingFrequency, weighingFactors); ArduinoFFT<float> FFT = ArduinoFFT<float>(vReal, vImag, samples, samplingFrequency, true);
#define SCL_INDEX 0x00 #define SCL_INDEX 0x00
#define SCL_TIME 0x01 #define SCL_TIME 0x01

113
README.md
Wyświetl plik

@ -4,8 +4,10 @@ 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 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 ## 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: 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` (this library's folder)
`Arduino\libraries\arduinoFTT\src\arduinoFTT.h` (the library header file. include this in your project) `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\keywords.txt` (the syntax coloring file)
`Arduino\libraries\arduinoFTT\Examples` (the examples in the "open" menu) `Arduino\libraries\arduinoFTT\Examples` (the examples in the "open" menu)
`Arduino\libraries\arduinoFTT\LICENSE` (GPL license file) `Arduino\libraries\arduinoFTT\LICENSE` (GPL license file)
`Arduino\libraries\arduinoFTT\README.md` (this file) `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. 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 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: 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 ## API
* ```ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, T * weighingFactors = nullptr);``` Documentation was moved to the project's [wiki](https://github.com/kosme/arduinoFFT/wiki).
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<float>(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<float>(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?~~

40
changeLog.txt
Wyświetl plik

@ -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.

19
keywords.txt
Wyświetl plik

@ -17,11 +17,11 @@ FFTWindow KEYWORD1
complexToMagnitude KEYWORD2 complexToMagnitude KEYWORD2
compute KEYWORD2 compute KEYWORD2
dcRemoval KEYWORD2 dcRemoval KEYWORD2
windowing KEYWORD2
exponent KEYWORD2
revision KEYWORD2
majorPeak KEYWORD2 majorPeak KEYWORD2
majorPeakParabola KEYWORD2
revision KEYWORD2
setArrays KEYWORD2 setArrays KEYWORD2
windowing KEYWORD2
####################################### #######################################
# Constants (LITERAL1) # Constants (LITERAL1)
@ -29,13 +29,14 @@ setArrays KEYWORD2
Forward LITERAL1 Forward LITERAL1
Reverse LITERAL1 Reverse LITERAL1
Rectangle LITERAL1
Blackman LITERAL1
Blackman_Harris LITERAL1
Blackman_Nuttall LITERAL1
Flat_top LITERAL1
Hamming LITERAL1 Hamming LITERAL1
Hann LITERAL1 Hann LITERAL1
Triangle LITERAL1
Nuttall LITERAL1 Nuttall LITERAL1
Blackman LITERAL1 Rectangle LITERAL1
Blackman_Nuttall LITERAL1 Triangle LITERAL1
Blackman_Harris LITERAL1
Flat_top LITERAL1
Welch LITERAL1 Welch LITERAL1

2
library.json
Wyświetl plik

@ -25,7 +25,7 @@
"email": "bim.overbohm@googlemail.com" "email": "bim.overbohm@googlemail.com"
} }
], ],
"version": "1.9.2", "version": "2.0",
"frameworks": ["arduino","mbed","espidf"], "frameworks": ["arduino","mbed","espidf"],
"platforms": "*" "platforms": "*"
} }

2
library.properties
Wyświetl plik

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

518
src/arduinoFFT.cpp 100644
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "arduinoFFT.h"
template <typename T> ArduinoFFT<T>::ArduinoFFT() {}
template <typename T>
ArduinoFFT<T>::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 <typename T> ArduinoFFT<T>::~ArduinoFFT(void) {
// Destructor
if (_precompiledWindowingFactors) {
delete [] _precompiledWindowingFactors;
}
}
template <typename T> void ArduinoFFT<T>::complexToMagnitude(void) const {
complexToMagnitude(this->_vReal, this->_vImag, this->_samples);
}
template <typename T>
void ArduinoFFT<T>::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 <typename T> void ArduinoFFT<T>::compute(FFTDirection dir) const {
compute(this->_vReal, this->_vImag, this->_samples, exponent(this->_samples),
dir);
}
template <typename T>
void ArduinoFFT<T>::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 <typename T>
void ArduinoFFT<T>::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 <typename T> void ArduinoFFT<T>::dcRemoval(void) const {
dcRemoval(this->_vReal, this->_samples);
}
template <typename T>
void ArduinoFFT<T>::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 <typename T> T ArduinoFFT<T>::majorPeak(void) const {
return majorPeak(this->_vReal, this->_samples, this->_samplingFrequency);
}
template <typename T> void ArduinoFFT<T>::majorPeak(T *f, T *v) const {
majorPeak(this->_vReal, this->_samples, this->_samplingFrequency, f, v);
}
template <typename T>
T ArduinoFFT<T>::majorPeak(T *vData, uint_fast16_t samples,
T samplingFrequency) const {
T frequency;
majorPeak(vData, samples, samplingFrequency, &frequency, nullptr);
return frequency;
}
template <typename T>
void ArduinoFFT<T>::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 <typename T> T ArduinoFFT<T>::majorPeakParabola(void) const {
T freq = 0;
majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency,
&freq, nullptr);
return freq;
}
template <typename T>
void ArduinoFFT<T>::majorPeakParabola(T *frequency, T *magnitude) const {
majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency,
frequency, magnitude);
}
template <typename T>
T ArduinoFFT<T>::majorPeakParabola(T *vData, uint_fast16_t samples,
T samplingFrequency) const {
T freq = 0;
majorPeakParabola(vData, samples, samplingFrequency, &freq, nullptr);
return freq;
}
template <typename T>
void ArduinoFFT<T>::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 <typename T> uint8_t ArduinoFFT<T>::revision(void) {
return (FFT_LIB_REV);
}
// Replace the data array pointers
template <typename T>
void ArduinoFFT<T>::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 <typename T>
void ArduinoFFT<T>::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 <typename T>
void ArduinoFFT<T>::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<uint_fast8_t>(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 <typename T>
uint_fast8_t ArduinoFFT<T>::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 <typename T>
void ArduinoFFT<T>::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 <typename T>
void ArduinoFFT<T>::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 <typename T> void ArduinoFFT<T>::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 <typename T> float ArduinoFFT<T>::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 <typename T> double ArduinoFFT<T>::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 <typename T>
const T ArduinoFFT<T>::_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<double>;
template class ArduinoFFT<float>;

572
src/arduinoFFT.h
Wyświetl plik

@ -1,22 +1,22 @@
/* /*
FFT library FFT library
Copyright (C) 2010 Didier Longueville Copyright (C) 2010 Didier Longueville
Copyright (C) 2014 Enrique Condes 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 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or the Free Software Foundation, either version 3 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. along with this program. If not, see <http://www.gnu.org/licenses/>.
*/ */
@ -29,472 +29,150 @@
#include "WProgram.h" /* This is where the standard Arduino code lies */ #include "WProgram.h" /* This is where the standard Arduino code lies */
#endif #endif
#else #else
#include <stdlib.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h>
#ifdef __AVR__ #ifdef __AVR__
#include <avr/io.h> #include <avr/io.h>
#include <avr/pgmspace.h> #include <avr/pgmspace.h>
#endif #endif
#include <math.h>
#include "defs.h" #include "defs.h"
#include "types.h" #include "types.h"
#include <math.h>
#include <stdint.h>
#endif #endif
// Define this to use reciprocal multiplication for division and some more speedups that might decrease precision // This definition uses a low-precision square root approximation instead of the
//#define FFT_SPEED_OVER_PRECISION // 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 #ifndef FFT_SQRT_APPROXIMATION
// This might only work for specific use cases, but is significantly faster. Only works for ArduinoFFT<float>. #define sqrt_internal sqrt
//#define FFT_SQRT_APPROXIMATION
#ifdef FFT_SQRT_APPROXIMATION
#include <type_traits>
#else
#ifndef sqrt_internal
#define sqrt_internal sqrt
#endif
#endif #endif
enum class FFTDirection enum class FFTDirection { Forward, Reverse };
{
Reverse,
Forward
};
enum class FFTWindow enum class FFTWindow {
{ Rectangle, // rectangle (Box car)
Rectangle, // rectangle (Box car) Hamming, // hamming
Hamming, // hamming Hann, // hann
Hann, // hann Triangle, // triangle (Bartlett)
Triangle, // triangle (Bartlett) Nuttall, // nuttall
Nuttall, // nuttall Blackman, // blackman
Blackman, //blackman Blackman_Nuttall, // blackman nuttall
Blackman_Nuttall, // blackman nuttall Blackman_Harris, // blackman harris
Blackman_Harris, // blackman harris Flat_top, // flat top
Flat_top, // flat top Welch, // welch
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 <typename T> /* Windowing type */
class ArduinoFFT #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 <typename T> class ArduinoFFT {
public: public:
// Constructor ArduinoFFT();
ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, T *windowWeighingFactors = nullptr) ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency,
: _vReal(vReal) bool windowingFactors = false);
, _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++;
}
}
// Destructor ~ArduinoFFT();
~ArduinoFFT()
{
}
// Get library revision void complexToMagnitude(void) const;
static uint8_t revision() void complexToMagnitude(T *vReal, T *vImag, uint_fast16_t samples) const;
{
return 0x19;
}
// Replace the data array pointers void compute(FFTDirection dir) const;
void setArrays(T *vReal, T *vImag) void compute(T *vReal, T *vImag, uint_fast16_t samples,
{ FFTDirection dir) const;
_vReal = vReal; void compute(T *vReal, T *vImag, uint_fast16_t samples, uint_fast8_t power,
_vImag = vImag; FFTDirection dir) const;
}
// Computes in-place complex-to-complex FFT void dcRemoval(void) const;
void compute(FFTDirection dir) const void dcRemoval(T *vData, uint_fast16_t samples) 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 complexToMagnitude() const T majorPeak(void) const;
{ void majorPeak(T *f, T *v) const;
// vM is half the size of vReal and vImag T majorPeak(T *vData, uint_fast16_t samples, T samplingFrequency) const;
for (uint_fast16_t i = 0; i < this->_samples; i++) void majorPeak(T *vData, uint_fast16_t samples, T samplingFrequency,
{ T *frequency, T *magnitude) const;
this->_vReal[i] = sqrt_internal(sq(this->_vReal[i]) + sq(this->_vImag[i]));
}
}
void dcRemoval() const T majorPeakParabola(void) const;
{ void majorPeakParabola(T *frequency, T *magnitude) const;
// calculate the mean of vData T majorPeakParabola(T *vData, uint_fast16_t samples,
T mean = 0; T samplingFrequency) const;
for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++) void majorPeakParabola(T *vData, uint_fast16_t samples, T samplingFrequency,
{ T *frequency, T *magnitude) const;
mean += this->_vReal[i];
}
mean /= this->_samples;
// Subtract the mean from vData
for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
this->_vReal[i] -= mean;
}
}
void windowing(FFTWindow windowType, FFTDirection dir, bool withCompensation = false) uint8_t revision(void);
{
// 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<uint_fast8_t>(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;
}
}
T majorPeak() const void setArrays(T *vReal, T *vImag, uint_fast16_t samples = 0);
{
T maxY = 0;
uint_fast16_t IndexOfMaxY = 0;
//If sampling_frequency = 2 * max_frequency in signal,
//value would be stored at position samples/2
for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
{
if ((this->_vReal[i - 1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i + 1]))
{
if (this->_vReal[i] > maxY)
{
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
T delta = 0.5 * ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) / (this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]));
T interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1);
if (IndexOfMaxY == (this->_samples >> 1))
{
//To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples);
}
// returned value: interpolated frequency peak apex
return interpolatedX;
}
void majorPeak(T &frequency, T &value) const void windowing(FFTWindow windowType, FFTDirection dir,
{ bool withCompensation = false);
T maxY = 0; void windowing(T *vData, uint_fast16_t samples, FFTWindow windowType,
uint_fast16_t IndexOfMaxY = 0; FFTDirection dir, T *windowingFactors = nullptr,
//If sampling_frequency = 2 * max_frequency in signal, bool withCompensation = false);
//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]);
}
private: private:
#ifdef __AVR__ /* Variables */
static const float _c1[] PROGMEM; static const T _WindowCompensationFactors[10];
static const float _c2[] PROGMEM; #ifdef FFT_SPEED_OVER_PRECISION
T _oneOverSamples = 0.0;
#endif #endif
static const T _WindowCompensationFactors[10]; bool _isPrecompiled = false;
bool _precompiledWithCompensation = false;
// Mathematial constants uint_fast8_t _power = 0;
#ifndef TWO_PI T *_precompiledWindowingFactors;
static constexpr T TWO_PI = 6.28318531; // might already be defined in Arduino.h uint_fast16_t _samples;
#endif T _samplingFrequency;
static constexpr T FOUR_PI = 12.56637061; T *_vImag;
static constexpr T SIX_PI = 18.84955593; T *_vReal;
FFTWindow _windowFunction;
static inline void Swap(T &x, T &y) /* Functions */
{ uint_fast8_t exponent(uint_fast16_t value) const;
T temp = x; void findMaxY(T *vData, uint_fast16_t length, T *maxY,
x = y; uint_fast16_t *index) const;
y = temp; 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 #ifdef FFT_SQRT_APPROXIMATION
// Fast inverse square root aka "Quake 3 fast inverse square root", multiplied by x. float sqrt_internal(float x) const;
// Uses one iteration of Halley's method for precision. double sqrt_internal(double x) const;
// See: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Iterative_methods_for_reciprocal_square_roots
// And: https://github.com/HorstBaerbel/approx
template <typename V = T>
static inline V sqrt_internal(typename std::enable_if<std::is_same<V, float>::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 <typename V = T>
static inline V sqrt_internal(typename std::enable_if<std::is_same<V, double>::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
}
#endif #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__ #if defined(__AVR__) && defined(USE_AVR_PROGMEM)
template <typename T> static const float _c1[] PROGMEM = {
const float ArduinoFFT<T>::_c1[] PROGMEM = { 0.0000000000, 0.7071067812, 0.9238795325, 0.9807852804, 0.9951847267,
0.0000000000, 0.7071067812, 0.9238795325, 0.9807852804, 0.9987954562, 0.9996988187, 0.9999247018, 0.9999811753, 0.9999952938,
0.9951847267, 0.9987954562, 0.9996988187, 0.9999247018, 0.9999988235, 0.9999997059, 0.9999999265, 0.9999999816, 0.9999999954,
0.9999811753, 0.9999952938, 0.9999988235, 0.9999997059, 0.9999999989, 0.9999999997};
0.9999999265, 0.9999999816, 0.9999999954, 0.9999999989, static const float _c2[] PROGMEM = {
0.9999999997}; 1.0000000000, 0.7071067812, 0.3826834324, 0.1950903220, 0.0980171403,
0.0490676743, 0.0245412285, 0.0122715383, 0.0061358846, 0.0030679568,
template <typename T> 0.0015339802, 0.0007669903, 0.0003834952, 0.0001917476, 0.0000958738,
const float ArduinoFFT<T>::_c2[] PROGMEM = { 0.0000479369, 0.0000239684};
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 #endif
template <typename T>
const T ArduinoFFT<T>::_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 #endif