diff --git a/README.md b/README.md index 1a93a22..577684e 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,10 @@ Tested on Arduino 1.6.11 ### 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 + To install this library, just place this entire folder as a subfolder in your Arduino installation When installed, this library should look like: @@ -37,23 +41,28 @@ select arduinoFTT. This will add a corresponding line to the top of your sketch * Document windowing functions advantages and disadvantages. * Optimize usage and arguments. * Add new windowing functions. -* Spectrum table? +* Spectrum table? ### API * **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); -Calculates the power value according to **Exponent** and calcuates the Fast Fourier Transform. * **Compute**(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir); +* **Compute**(uint8_t dir); Calcuates the Fast Fourier Transform. * **MajorPeak**(double *vD, uint16_t samples, double samplingFrequency); +* **MajorPeak**(); +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 diff --git a/changeLog.txt b/changeLog.txt index 3585e25..287cbb6 100644 --- a/changeLog.txt +++ b/changeLog.txt @@ -1,3 +1,9 @@ +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. diff --git a/library.json b/library.json new file mode 100644 index 0000000..cf9dd7a --- /dev/null +++ b/library.json @@ -0,0 +1,26 @@ +{ + "name": "arduinoFFT", + "keywords": "FFT, Fourier, FDT, frequency", + "description": "A library for implementing floating point Fast Fourier Transform calculations.", + "repository": + { + "type": "git", + "url": "https://github.com/kosme/arduinoFFT.git" + }, + "authors": + [ + { + "name": "Enrique Condes", + "email": "enrique@shapeoko.com", + "maintainer": true + }, + { + "name": "Didier Longueville", + "url": "http://www.arduinoos.com/", + "email": "contact@arduinoos.com" + } + ], + "version": "1.4", + "frameworks": ["arduino","mbed","espidf"], + "platforms": "*" +} diff --git a/library.properties b/library.properties index 13c1da4..4342ee6 100644 --- a/library.properties +++ b/library.properties @@ -1,6 +1,6 @@ name=arduinoFFT -version=1.3 -author=kosme +version=1.4 +author=Enrique Condes maintainer=Enrique Condes sentence=A library for implementing floating point Fast Fourier Transform calculations on Arduino. paragraph=With this library you can calculate the frequency of a sampled signal. diff --git a/src/arduinoFFT.cpp b/src/arduinoFFT.cpp index 29db02e..37d47d9 100644 --- a/src/arduinoFFT.cpp +++ b/src/arduinoFFT.cpp @@ -22,13 +22,22 @@ #include "arduinoFFT.h" arduinoFFT::arduinoFFT(void) -{ -/* Constructor */ +{ // Constructor + #warning("This method is deprecated and will be removed on future revisions.") +} + +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) { -/* Destructor */ +// Destructor } uint8_t arduinoFFT::Revision(void) @@ -38,18 +47,75 @@ uint8_t arduinoFFT::Revision(void) void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t dir) { + #warning("This method is deprecated and will be removed on future revisions.") Compute(vReal, vImag, samples, Exponent(samples), dir); } +void arduinoFFT::Compute(uint8_t 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 / + 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; + } + c2 = sqrt((1.0 - c1) / 2.0); + if (dir == FFT_FORWARD) { + c2 = -c2; + } + c1 = sqrt((1.0 + c1) / 2.0); + } + // Scaling for reverse transform / + if (dir != FFT_FORWARD) { + for (uint16_t i = 0; i < this->_samples; i++) { + this->_vReal[i] /= this->_samples; + this->_vImag[i] /= this->_samples; + } + } +} + void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir) -{ -/* Computes in-place complex-to-complex FFT */ - /* Reverse bits */ +{ // Computes in-place complex-to-complex FFT + // Reverse bits + #warning("This method is deprecated and will be removed on future revisions.") uint16_t j = 0; for (uint16_t i = 0; i < (samples - 1); i++) { if (i < j) { Swap(&vReal[i], &vReal[j]); - Swap(&vImag[i], &vImag[j]); + if(dir==FFT_REVERSE) + Swap(&vImag[i], &vImag[j]); } uint16_t k = (samples >> 1); while (k <= j) { @@ -58,7 +124,7 @@ void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t } j += k; } - /* Compute the FFT */ + // Compute the FFT double c1 = -1.0; double c2 = 0.0; uint16_t l2 = 1; @@ -87,7 +153,7 @@ void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t } c1 = sqrt((1.0 + c1) / 2.0); } - /* Scaling for reverse transform */ + // Scaling for reverse transform if (dir != FFT_FORWARD) { for (uint16_t i = 0; i < samples; i++) { vReal[i] /= samples; @@ -96,44 +162,95 @@ void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t } } +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])); + } +} + void arduinoFFT::ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples) -{ -/* vM is half the size of vReal and vImag */ +{ // vM is half the size of vReal and vImag + #warning("This method is deprecated and will be removed on future revisions.") for (uint16_t i = 0; i < samples; i++) { vReal[i] = sqrt(sq(vReal[i]) + sq(vImag[i])); } } +void arduinoFFT::Windowing(uint8_t windowType, uint8_t dir) +{// Weighing factors are computed once before multiple use of FFT +// The weighing function is symetric; 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) + weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne); + break; + case FFT_WIN_TYP_BLACKMAN: // blackmann + weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * 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; + } + } +} + + void arduinoFFT::Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir) -{ -/* Weighing factors are computed once before multiple use of FFT */ -/* The weighing function is symetric; half the weighs are recorded */ +{// Weighing factors are computed once before multiple use of FFT +// The weighing function is symetric; half the weighs are recorded + #warning("This method is deprecated and will 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 */ + // Compute and record weighting factor switch (windowType) { - case FFT_WIN_TYP_RECTANGLE: /* rectangle (box car) */ + case FFT_WIN_TYP_RECTANGLE: // rectangle (box car) weighingFactor = 1.0; break; - case FFT_WIN_TYP_HAMMING: /* hamming */ + case FFT_WIN_TYP_HAMMING: // hamming weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio)); break; - case FFT_WIN_TYP_HANN: /* hann */ + case FFT_WIN_TYP_HANN: // hann weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio)); break; - case FFT_WIN_TYP_TRIANGLE: /* triangle (Bartlett) */ + case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett) weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne); break; - case FFT_WIN_TYP_BLACKMAN: /* blackmann */ + case FFT_WIN_TYP_BLACKMAN: // blackmann weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * ratio))); 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)); break; - case FFT_WIN_TYP_WELCH: /* welch */ + case FFT_WIN_TYP_WELCH: // welch weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0)); break; } @@ -148,31 +265,36 @@ void arduinoFFT::Windowing(double *vData, uint16_t samples, uint8_t windowType, } } -void arduinoFFT::PrintVector(double *vData, uint16_t samples, double samplingFrequency) +double arduinoFFT::MajorPeak() { - PrintArray(vData,samples, samplingFrequency, SCL_INDEX); -} - -void arduinoFFT::PrintSignal(double *vData, uint16_t samples, double samplingFrequency) -{ - PrintArray(vData,samples, samplingFrequency, SCL_TIME); -} - -void arduinoFFT::PrintSpectrum(double *vData, uint16_t samples, double samplingFrequency) -{ - PrintArray(vData,samples, samplingFrequency, SCL_FREQUENCY); -} - -void arduinoFFT::PlotSpectrum(double *vData, uint16_t samples, double samplingFrequency) -{ - PrintArray(vData,samples, samplingFrequency, SCL_PLOT); + 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); + // retuned value: interpolated frequency peak apex + return(interpolatedX); } double arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequency) { + #warning("This method is deprecated and will be removed on future revisions.") double maxY = 0; uint16_t IndexOfMaxY = 0; - for (uint16_t i = 1; i < ((samples >> 1) - 1); i++) { + //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]; @@ -182,19 +304,22 @@ 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 interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples-1); - /* retuned value: interpolated frequency peak apex */ + if(IndexOfMaxY==(samples >> 1)) //To improve calculation on edge values + interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples); + // returned value: interpolated frequency peak apex return(interpolatedX); } uint8_t arduinoFFT::Exponent(uint16_t value) { - /* Calculates the base 2 logarithm of a value */ + #warning("This method will not be accessible on future revisions.") + // Calculates the base 2 logarithm of a value uint8_t result = 0; while (((value >> result) & 1) != 1) result++; return(result); } -/* Private functions */ +// Private functions void arduinoFFT::Swap(double *x, double *y) { @@ -202,35 +327,3 @@ void arduinoFFT::Swap(double *x, double *y) *x = *y; *y = temp; } - -void arduinoFFT::PrintArray(double *vData, uint16_t samples, double samplingFrequency, uint8_t scaleType) -{ - uint16_t bufferSize = samples; - if((scaleType == SCL_FREQUENCY)||(scaleType == SCL_PLOT)) - bufferSize = bufferSize>>1; - for (uint16_t i = 0; i < bufferSize; i++) - { - double abscissa; - switch (scaleType) - { - case SCL_INDEX: - abscissa = (i * 1.0); - break; - case SCL_TIME: - abscissa = ((i * 1.0) / samplingFrequency); - break; - case SCL_FREQUENCY: - case SCL_PLOT: - abscissa = ((i * 1.0 * samplingFrequency) / samples); - break; - } - if(scaleType!=SCL_PLOT){ - Serial.print(abscissa, 6); - if(scaleType==SCL_FREQUENCY) - Serial.print(" Hz"); - Serial.print(" "); - } - Serial.println(vData[i], 4); - } - Serial.println(); -} diff --git a/src/arduinoFFT.h b/src/arduinoFFT.h index 61aefde..e65eb17 100644 --- a/src/arduinoFFT.h +++ b/src/arduinoFFT.h @@ -42,10 +42,7 @@ /* Custom constants */ #define FFT_FORWARD 0x01 #define FFT_REVERSE 0x00 -#define SCL_INDEX 0x00 -#define SCL_TIME 0x01 -#define SCL_FREQUENCY 0x02 -#define SCL_PLOT 0x03 + /* Windowing type */ #define FFT_WIN_TYP_RECTANGLE 0x00 /* rectangle (Box car) */ #define FFT_WIN_TYP_HAMMING 0x01 /* hamming */ @@ -62,26 +59,31 @@ 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, uint8_t dir); void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir); - void PrintVector(double *vData, uint16_t samples, double samplingFrequency); - void PrintSignal(double *vData, uint16_t samples, double samplingFrequency); - void PrintSpectrum(double *vData, uint16_t samples, double samplingFrequency); - void PlotSpectrum(double *vData, uint16_t samples, double samplingFrequency); double MajorPeak(double *vD, uint16_t samples, double samplingFrequency); - uint8_t Revision(void); void Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir); - uint8_t Exponent(uint16_t value); + void ComplexToMagnitude(); + void Compute(uint8_t dir); + double MajorPeak(); + void Windowing(uint8_t windowType, uint8_t dir); private: + /* Variables */ + uint16_t _samples; + double _samplingFrequency; + double *_vReal; + double *_vImag; + uint8_t _power; /* Functions */ void Swap(double *x, double *y); - void PrintArray(double *vData, uint16_t samples, double samplingFrequency, uint8_t scaleType); - }; #endif