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

Bump to version 1.6

pull/81/head v1.6
kosme 2023-03-09 22:28:40 -06:00
rodzic 11b7937333
commit 94453e54ac
5 zmienionych plików z 579 dodań i 500 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

3
changeLog.txt
Wyświetl plik

@ -1,3 +1,6 @@
03/09/23 v1.6
Include some functionality from development branch.
01/27/20 v1.5.5 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%. Lookup table for constants c1 and c2 used during FFT comupting. This increases the FFT computing speed in around 5%.

2
library.json
Wyświetl plik

@ -20,7 +20,7 @@
"email": "contact@arduinoos.com" "email": "contact@arduinoos.com"
} }
], ],
"version": "1.5.6", "version": "1.6",
"frameworks": ["arduino","mbed","espidf"], "frameworks": ["arduino","mbed","espidf"],
"platforms": "*", "platforms": "*",
"headers": "arduinoFFT.h" "headers": "arduinoFFT.h"

2
library.properties
Wyświetl plik

@ -1,5 +1,5 @@
name=arduinoFFT name=arduinoFFT
version=1.5.6 version=1.6
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.

246
src/arduinoFFT.cpp
Wyświetl plik

@ -21,13 +21,12 @@
#include "arduinoFFT.h" #include "arduinoFFT.h"
arduinoFFT::arduinoFFT(void) arduinoFFT::arduinoFFT(void) { // Constructor
{ // Constructor
#warning("This method is deprecated and may be removed on future revisions.") #warning("This method is deprecated and may be removed on future revisions.")
} }
arduinoFFT::arduinoFFT(double *vReal, double *vImag, uint16_t samples, double samplingFrequency) arduinoFFT::arduinoFFT(double *vReal, double *vImag, uint16_t samples,
{// Constructor double samplingFrequency) { // Constructor
this->_vReal = vReal; this->_vReal = vReal;
this->_vImag = vImag; this->_vImag = vImag;
this->_samples = samples; this->_samples = samples;
@ -35,24 +34,20 @@ arduinoFFT::arduinoFFT(double *vReal, double *vImag, uint16_t samples, double sa
this->_power = Exponent(samples); this->_power = Exponent(samples);
} }
arduinoFFT::~arduinoFFT(void) arduinoFFT::~arduinoFFT(void) {
{
// Destructor // Destructor
} }
uint8_t arduinoFFT::Revision(void) uint8_t arduinoFFT::Revision(void) { return (FFT_LIB_REV); }
{
return(FFT_LIB_REV);
}
void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t dir) void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples,
{ FFTDirection dir) {
#warning("This method is deprecated and may be removed on future revisions.") #warning("This method is deprecated and may be removed on future revisions.")
Compute(vReal, vImag, samples, Exponent(samples), dir); Compute(vReal, vImag, samples, Exponent(samples), dir);
} }
void arduinoFFT::Compute(uint8_t dir) void arduinoFFT::Compute(FFTDirection dir) {
{// Computes in-place complex-to-complex FFT / // Computes in-place complex-to-complex FFT /
// Reverse bits / // Reverse bits /
uint16_t j = 0; uint16_t j = 0;
for (uint16_t i = 0; i < (this->_samples - 1); i++) { for (uint16_t i = 0; i < (this->_samples - 1); i++) {
@ -68,7 +63,7 @@ void arduinoFFT::Compute(uint8_t dir)
} }
j += k; j += k;
} }
// Compute the FFT / // Compute the FFT
#ifdef __AVR__ #ifdef __AVR__
uint8_t index = 0; uint8_t index = 0;
#endif #endif
@ -115,8 +110,9 @@ void arduinoFFT::Compute(uint8_t dir)
} }
} }
void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir) void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples,
{ // Computes in-place complex-to-complex FFT uint8_t power, FFTDirection dir) {
// Computes in-place complex-to-complex FFT
// Reverse bits // Reverse bits
#warning("This method is deprecated and may be removed on future revisions.") #warning("This method is deprecated and may be removed on future revisions.")
uint16_t j = 0; uint16_t j = 0;
@ -180,56 +176,51 @@ void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t
} }
} }
void arduinoFFT::ComplexToMagnitude() void arduinoFFT::ComplexToMagnitude() {
{ // vM is half the size of vReal and vImag // vM is half the size of vReal and vImag
for (uint16_t i = 0; i < this->_samples; i++) { for (uint16_t i = 0; i < this->_samples; i++) {
this->_vReal[i] = sqrt(sq(this->_vReal[i]) + sq(this->_vImag[i])); this->_vReal[i] = sqrt(sq(this->_vReal[i]) + sq(this->_vImag[i]));
} }
} }
void arduinoFFT::ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples) void arduinoFFT::ComplexToMagnitude(double *vReal, double *vImag,
{ // vM is half the size of vReal and 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.") #warning("This method is deprecated and may be removed on future revisions.")
for (uint16_t i = 0; i < samples; i++) { for (uint16_t i = 0; i < samples; i++) {
vReal[i] = sqrt(sq(vReal[i]) + sq(vImag[i])); vReal[i] = sqrt(sq(vReal[i]) + sq(vImag[i]));
} }
} }
void arduinoFFT::DCRemoval() void arduinoFFT::DCRemoval() {
{
// calculate the mean of vData // calculate the mean of vData
double mean = 0; double mean = 0;
for (uint16_t i = 0; i < this->_samples; i++) for (uint16_t i = 0; i < this->_samples; i++) {
{
mean += this->_vReal[i]; mean += this->_vReal[i];
} }
mean /= this->_samples; mean /= this->_samples;
// Subtract the mean from vData // Subtract the mean from vData
for (uint16_t i = 0; i < this->_samples; i++) for (uint16_t i = 0; i < this->_samples; i++) {
{
this->_vReal[i] -= mean; this->_vReal[i] -= mean;
} }
} }
void arduinoFFT::DCRemoval(double *vData, uint16_t samples) void arduinoFFT::DCRemoval(double *vData, uint16_t samples) {
{
// calculate the mean of vData // calculate the mean of vData
#warning("This method is deprecated and may be removed on future revisions.") #warning("This method is deprecated and may be removed on future revisions.")
double mean = 0; double mean = 0;
for (uint16_t i = 0; i < samples; i++) for (uint16_t i = 0; i < samples; i++) {
{
mean += vData[i]; mean += vData[i];
} }
mean /= samples; mean /= samples;
// Subtract the mean from vData // Subtract the mean from vData
for (uint16_t i = 0; i < samples; i++) for (uint16_t i = 0; i < samples; i++) {
{
vData[i] -= mean; vData[i] -= mean;
} }
} }
void arduinoFFT::Windowing(uint8_t windowType, uint8_t dir) void arduinoFFT::Windowing(FFTWindow windowType, FFTDirection dir) {
{// Weighing factors are computed once before multiple use of FFT // Weighing factors are computed once before multiple use of FFT
// The weighing function is symmetric; half the weighs are recorded // The weighing function is symmetric; half the weighs are recorded
double samplesMinusOne = (double(this->_samples) - 1.0); double samplesMinusOne = (double(this->_samples) - 1.0);
for (uint16_t i = 0; i < (this->_samples >> 1); i++) { for (uint16_t i = 0; i < (this->_samples >> 1); i++) {
@ -249,44 +240,56 @@ void arduinoFFT::Windowing(uint8_t windowType, uint8_t dir)
break; break;
case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett) case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
#if defined(ESP8266) || defined(ESP32) #if defined(ESP8266) || defined(ESP32)
weighingFactor = 1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne); weighingFactor =
1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#else #else
weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne); weighingFactor =
1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#endif #endif
break; break;
case FFT_WIN_TYP_NUTTALL: // nuttall case FFT_WIN_TYP_NUTTALL: // nuttall
weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) + (0.144232 * (cos(fourPi * ratio))) - (0.012604 * (cos(sixPi * ratio))); weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) +
(0.144232 * (cos(fourPi * ratio))) -
(0.012604 * (cos(sixPi * ratio)));
break; break;
case FFT_WIN_TYP_BLACKMAN: // blackman case FFT_WIN_TYP_BLACKMAN: // blackman
weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * ratio))); weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) +
(0.07922 * (cos(fourPi * ratio)));
break; break;
case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall 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))); weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) +
(0.1365995 * (cos(fourPi * ratio))) -
(0.0106411 * (cos(sixPi * ratio)));
break; break;
case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris 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))); weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) +
(0.14128 * (cos(fourPi * ratio))) -
(0.01168 * (cos(sixPi * ratio)));
break; break;
case FFT_WIN_TYP_FLT_TOP: // flat top case FFT_WIN_TYP_FLT_TOP: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) + (0.1980399 * cos(fourPi * ratio)); weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) +
(0.1980399 * cos(fourPi * ratio));
break; break;
case FFT_WIN_TYP_WELCH: // welch case FFT_WIN_TYP_WELCH: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0)); weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) /
(samplesMinusOne / 2.0));
break; break;
} }
if (dir == FFT_FORWARD) { if (dir == FFT_FORWARD) {
this->_vReal[i] *= weighingFactor; this->_vReal[i] *= weighingFactor;
this->_vReal[this->_samples - (i + 1)] *= weighingFactor; this->_vReal[this->_samples - (i + 1)] *= weighingFactor;
} } else {
else {
this->_vReal[i] /= weighingFactor; this->_vReal[i] /= weighingFactor;
this->_vReal[this->_samples - (i + 1)] /= weighingFactor; this->_vReal[this->_samples - (i + 1)] /= weighingFactor;
} }
} }
} }
void arduinoFFT::Windowing(double *vData, uint16_t samples,
void arduinoFFT::Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir) FFTWindow windowType, FFTDirection dir) {
{// Weighing factors are computed once before multiple use of FFT // Weighing factors are computed once before multiple use of FFT
// The weighing function is symetric; half the weighs are recorded // The weighing function is symetric; half the weighs are recorded
#warning("This method is deprecated and may be removed on future revisions.") #warning("This method is deprecated and may be removed on future revisions.")
double samplesMinusOne = (double(samples) - 1.0); double samplesMinusOne = (double(samples) - 1.0);
@ -307,92 +310,120 @@ void arduinoFFT::Windowing(double *vData, uint16_t samples, uint8_t windowType,
break; break;
case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett) case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
#if defined(ESP8266) || defined(ESP32) #if defined(ESP8266) || defined(ESP32)
weighingFactor = 1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne); weighingFactor =
1.0 - ((2.0 * fabs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#else #else
weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne); weighingFactor =
1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) /
samplesMinusOne);
#endif #endif
break; break;
case FFT_WIN_TYP_NUTTALL: // nuttall case FFT_WIN_TYP_NUTTALL: // nuttall
weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) + (0.144232 * (cos(fourPi * ratio))) - (0.012604 * (cos(sixPi * ratio))); weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) +
(0.144232 * (cos(fourPi * ratio))) -
(0.012604 * (cos(sixPi * ratio)));
break; break;
case FFT_WIN_TYP_BLACKMAN: // blackman case FFT_WIN_TYP_BLACKMAN: // blackman
weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * ratio))); weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) +
(0.07922 * (cos(fourPi * ratio)));
break; break;
case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall 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))); weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) +
(0.1365995 * (cos(fourPi * ratio))) -
(0.0106411 * (cos(sixPi * ratio)));
break; break;
case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris 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))); weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) +
(0.14128 * (cos(fourPi * ratio))) -
(0.01168 * (cos(sixPi * ratio)));
break; break;
case FFT_WIN_TYP_FLT_TOP: // flat top case FFT_WIN_TYP_FLT_TOP: // flat top
weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) + (0.1980399 * cos(fourPi * ratio)); weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) +
(0.1980399 * cos(fourPi * ratio));
break; break;
case FFT_WIN_TYP_WELCH: // welch case FFT_WIN_TYP_WELCH: // welch
weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0)); weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) /
(samplesMinusOne / 2.0));
break; break;
} }
if (dir == FFT_FORWARD) { if (dir == FFT_FORWARD) {
vData[i] *= weighingFactor; vData[i] *= weighingFactor;
vData[samples - (i + 1)] *= weighingFactor; vData[samples - (i + 1)] *= weighingFactor;
} } else {
else {
vData[i] /= weighingFactor; vData[i] /= weighingFactor;
vData[samples - (i + 1)] /= weighingFactor; vData[samples - (i + 1)] /= weighingFactor;
} }
} }
} }
double arduinoFFT::MajorPeak() double arduinoFFT::MajorPeak() {
{
double maxY = 0; double maxY = 0;
uint16_t IndexOfMaxY = 0; uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal, // If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2 // value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) { 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 - 1] < this->_vReal[i]) &&
(this->_vReal[i] > this->_vReal[i + 1])) {
if (this->_vReal[i] > maxY) { if (this->_vReal[i] > maxY) {
maxY = this->_vReal[i]; maxY = this->_vReal[i];
IndexOfMaxY = 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 delta =
double interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples-1); 0.5 *
if(IndexOfMaxY==(this->_samples >> 1)) //To improve calculation on edge values ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) /
interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples); (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 // returned value: interpolated frequency peak apex
return (interpolatedX); return (interpolatedX);
} }
void arduinoFFT::MajorPeak(double *f, double *v) void arduinoFFT::MajorPeak(double *f, double *v) {
{
double maxY = 0; double maxY = 0;
uint16_t IndexOfMaxY = 0; uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal, // If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2 // value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) { 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 - 1] < this->_vReal[i]) &&
(this->_vReal[i] > this->_vReal[i + 1])) {
if (this->_vReal[i] > maxY) { if (this->_vReal[i] > maxY) {
maxY = this->_vReal[i]; maxY = this->_vReal[i];
IndexOfMaxY = 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 delta =
double interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1); 0.5 *
if (IndexOfMaxY == (this->_samples >> 1)) //To improve calculation on edge values ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) /
interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples); (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 // returned value: interpolated frequency peak apex
*f = interpolatedX; *f = interpolatedX;
#if defined(ESP8266) || defined(ESP32) #if defined(ESP8266) || defined(ESP32)
*v = fabs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]); *v = fabs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) +
this->_vReal[IndexOfMaxY + 1]);
#else #else
*v = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]); *v = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) +
this->_vReal[IndexOfMaxY + 1]);
#endif #endif
} }
double arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequency) double arduinoFFT::MajorPeak(double *vD, uint16_t samples,
{ double samplingFrequency) {
#warning("This method is deprecated and may be removed on future revisions.") #warning("This method is deprecated and may be removed on future revisions.")
double maxY = 0; double maxY = 0;
uint16_t IndexOfMaxY = 0; uint16_t IndexOfMaxY = 0;
@ -406,16 +437,20 @@ double arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFreque
} }
} }
} }
double delta = 0.5 * ((vD[IndexOfMaxY-1] - vD[IndexOfMaxY+1]) / (vD[IndexOfMaxY-1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY+1])); double delta =
double interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples-1); 0.5 *
((vD[IndexOfMaxY - 1] - vD[IndexOfMaxY + 1]) /
(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]));
double interpolatedX =
((IndexOfMaxY + delta) * samplingFrequency) / (samples - 1);
if (IndexOfMaxY == (samples >> 1)) // To improve calculation on edge values if (IndexOfMaxY == (samples >> 1)) // To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples); interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples);
// returned value: interpolated frequency peak apex // returned value: interpolated frequency peak apex
return (interpolatedX); return (interpolatedX);
} }
void arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v) 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.") #warning("This method is deprecated and may be removed on future revisions.")
double maxY = 0; double maxY = 0;
uint16_t IndexOfMaxY = 0; uint16_t IndexOfMaxY = 0;
@ -429,32 +464,34 @@ void arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequenc
} }
} }
} }
double delta = 0.5 * ((vD[IndexOfMaxY - 1] - vD[IndexOfMaxY + 1]) / (vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1])); double delta =
double interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples - 1); 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 = // double popo =
if (IndexOfMaxY == (samples >> 1)) // To improve calculation on edge values if (IndexOfMaxY == (samples >> 1)) // To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples); interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples);
// returned value: interpolated frequency peak apex // returned value: interpolated frequency peak apex
*f = interpolatedX; *f = interpolatedX;
#if defined(ESP8266) || defined(ESP32) #if defined(ESP8266) || defined(ESP32)
*v = fabs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]); *v =
fabs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]);
#else #else
*v = abs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]); *v = abs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]);
#endif #endif
} }
double arduinoFFT::MajorPeakParabola() double arduinoFFT::MajorPeakParabola() {
{
double maxY = 0; double maxY = 0;
uint16_t IndexOfMaxY = 0; uint16_t IndexOfMaxY = 0;
// If sampling_frequency = 2 * max_frequency in signal, // If sampling_frequency = 2 * max_frequency in signal,
// value would be stored at position samples/2 // value would be stored at position samples/2
for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
{ if ((this->_vReal[i - 1] < this->_vReal[i]) &&
if ((this->_vReal[i-1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i+1])) (this->_vReal[i] > this->_vReal[i + 1])) {
{ if (this->_vReal[i] > maxY) {
if (this->_vReal[i] > maxY)
{
maxY = this->_vReal[i]; maxY = this->_vReal[i];
IndexOfMaxY = i; IndexOfMaxY = i;
} }
@ -462,11 +499,12 @@ double arduinoFFT::MajorPeakParabola()
} }
double freq = 0; double freq = 0;
if( IndexOfMaxY>0 ) if (IndexOfMaxY > 0) {
{
// Assume the three points to be on a parabola // Assume the three points to be on a parabola
double a, b, c; double a, b, c;
Parabola(IndexOfMaxY-1, this->_vReal[IndexOfMaxY-1], IndexOfMaxY, this->_vReal[IndexOfMaxY], IndexOfMaxY+1, this->_vReal[IndexOfMaxY+1], &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 // Peak is at the middle of the parabola
double x = -b / (2 * a); double x = -b / (2 * a);
@ -481,28 +519,30 @@ double arduinoFFT::MajorPeakParabola()
return freq; return freq;
} }
void arduinoFFT::Parabola(double x1, double y1, double x2, double y2, double x3, double y3, double *a, double *b, double *c) 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)); double reversed_denom = 1 / ((x1 - x2) * (x1 - x3) * (x2 - x3));
*a = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) * reversed_denom; *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; *b = (x3 * x3 * (y1 - y2) + x2 * x2 * (y3 - y1) + x1 * x1 * (y2 - y3)) *
*c = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3) *reversed_denom; reversed_denom;
*c = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 +
x1 * x2 * (x1 - x2) * y3) *
reversed_denom;
} }
uint8_t arduinoFFT::Exponent(uint16_t value) uint8_t arduinoFFT::Exponent(uint16_t value) {
{
#warning("This method may not be accessible on future revisions.") #warning("This method may not be accessible on future revisions.")
// Calculates the base 2 logarithm of a value // Calculates the base 2 logarithm of a value
uint8_t result = 0; uint8_t result = 0;
while (((value >> result) & 1) != 1) result++; while (((value >> result) & 1) != 1)
result++;
return (result); return (result);
} }
// Private functions // Private functions
void arduinoFFT::Swap(double *x, double *y) void arduinoFFT::Swap(double *x, double *y) {
{
double temp = *x; double temp = *x;
*x = *y; *x = *y;
*y = temp; *y = temp;

106
src/arduinoFFT.h
Wyświetl plik

@ -19,8 +19,8 @@
*/ */
#ifndef arduinoFFT_h /* Prevent loading library twice */ #ifndef ArduinoFFT_h /* Prevent loading library twice */
#define arduinoFFT_h #define ArduinoFFT_h
#ifdef ARDUINO #ifdef ARDUINO
#if ARDUINO >= 100 #if ARDUINO >= 100
#include "Arduino.h" #include "Arduino.h"
@ -28,55 +28,86 @@
#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>
#endif #endif
#define FFT_LIB_REV 0x14 // 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
#ifdef FFT_SQRT_APPROXIMATION
#include <type_traits>
#else
#define sqrt_internal sqrt
#endif
enum class FFTDirection { Reverse, Forward };
enum class FFTWindow {
Rectangle, // rectangle (Box car)
Hamming, // hamming
Hann, // hann
Triangle, // triangle (Bartlett)
Nuttall, // nuttall
Blackman, // blackman
Blackman_Nuttall, // blackman nuttall
Blackman_Harris, // blackman harris
Flat_top, // flat top
Welch // welch
};
#define FFT_LIB_REV 0x15
/* Custom constants */ /* Custom constants */
#define FFT_FORWARD 0x01 #define FFT_FORWARD FFTDirection::Forward
#define FFT_REVERSE 0x00 #define FFT_REVERSE FFTDirection::Reverse
/* Windowing type */ /* Windowing type */
#define FFT_WIN_TYP_RECTANGLE 0x00 /* rectangle (Box car) */ #define FFT_WIN_TYP_RECTANGLE FFTWindow::Rectangle /* rectangle (Box car) */
#define FFT_WIN_TYP_HAMMING 0x01 /* hamming */ #define FFT_WIN_TYP_HAMMING FFTWindow::Hamming /* hamming */
#define FFT_WIN_TYP_HANN 0x02 /* hann */ #define FFT_WIN_TYP_HANN FFTWindow::Hann /* hann */
#define FFT_WIN_TYP_TRIANGLE 0x03 /* triangle (Bartlett) */ #define FFT_WIN_TYP_TRIANGLE FFTWindow::Triangle /* triangle (Bartlett) */
#define FFT_WIN_TYP_NUTTALL 0x04 /* nuttall */ #define FFT_WIN_TYP_NUTTALL FFTWindow::Nuttall /* nuttall */
#define FFT_WIN_TYP_BLACKMAN 0x05 /* blackman */ #define FFT_WIN_TYP_BLACKMAN FFTWindow::Blackman /* blackman */
#define FFT_WIN_TYP_BLACKMAN_NUTTALL 0x06 /* blackman nuttall */ #define FFT_WIN_TYP_BLACKMAN_NUTTALL \
#define FFT_WIN_TYP_BLACKMAN_HARRIS 0x07 /* blackman harris*/ FFTWindow::Blackman_Nuttall /* blackman nuttall */
#define FFT_WIN_TYP_FLT_TOP 0x08 /* flat top */ #define FFT_WIN_TYP_BLACKMAN_HARRIS \
#define FFT_WIN_TYP_WELCH 0x09 /* welch */ 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*/ /*Mathematial constants*/
#define twoPi 6.28318531 #define twoPi 6.28318531
#define fourPi 12.56637061 #define fourPi 12.56637061
#define sixPi 18.84955593 #define sixPi 18.84955593
#ifdef __AVR__ #ifdef __AVR__
static const double _c1[]PROGMEM = {0.0000000000, 0.7071067812, 0.9238795325, 0.9807852804, static const double _c1[] PROGMEM = {
0.9951847267, 0.9987954562, 0.9996988187, 0.9999247018, 0.0000000000, 0.7071067812, 0.9238795325, 0.9807852804, 0.9951847267,
0.9999811753, 0.9999952938, 0.9999988235, 0.9999997059, 0.9987954562, 0.9996988187, 0.9999247018, 0.9999811753, 0.9999952938,
0.9999999265, 0.9999999816, 0.9999999954, 0.9999999989, 0.9999988235, 0.9999997059, 0.9999999265, 0.9999999816, 0.9999999954,
0.9999999997}; 0.9999999989, 0.9999999997};
static const double _c2[]PROGMEM = {1.0000000000, 0.7071067812, 0.3826834324, 0.1950903220, static const double _c2[] PROGMEM = {
0.0980171403, 0.0490676743, 0.0245412285, 0.0122715383, 1.0000000000, 0.7071067812, 0.3826834324, 0.1950903220, 0.0980171403,
0.0061358846, 0.0030679568, 0.0015339802, 0.0007669903, 0.0490676743, 0.0245412285, 0.0122715383, 0.0061358846, 0.0030679568,
0.0003834952, 0.0001917476, 0.0000958738, 0.0000479369, 0.0015339802, 0.0007669903, 0.0003834952, 0.0001917476, 0.0000958738,
0.0000239684}; 0.0000479369, 0.0000239684};
#endif #endif
class arduinoFFT { class arduinoFFT {
public: public:
/* Constructor */ /* Constructor */
arduinoFFT(void); arduinoFFT(void);
arduinoFFT(double *vReal, double *vImag, uint16_t samples, double samplingFrequency); arduinoFFT(double *vReal, double *vImag, uint16_t samples,
double samplingFrequency);
/* Destructor */ /* Destructor */
~arduinoFFT(void); ~arduinoFFT(void);
/* Functions */ /* Functions */
@ -84,19 +115,23 @@ public:
uint8_t Exponent(uint16_t value); uint8_t Exponent(uint16_t value);
void ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples); void ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples);
void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t dir); void Compute(double *vReal, double *vImag, uint16_t samples,
void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir); FFTDirection dir);
void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power,
FFTDirection dir);
void DCRemoval(double *vData, uint16_t samples); void DCRemoval(double *vData, uint16_t samples);
double MajorPeak(double *vD, uint16_t samples, double samplingFrequency); double MajorPeak(double *vD, uint16_t samples, double samplingFrequency);
void MajorPeak(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v); void MajorPeak(double *vD, uint16_t samples, double samplingFrequency,
void Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir); double *f, double *v);
void Windowing(double *vData, uint16_t samples, FFTWindow windowType,
FFTDirection dir);
void ComplexToMagnitude(); void ComplexToMagnitude();
void Compute(uint8_t dir); void Compute(FFTDirection dir);
void DCRemoval(); void DCRemoval();
double MajorPeak(); double MajorPeak();
void MajorPeak(double *f, double *v); void MajorPeak(double *f, double *v);
void Windowing(uint8_t windowType, uint8_t dir); void Windowing(FFTWindow windowType, FFTDirection dir);
double MajorPeakParabola(); double MajorPeakParabola();
@ -109,7 +144,8 @@ private:
uint8_t _power; uint8_t _power;
/* Functions */ /* Functions */
void Swap(double *x, double *y); 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); void Parabola(double x1, double y1, double x2, double y2, double x3,
double y3, double *a, double *b, double *c);
}; };
#endif #endif