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

Merge pull request #91 from kosme/develop 

Version 2.0 from develop branch
master v2.0
Enrique Condes 2024-03-06 14:29:10 +08:00 zatwierdzone przez GitHub
rodzic bb90780ace 909240fe79
commit 3d17cac4c9
Nie znaleziono w bazie danych klucza dla tego podpisu
ID klucza GPG: B5690EEEBB952194
14 zmienionych plików z 765 dodań i 711 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

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

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

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

@ -1,10 +1,12 @@
/*
Example of use of the FFT library to compute FFT for several signals over a range of frequencies.
The exponent is calculated once before the execution since it is a constant.
This saves resources during the execution of the sketch and reduces the compiled size.
The sketch shows the time that the computation takes.
Copyright (C) 2014 Enrique Condes
Example of use of the FFT libray to compute FFT for several signals over a range of frequencies.
The exponent is calculated once before the excecution since it is a constant.
This saves resources during the excecution of the sketch and reduces the compiled size.
The sketch shows the time that the computing is taking.
Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@ -23,11 +25,9 @@
#include "arduinoFFT.h"
arduinoFFT FFT;
/*
These values can be changed in order to evaluate the functions
*/
const uint16_t samples = 64;
const double sampling = 40;
const uint8_t amplitude = 4;
@ -42,7 +42,10 @@ Input vectors receive computed results from FFT
double vReal[samples];
double vImag[samples];
unsigned long time;
/* Create FFT object */
ArduinoFFT<double> FFT = ArduinoFFT<double>(vReal, vImag, samples, sampling);
unsigned long startTime;
#define SCL_INDEX 0x00
#define SCL_TIME 0x01
@ -71,25 +74,24 @@ void loop()
}
/*Serial.println("Data:");
PrintVector(vReal, samples, SCL_TIME);*/
time=millis();
FFT = arduinoFFT(vReal, vImag, samples, sampling); /* Create FFT object */
FFT.Windowing(FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */
startTime=millis();
FFT.windowing(FFTWindow::Hamming, FFTDirection::Forward); /* Weigh data */
/*Serial.println("Weighed data:");
PrintVector(vReal, samples, SCL_TIME);*/
FFT.Compute(FFT_FORWARD); /* Compute FFT */
FFT.compute(FFTDirection::Forward); /* Compute FFT */
/*Serial.println("Computed Real values:");
PrintVector(vReal, samples, SCL_INDEX);
Serial.println("Computed Imaginary values:");
PrintVector(vImag, samples, SCL_INDEX);*/
FFT.ComplexToMagnitude(); /* Compute magnitudes */
FFT.complexToMagnitude(); /* Compute magnitudes */
/*Serial.println("Computed magnitudes:");
PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);*/
double x = FFT.MajorPeak();
double x = FFT.majorPeak();
Serial.print(frequency);
Serial.print(": \t\t");
Serial.print(x, 4);
Serial.print("\t\t");
Serial.print(millis()-time);
Serial.print(millis()-startTime);
Serial.println(" ms");
// delay(2000); /* Repeat after delay */
}

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

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

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

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

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

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

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

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

79
README.md
Wyświetl plik

@ -1,87 +1,38 @@
arduinoFFT
==========
Fast Fourier Transform for Arduino
# Fast Fourier Transform for Arduino
This is a fork from https://code.google.com/p/makefurt/ which has been abandoned since 2011.
<del>This is a C++ library for Arduino for computing FFT.</del> Now it works both on Arduino and C projects.
This is version 2.0 of the library, which has a different [API](#api).
Tested on Arduino 1.6.11
Tested on Arduino 1.8.19 and 2.3.2.
### Installation on Arduino
## Installation on Arduino
Use the Arduino Library Manager to install and keep it updated. Just look for arduinoFFT. Only for Arduino 1.5+
### Manual installation on Arduino
## Manual installation on Arduino
To install this library, just place this entire folder as a subfolder in your Arduino installation
To install this library, just place this entire folder as a subfolder in your Arduino installation. When installed, this library should look like:
When installed, this library should look like:
`Arduino\libraries\arduinoFTT` (this library's folder)
`Arduino\libraries\arduinoFTT\src\arduinoFTT.h` (the library header file. include this in your project)
`Arduino\libraries\arduinoFTT\keywords.txt` (the syntax coloring file)
`Arduino\libraries\arduinoFTT\Examples` (the examples in the "open" menu)
`Arduino\libraries\arduinoFTT\LICENSE` (GPL license file)
`Arduino\libraries\arduinoFTT\README.md` (this file)
```
Arduino\libraries\arduinoFTT (this library's folder)
Arduino\libraries\arduinoFTT\arduinoFTT.cpp (the library implementation file, uses 32 bits floats vectors)
Arduino\libraries\arduinoFTT\arduinoFTT.h (the library header file, uses 32 bits floats vectors)
Arduino\libraries\arduinoFTT\keywords.txt (the syntax coloring file)
Arduino\libraries\arduinoFTT\examples (the examples in the "open" menu)
Arduino\libraries\arduinoFTT\readme.md (this file)
```
### Building on Arduino
## Building on Arduino
After this library is installed, you just have to start the Arduino application.
You may see a few warning messages as it's built.
To use this library in a sketch, go to the Sketch | Import Library menu and
select arduinoFTT. This will add a corresponding line to the top of your sketch:
`#include <arduinoFTT.h>`
### TODO
* Ratio table for windowing function.
* Document windowing functions advantages and disadvantages.
* Optimize usage and arguments.
* Add new windowing functions.
<del>* Spectrum table? </del>
## API
### API
The exclamation mark `!` denotes that this method is deprecated and may be removed on future revisions.
* **!arduinoFFT**(void);
* **arduinoFFT**(double *vReal, double *vImag, uint16_t samples, double samplingFrequency);
Constructor
* **~arduinoFFT**(void);
Destructor
* **!ComplexToMagnitude**(double *vReal, double *vImag, uint16_t samples);
* **ComplexToMagnitude**();
* **!Compute**(double *vReal, double *vImag, uint16_t samples, uint8_t dir);
* **!Compute**(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir);
* **Compute**(uint8_t dir);
Calculates the Fast Fourier Transform.
* **!DCRemoval**(double *vData, uint16_t samples);
* **DCRemoval**();
Removes the DC component from the sample data.
* **!MajorPeak**(double *vD, uint16_t samples, double samplingFrequency);
* **!MajorPeak**(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v);
* **MajorPeak**();
* **MajorPeak**(double *f, double *v);
* **MajorPeakParabola**();
Looks for and returns the frequency of the biggest spike in the analyzed signal.
* **Revision**(void);
Returns the library revision.
* **!Windowing**(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir);
* **Windowing**(uint8_t windowType, uint8_t dir);
Performs a windowing function on the values array. The possible windowing options are:
* FFT_WIN_TYP_RECTANGLE
* FFT_WIN_TYP_HAMMING
* FFT_WIN_TYP_HANN
* FFT_WIN_TYP_TRIANGLE
* FFT_WIN_TYP_NUTTALL
* FFT_WIN_TYP_BLACKMAN
* FFT_WIN_TYP_BLACKMAN_NUTTALL
* FFT_WIN_TYP_BLACKMAN_HARRIS
* FFT_WIN_TYP_FLT_TOP
* FFT_WIN_TYP_WELCH
* **!Exponent**(uint16_t value);
Calculates and returns the base 2 logarithm of the given value.
Documentation was moved to the project's [wiki](https://github.com/kosme/arduinoFFT/wiki).

35
changeLog.txt
Wyświetl plik

@ -1,35 +0,0 @@
28/12/23 v1.6.2
Fix issue 52 and symplify calculation of synthetic data for examples 1, 2, 4, and 5.
25/07/23 v1.6.1
Code optimization for speed improvements. See issue 84 for more details.
03/09/23 v1.6
Include some functionality from development branch.
01/27/20 v1.5.5
Lookup table for constants c1 and c2 used during FFT comupting. This increases the FFT computing speed in around 5%.
02/10/18 v1.4
Transition version. Minor optimization to functions. New API. Deprecation of old functions.
12/06/18 v1.3
Add support for mbed development boards.
09/04/17 v1.2.3
Finally solves the issue of Arduino IDE not correctly detecting and highlighting the keywords.
09/03/17 v1.2.2
Solves a format issue in keywords.txt that prevented keywords from being detected.
08/28/17 v1.2.1
Fix to issues 6 and 7. Not cleaning the imaginary vector after each cycle leaded to erroneous calculations and could cause buffer overflows.
08/04/17 v1.2
Fix to bug preventing the number of samples to be greater than 128. New logical limit is 32768 samples but it is bound to the RAM on the chip.
05/12/17 v1.1
Fix issue that prevented installation through the Arduino Library Manager interface.
05/11/17 v1.0
Initial commit to Arduino Library Manager.

47
keywords.txt
Wyświetl plik

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

7
library.json
Wyświetl plik

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

2
library.properties
Wyświetl plik

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

1
src/.gitignore vendored 100644
Wyświetl plik

@ -0,0 +1 @@
/arduinoFFT.h.gch

899
src/arduinoFFT.cpp
Wyświetl plik

@ -1,529 +1,448 @@
/*
FFT library
Copyright (C) 2010 Didier Longueville
Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (template, speed improvements)
FFT library
Copyright (C) 2010 Didier Longueville
Copyright (C) 2014 Enrique Condes
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "arduinoFFT.h"
arduinoFFT::arduinoFFT(void) { // Constructor
#warning("This method is deprecated and may be removed on future revisions.")
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
}
arduinoFFT::arduinoFFT(double *vReal, double *vImag, uint16_t samples,
double samplingFrequency) { // Constructor
this->_vReal = vReal;
this->_vImag = vImag;
this->_samples = samples;
this->_samplingFrequency = samplingFrequency;
this->_power = Exponent(samples);
}
arduinoFFT::~arduinoFFT(void) {
template <typename T> ArduinoFFT<T>::~ArduinoFFT(void) {
// Destructor
}
uint8_t arduinoFFT::Revision(void) { return (FFT_LIB_REV); }
void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples,
FFTDirection dir) {
#warning("This method is deprecated and may be removed on future revisions.")
Compute(vReal, vImag, samples, Exponent(samples), dir);
}
void arduinoFFT::Compute(FFTDirection dir) {
// Computes in-place complex-to-complex FFT /
// Reverse bits /
uint16_t j = 0;
for (uint16_t i = 0; i < (this->_samples - 1); i++) {
if (i < j) {
Swap(&this->_vReal[i], &this->_vReal[j]);
if (dir == FFT_REVERSE)
Swap(&this->_vImag[i], &this->_vImag[j]);
}
uint16_t k = (this->_samples >> 1);
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
// Compute the FFT
#ifdef __AVR__
uint8_t index = 0;
#endif
double c1 = -1.0;
double c2 = 0.0;
uint16_t l2 = 1;
for (uint8_t l = 0; (l < this->_power); l++) {
uint16_t l1 = l2;
l2 <<= 1;
double u1 = 1.0;
double u2 = 0.0;
for (j = 0; j < l1; j++) {
for (uint16_t i = j; i < this->_samples; i += l2) {
uint16_t i1 = i + l1;
double t1 = u1 * this->_vReal[i1] - u2 * this->_vImag[i1];
double t2 = u1 * this->_vImag[i1] + u2 * this->_vReal[i1];
this->_vReal[i1] = this->_vReal[i] - t1;
this->_vImag[i1] = this->_vImag[i] - t2;
this->_vReal[i] += t1;
this->_vImag[i] += t2;
}
double z = ((u1 * c1) - (u2 * c2));
u2 = ((u1 * c2) + (u2 * c1));
u1 = z;
}
#ifdef __AVR__
c2 = pgm_read_float_near(&(_c2[index]));
c1 = pgm_read_float_near(&(_c1[index]));
index++;
#else
c2 = sqrt((1.0 - c1) / 2.0);
c1 = sqrt((1.0 + c1) / 2.0);
#endif
if (dir == FFT_FORWARD) {
c2 = -c2;
}
}
// Scaling for reverse transform /
if (dir != FFT_FORWARD) {
double reciprocal = 1.0 / this->_samples;
for (uint16_t i = 0; i < this->_samples; i++) {
this->_vReal[i] *= reciprocal;
this->_vImag[i] *= reciprocal;
}
if (_precompiledWindowingFactors) {
delete [] _precompiledWindowingFactors;
}
}
void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples,
uint8_t power, FFTDirection dir) {
template <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
// Reverse bits
#warning("This method is deprecated and may be removed on future revisions.")
uint16_t j = 0;
for (uint16_t i = 0; i < (samples - 1); i++) {
template <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 == FFT_REVERSE)
Swap(&vImag[i], &vImag[j]);
swap(&vReal[i], &vReal[j]);
if (dir == FFTDirection::Reverse)
swap(&vImag[i], &vImag[j]);
}
uint16_t k = (samples >> 1);
uint_fast16_t k = (samples >> 1);
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
// Compute the FFT
#ifdef __AVR__
uint8_t index = 0;
#endif
double c1 = -1.0;
double c2 = 0.0;
uint16_t l2 = 1;
for (uint8_t l = 0; (l < power); l++) {
uint16_t l1 = l2;
// Compute the FFT
T c1 = -1.0;
T c2 = 0.0;
uint_fast16_t l2 = 1;
for (uint_fast8_t l = 0; (l < power); l++) {
uint_fast16_t l1 = l2;
l2 <<= 1;
double u1 = 1.0;
double u2 = 0.0;
T u1 = 1.0;
T u2 = 0.0;
for (j = 0; j < l1; j++) {
for (uint16_t i = j; i < samples; i += l2) {
uint16_t i1 = i + l1;
double t1 = u1 * vReal[i1] - u2 * vImag[i1];
double t2 = u1 * vImag[i1] + u2 * vReal[i1];
for (uint_fast16_t i = j; i < samples; i += l2) {
uint_fast16_t i1 = i + l1;
T t1 = u1 * vReal[i1] - u2 * vImag[i1];
T t2 = u1 * vImag[i1] + u2 * vReal[i1];
vReal[i1] = vReal[i] - t1;
vImag[i1] = vImag[i] - t2;
vReal[i] += t1;
vImag[i] += t2;
}
double z = ((u1 * c1) - (u2 * c2));
T z = ((u1 * c1) - (u2 * c2));
u2 = ((u1 * c2) + (u2 * c1));
u1 = z;
}
#ifdef __AVR__
c2 = pgm_read_float_near(&(_c2[index]));
c1 = pgm_read_float_near(&(_c1[index]));
index++;
#if defined(__AVR__) && defined(USE_AVR_PROGMEM)
c2 = pgm_read_float_near(&(_c2[l]));
c1 = pgm_read_float_near(&(_c1[l]));
#else
c2 = sqrt((1.0 - c1) / 2.0);
c1 = sqrt((1.0 + c1) / 2.0);
T cTemp = 0.5 * c1;
c2 = sqrt_internal(0.5 - cTemp);
c1 = sqrt_internal(0.5 + cTemp);
#endif
if (dir == FFT_FORWARD) {
if (dir == FFTDirection::Forward) {
c2 = -c2;
}
}
// Scaling for reverse transform
if (dir != FFT_FORWARD) {
double reciprocal = 1.0 / samples;
for (uint16_t i = 0; i < samples; i++) {
vReal[i] *= reciprocal;
vImag[i] *= reciprocal;
if (dir == FFTDirection::Reverse) {
for (uint_fast16_t i = 0; i < samples; i++) {
#ifdef FFT_SPEED_OVER_PRECISION
vReal[i] *= oneOverSamples;
vImag[i] *= oneOverSamples;
#else
vReal[i] /= samples;
vImag[i] /= samples;
#endif
}
}
}
void arduinoFFT::ComplexToMagnitude() {
// vM is half the size of vReal and vImag
for (uint16_t i = 0; i < this->_samples; i++) {
this->_vReal[i] = sqrt(sq(this->_vReal[i]) + sq(this->_vImag[i]));
}
template <typename T> void ArduinoFFT<T>::dcRemoval(void) const {
dcRemoval(this->_vReal, this->_samples);
}
void arduinoFFT::ComplexToMagnitude(double *vReal, double *vImag,
uint16_t samples) {
// vM is half the size of vReal and vImag
#warning("This method is deprecated and may be removed on future revisions.")
for (uint16_t i = 0; i < samples; i++) {
vReal[i] = sqrt(sq(vReal[i]) + sq(vImag[i]));
}
}
void arduinoFFT::DCRemoval() {
template <typename T>
void ArduinoFFT<T>::dcRemoval(T *vData, uint_fast16_t samples) const {
// calculate the mean of vData
double mean = 0;
for (uint16_t i = 0; i < this->_samples; i++) {
mean += this->_vReal[i];
}
mean /= this->_samples;
// Subtract the mean from vData
for (uint16_t i = 0; i < this->_samples; i++) {
this->_vReal[i] -= mean;
}
}
void arduinoFFT::DCRemoval(double *vData, uint16_t samples) {
// calculate the mean of vData
#warning("This method is deprecated and may be removed on future revisions.")
double mean = 0;
for (uint16_t i = 0; i < samples; i++) {
T mean = 0;
for (uint_fast16_t i = 0; i < samples; i++) {
mean += vData[i];
}
mean /= samples;
// Subtract the mean from vData
for (uint16_t i = 0; i < samples; i++) {
for (uint_fast16_t i = 0; i < samples; i++) {
vData[i] -= mean;
}
}
void arduinoFFT::Windowing(FFTWindow windowType, FFTDirection dir) {
// Weighing factors are computed once before multiple use of FFT
// The weighing function is symmetric; half the weighs are recorded
double samplesMinusOne = (double(this->_samples) - 1.0);
for (uint16_t i = 0; i < (this->_samples >> 1); i++) {
double indexMinusOne = double(i);
double ratio = (indexMinusOne / samplesMinusOne);
double weighingFactor = 1.0;
// Compute and record weighting factor
switch (windowType) {
case FFT_WIN_TYP_RECTANGLE: // rectangle (box car)
weighingFactor = 1.0;
break;
case FFT_WIN_TYP_HAMMING: // hamming
weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio));
break;
case FFT_WIN_TYP_HANN: // hann
weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio));
break;
case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
#if defined(ESP8266) || defined(ESP32)
weighingFactor =
1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#else
weighingFactor =
1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#endif
break;
case FFT_WIN_TYP_NUTTALL: // nuttall
weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) +
(0.144232 * (cos(fourPi * ratio))) -
(0.012604 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN: // blackman
weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) +
(0.07922 * (cos(fourPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall
weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) +
(0.1365995 * (cos(fourPi * ratio))) -
(0.0106411 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris
weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) +
(0.14128 * (cos(fourPi * ratio))) -
(0.01168 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_FLT_TOP: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) +
(0.1980399 * cos(fourPi * ratio));
break;
case FFT_WIN_TYP_WELCH: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) /
(samplesMinusOne / 2.0));
break;
}
if (dir == FFT_FORWARD) {
this->_vReal[i] *= weighingFactor;
this->_vReal[this->_samples - (i + 1)] *= weighingFactor;
} else {
this->_vReal[i] /= weighingFactor;
this->_vReal[this->_samples - (i + 1)] /= weighingFactor;
}
}
template <typename T> T ArduinoFFT<T>::majorPeak(void) const {
return majorPeak(this->_vReal, this->_samples, this->_samplingFrequency);
}
void arduinoFFT::Windowing(double *vData, uint16_t samples,
FFTWindow windowType, FFTDirection dir) {
// Weighing factors are computed once before multiple use of FFT
// The weighing function is symmetric; half the weighs are recorded
#warning("This method is deprecated and may be removed on future revisions.")
double samplesMinusOne = (double(samples) - 1.0);
for (uint16_t i = 0; i < (samples >> 1); i++) {
double indexMinusOne = double(i);
double ratio = (indexMinusOne / samplesMinusOne);
double weighingFactor = 1.0;
// Compute and record weighting factor
switch (windowType) {
case FFT_WIN_TYP_RECTANGLE: // rectangle (box car)
weighingFactor = 1.0;
break;
case FFT_WIN_TYP_HAMMING: // hamming
weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio));
break;
case FFT_WIN_TYP_HANN: // hann
weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio));
break;
case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
#if defined(ESP8266) || defined(ESP32)
weighingFactor =
1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#else
weighingFactor =
1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#endif
break;
case FFT_WIN_TYP_NUTTALL: // nuttall
weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) +
(0.144232 * (cos(fourPi * ratio))) -
(0.012604 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN: // blackman
weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) +
(0.07922 * (cos(fourPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall
weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) +
(0.1365995 * (cos(fourPi * ratio))) -
(0.0106411 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris
weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) +
(0.14128 * (cos(fourPi * ratio))) -
(0.01168 * (cos(sixPi * ratio)));
break;
case FFT_WIN_TYP_FLT_TOP: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) +
(0.1980399 * cos(fourPi * ratio));
break;
case FFT_WIN_TYP_WELCH: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) /
(samplesMinusOne / 2.0));
break;
}
if (dir == FFT_FORWARD) {
vData[i] *= weighingFactor;
vData[samples - (i + 1)] *= weighingFactor;
} else {
vData[i] /= weighingFactor;
vData[samples - (i + 1)] /= weighingFactor;
}
}
template <typename T> void ArduinoFFT<T>::majorPeak(T *f, T *v) const {
majorPeak(this->_vReal, this->_samples, this->_samplingFrequency, f, v);
}
double arduinoFFT::MajorPeak() {
double maxY = 0;
uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
if ((this->_vReal[i - 1] < this->_vReal[i]) &&
(this->_vReal[i] > this->_vReal[i + 1])) {
if (this->_vReal[i] > maxY) {
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
double delta =
0.5 *
((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) /
(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) +
this->_vReal[IndexOfMaxY + 1]));
double interpolatedX =
((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1);
if (IndexOfMaxY ==
(this->_samples >> 1)) // To improve calculation on edge values
interpolatedX =
((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples);
// returned value: interpolated frequency peak apex
return (interpolatedX);
template <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;
}
void arduinoFFT::MajorPeak(double *f, double *v) {
double maxY = 0;
uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
if ((this->_vReal[i - 1] < this->_vReal[i]) &&
(this->_vReal[i] > this->_vReal[i + 1])) {
if (this->_vReal[i] > maxY) {
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
double delta =
0.5 *
((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) /
(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) +
this->_vReal[IndexOfMaxY + 1]));
double interpolatedX =
((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1);
if (IndexOfMaxY ==
(this->_samples >> 1)) // To improve calculation on edge values
interpolatedX =
((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples);
// returned value: interpolated frequency peak apex
*f = interpolatedX;
#if defined(ESP8266) || defined(ESP32)
*v = fabs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) +
this->_vReal[IndexOfMaxY + 1]);
#else
*v = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) +
this->_vReal[IndexOfMaxY + 1]);
#endif
}
template <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);
double arduinoFFT::MajorPeak(double *vD, uint16_t samples,
double samplingFrequency) {
#warning("This method is deprecated and may be removed on future revisions.")
double maxY = 0;
uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2
for (uint16_t i = 1; i < ((samples >> 1) + 1); i++) {
if ((vD[i - 1] < vD[i]) && (vD[i] > vD[i + 1])) {
if (vD[i] > maxY) {
maxY = vD[i];
IndexOfMaxY = i;
}
}
}
double delta =
0.5 *
((vD[IndexOfMaxY - 1] - vD[IndexOfMaxY + 1]) /
(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]));
double interpolatedX =
((IndexOfMaxY + delta) * samplingFrequency) / (samples - 1);
T delta = 0.5 * ((vData[IndexOfMaxY - 1] - vData[IndexOfMaxY + 1]) /
(vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) +
vData[IndexOfMaxY + 1]));
T interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples - 1);
if (IndexOfMaxY == (samples >> 1)) // To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples);
// returned value: interpolated frequency peak apex
return (interpolatedX);
}
void arduinoFFT::MajorPeak(double *vD, uint16_t samples,
double samplingFrequency, double *f, double *v) {
#warning("This method is deprecated and may be removed on future revisions.")
double maxY = 0;
uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2
for (uint16_t i = 1; i < ((samples >> 1) + 1); i++) {
if ((vD[i - 1] < vD[i]) && (vD[i] > vD[i + 1])) {
if (vD[i] > maxY) {
maxY = vD[i];
IndexOfMaxY = i;
}
}
}
double delta =
0.5 *
((vD[IndexOfMaxY - 1] - vD[IndexOfMaxY + 1]) /
(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]));
double interpolatedX =
((IndexOfMaxY + delta) * samplingFrequency) / (samples - 1);
// double popo =
if (IndexOfMaxY == (samples >> 1)) // To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples);
// returned value: interpolated frequency peak apex
*f = interpolatedX;
*frequency = interpolatedX;
if (magnitude != nullptr) {
#if defined(ESP8266) || defined(ESP32)
*v =
fabs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]);
*magnitude = fabs(vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) +
vData[IndexOfMaxY + 1]);
#else
*v = abs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]);
*magnitude = abs(vData[IndexOfMaxY - 1] - (2.0 * vData[IndexOfMaxY]) +
vData[IndexOfMaxY + 1]);
#endif
}
}
double arduinoFFT::MajorPeakParabola() {
double maxY = 0;
uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
if ((this->_vReal[i - 1] < this->_vReal[i]) &&
(this->_vReal[i] > this->_vReal[i + 1])) {
if (this->_vReal[i] > maxY) {
maxY = this->_vReal[i];
IndexOfMaxY = i;
}
}
}
double freq = 0;
if (IndexOfMaxY > 0) {
// Assume the three points to be on a parabola
double a, b, c;
Parabola(IndexOfMaxY - 1, this->_vReal[IndexOfMaxY - 1], IndexOfMaxY,
this->_vReal[IndexOfMaxY], IndexOfMaxY + 1,
this->_vReal[IndexOfMaxY + 1], &a, &b, &c);
// Peak is at the middle of the parabola
double x = -b / (2 * a);
// And magnitude is at the extrema of the parabola if you want It...
// double y = a*x*x+b*x+c;
// Convert to frequency
freq = (x * this->_samplingFrequency) / (this->_samples);
}
template <typename T> T ArduinoFFT<T>::majorPeakParabola(void) const {
T freq = 0;
majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency,
&freq, nullptr);
return freq;
}
void arduinoFFT::Parabola(double x1, double y1, double x2, double y2, double x3,
double y3, double *a, double *b, double *c) {
double reversed_denom = 1 / ((x1 - x2) * (x1 - x3) * (x2 - x3));
template <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)) *
@ -533,19 +452,67 @@ void arduinoFFT::Parabola(double x1, double y1, double x2, double y2, double x3,
reversed_denom;
}
uint8_t arduinoFFT::Exponent(uint16_t value) {
#warning("This method may not be accessible on future revisions.")
// Calculates the base 2 logarithm of a value
uint8_t result = 0;
while (value >>= 1)
result++;
return (result);
template <typename T> void ArduinoFFT<T>::swap(T *a, T *b) const {
T temp = *a;
*a = *b;
*b = temp;
}
// Private functions
void arduinoFFT::Swap(double *x, double *y) {
double temp = *x;
*x = *y;
*y = temp;
#ifdef FFT_SQRT_APPROXIMATION
// Fast inverse square root aka "Quake 3 fast inverse square root", multiplied
// by x. Uses one iteration of Halley's method for precision. See:
// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Iterative_methods_for_reciprocal_square_roots
// And: https://github.com/HorstBaerbel/approx
template <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>;

171
src/arduinoFFT.h
Wyświetl plik

@ -1,21 +1,22 @@
/*
FFT library
Copyright (C) 2010 Didier Longueville
Copyright (C) 2014 Enrique Condes
FFT library
Copyright (C) 2010 Didier Longueville
Copyright (C) 2014 Enrique Condes
Copyright (C) 2020 Bim Overbohm (template, speed improvements)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
@ -38,22 +39,18 @@
#include "defs.h"
#include "types.h"
#include <math.h>
#include <stdint.h>
#endif
// Define this to use a low-precision square root approximation instead of the
// This definition uses a low-precision square root approximation instead of the
// regular sqrt() call
// This might only work for specific use cases, but is significantly faster.
// Only works for ArduinoFFT<float>.
// #define FFT_SQRT_APPROXIMATION
#ifdef FFT_SQRT_APPROXIMATION
#include <type_traits>
#else
#ifndef FFT_SQRT_APPROXIMATION
#define sqrt_internal sqrt
#endif
enum class FFTDirection { Reverse, Forward };
enum class FFTDirection { Forward, Reverse };
enum class FFTWindow {
Rectangle, // rectangle (Box car)
@ -65,10 +62,12 @@ enum class FFTWindow {
Blackman_Nuttall, // blackman nuttall
Blackman_Harris, // blackman harris
Flat_top, // flat top
Welch // welch
Welch, // welch
Precompiled // Placeholder for using custom or precompiled window values
};
#define FFT_LIB_REV 0x15
#define FFT_LIB_REV 0x20
/* Custom constants */
/* These defines keep compatibility with pre 2.0 code */
#define FFT_FORWARD FFTDirection::Forward
#define FFT_REVERSE FFTDirection::Reverse
@ -85,67 +84,95 @@ enum class FFTWindow {
FFTWindow::Blackman_Harris /* blackman harris*/
#define FFT_WIN_TYP_FLT_TOP FFTWindow::Flat_top /* flat top */
#define FFT_WIN_TYP_WELCH FFTWindow::Welch /* welch */
/*Mathematial constants*/
/* End of compatibility defines */
/* Mathematial constants */
#define twoPi 6.28318531
#define fourPi 12.56637061
#define sixPi 18.84955593
#ifdef __AVR__
static const double _c1[] PROGMEM = {
template <typename T> class ArduinoFFT {
public:
ArduinoFFT();
ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency,
bool windowingFactors = false);
~ArduinoFFT();
void complexToMagnitude(void) const;
void complexToMagnitude(T *vReal, T *vImag, uint_fast16_t samples) const;
void compute(FFTDirection dir) const;
void compute(T *vReal, T *vImag, uint_fast16_t samples,
FFTDirection dir) const;
void compute(T *vReal, T *vImag, uint_fast16_t samples, uint_fast8_t power,
FFTDirection dir) const;
void dcRemoval(void) const;
void dcRemoval(T *vData, uint_fast16_t samples) const;
T majorPeak(void) const;
void majorPeak(T *f, T *v) const;
T majorPeak(T *vData, uint_fast16_t samples, T samplingFrequency) const;
void majorPeak(T *vData, uint_fast16_t samples, T samplingFrequency,
T *frequency, T *magnitude) const;
T majorPeakParabola(void) const;
void majorPeakParabola(T *frequency, T *magnitude) const;
T majorPeakParabola(T *vData, uint_fast16_t samples,
T samplingFrequency) const;
void majorPeakParabola(T *vData, uint_fast16_t samples, T samplingFrequency,
T *frequency, T *magnitude) const;
uint8_t revision(void);
void setArrays(T *vReal, T *vImag, uint_fast16_t samples = 0);
void windowing(FFTWindow windowType, FFTDirection dir,
bool withCompensation = false);
void windowing(T *vData, uint_fast16_t samples, FFTWindow windowType,
FFTDirection dir, T *windowingFactors = nullptr,
bool withCompensation = false);
private:
/* Variables */
static const T _WindowCompensationFactors[10];
#ifdef FFT_SPEED_OVER_PRECISION
T _oneOverSamples = 0.0;
#endif
bool _isPrecompiled = false;
bool _precompiledWithCompensation = false;
uint_fast8_t _power = 0;
T *_precompiledWindowingFactors;
uint_fast16_t _samples;
T _samplingFrequency;
T *_vImag;
T *_vReal;
FFTWindow _windowFunction;
/* Functions */
uint_fast8_t exponent(uint_fast16_t value) const;
void findMaxY(T *vData, uint_fast16_t length, T *maxY,
uint_fast16_t *index) const;
void parabola(T x1, T y1, T x2, T y2, T x3, T y3, T *a, T *b, T *c) const;
void swap(T *a, T *b) const;
#ifdef FFT_SQRT_APPROXIMATION
float sqrt_internal(float x) const;
double sqrt_internal(double x) const;
#endif
};
#if defined(__AVR__) && defined(USE_AVR_PROGMEM)
static const float _c1[] PROGMEM = {
0.0000000000, 0.7071067812, 0.9238795325, 0.9807852804, 0.9951847267,
0.9987954562, 0.9996988187, 0.9999247018, 0.9999811753, 0.9999952938,
0.9999988235, 0.9999997059, 0.9999999265, 0.9999999816, 0.9999999954,
0.9999999989, 0.9999999997};
static const double _c2[] PROGMEM = {
static const float _c2[] PROGMEM = {
1.0000000000, 0.7071067812, 0.3826834324, 0.1950903220, 0.0980171403,
0.0490676743, 0.0245412285, 0.0122715383, 0.0061358846, 0.0030679568,
0.0015339802, 0.0007669903, 0.0003834952, 0.0001917476, 0.0000958738,
0.0000479369, 0.0000239684};
#endif
class arduinoFFT {
public:
/* Constructor */
arduinoFFT(void);
arduinoFFT(double *vReal, double *vImag, uint16_t samples,
double samplingFrequency);
/* Destructor */
~arduinoFFT(void);
/* Functions */
uint8_t Revision(void);
uint8_t Exponent(uint16_t value);
void ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples);
void Compute(double *vReal, double *vImag, uint16_t samples,
FFTDirection dir);
void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power,
FFTDirection dir);
void DCRemoval(double *vData, uint16_t samples);
double MajorPeak(double *vD, uint16_t samples, double samplingFrequency);
void MajorPeak(double *vD, uint16_t samples, double samplingFrequency,
double *f, double *v);
void Windowing(double *vData, uint16_t samples, FFTWindow windowType,
FFTDirection dir);
void ComplexToMagnitude();
void Compute(FFTDirection dir);
void DCRemoval();
double MajorPeak();
void MajorPeak(double *f, double *v);
void Windowing(FFTWindow windowType, FFTDirection dir);
double MajorPeakParabola();
private:
/* Variables */
uint16_t _samples;
double _samplingFrequency;
double *_vReal;
double *_vImag;
uint8_t _power;
/* Functions */
void Swap(double *x, double *y);
void Parabola(double x1, double y1, double x2, double y2, double x3,
double y3, double *a, double *b, double *c);
};
#endif