Tinkercad
https://www.tinkercad.com/things/h9WtBi4SMII-fft
Código fonte com biblioteca FFT preparada para Tinkercad
// 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
//#define FFT_SQRT_APPROXIMATION
/*
These values can be changed in order to evaluate the functions
*/
#ifndef ArduinoFFT_h
#define ArduinoFFT_h
#ifdef ARDUINO
#if ARDUINO >= 100
#else
#include “WProgram.h”
#endif
#else
#include
#include
#ifdef __AVR__
#include
#include
#endif
#include “defs.h”
#include “types.h”
#include
#include
#endif
#ifndef TYPES_H
#define TYPES_H
#ifndef WIN32
#define FALSE 0
#define TRUE -1
#endif
typedef unsigned char u08;
typedef signed char s08;
typedef unsigned short u16;
typedef signed short s16;
typedef unsigned long u32;
typedef signed long s32;
typedef unsigned long long u64;
typedef signed long long s64;
#ifdef __MBED__
typedef unsigned char uint8_t;
typedef signed char int8_t;
typedef unsigned short uint16_t;
typedef signed short int16_t;
#endif
#define MAX_U08 255
#define MAX_U16 65535
#define MAX_U32 4294967295
#define MIN_S08 -128
#define MAX_S08 127
#define MIN_S16 -32768
#define MAX_S16 32767
#define MIN_S32 -2147483648
#define MAX_S32 2147483647
#ifndef WIN32
typedef unsigned char BOOL;
typedef unsigned char BYTE;
typedef unsigned int WORD;
typedef unsigned long DWORD;
typedef unsigned char UCHAR;
typedef unsigned int UINT;
typedef unsigned short USHORT;
typedef unsigned long ULONG;
typedef char CHAR;
typedef int INT;
typedef long LONG;
#endif
#endif
#ifndef enumsFFT_h
#define enumsFFT_h
#define FFT_FORWARD FFTDirection::Forward
#define FFT_REVERSE FFTDirection::Reverse
#define FFT_WIN_TYP_RECTANGLE FFTWindow::Rectangle /* rectangle (Box car) */
#define FFT_WIN_TYP_HAMMING FFTWindow::Hamming /* hamming */
#define FFT_WIN_TYP_HANN FFTWindow::Hann /* hann */
#define FFT_WIN_TYP_TRIANGLE FFTWindow::Triangle /* triangle (Bartlett) */
#define FFT_WIN_TYP_NUTTALL FFTWindow::Nuttall /* nuttall */
#define FFT_WIN_TYP_BLACKMAN FFTWindow::Blackman /* blackman */
#define FFT_WIN_TYP_BLACKMAN_NUTTALL \
FFTWindow::Blackman_Nuttall /* blackman nuttall */
#define FFT_WIN_TYP_BLACKMAN_HARRIS \
FFTWindow::Blackman_Harris /* blackman harris*/
#define FFT_WIN_TYP_FLT_TOP FFTWindow::Flat_top /* flat top */
#define FFT_WIN_TYP_WELCH FFTWindow::Welch /* welch */
#define twoPi 6.28318531
#define fourPi 12.56637061
#define sixPi 18.84955593
enum class FFTWindow {
Rectangle, // rectangle (Box car)
Hamming, // hamming
Hann, // hann
Triangle, // triangle (Bartlett)
Nuttall, // nuttall
Blackman, // blackman
Blackman_Nuttall, // blackman nuttall
Blackman_Harris, // blackman harris
Flat_top, // flat top
Welch, // welch
Precompiled // Placeholder for using custom or precompiled window values
};
enum class FFTDirection { Forward, Reverse };
#endif
#ifndef AVRLIBDEFS_H
#define AVRLIBDEFS_H
#define MEM_TYPE 1
#ifndef outb
#define outb(addr, data) addr = (data)
#endif
#ifndef inb
#define inb(addr) (addr)
#endif
#ifndef outw
#define outw(addr, data) addr = (data)
#endif
#ifndef inw
#define inw(addr) (addr)
#endif
#ifndef BV
#define BV(bit) (1<<(bit))
#endif
#ifndef cli
#define cli() __asm__ __volatile__ ("cli" ::)
#endif
#ifndef sei
#define sei() __asm__ __volatile__ ("sei" ::)
#endif
#ifdef __AVR_ATmega128__
#ifndef PD0
#endif
#endif
#define GNUC_PACKED __attribute__((packed))
#define DDR(x) ((x)-1) // address of data direction register of port x
#define PIN(x) ((x)-2) // address of input register of port x
#define MIN(a,b) ((a
#define ABS(x) ((x>0)?(x):(-x))
#define PI 3.14159265359
#define sq(x) ((x)*(x))
#define constrain(amt,low,high) ((amt)<(low)?(low):((amt)>(high)?(high):(amt)))
#endif
#ifndef FFT_SQRT_APPROXIMATION
#ifndef sqrt_internal
#define sqrt_internal sqrt
#endif
#endif
#define FFT_LIB_REV 0x20
template
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:
static const T _WindowCompensationFactors[11];
#ifdef FFT_SPEED_OVER_PRECISION
T _oneOverSamples = 0.0;
#endif
bool _isPrecompiled = false;
bool _precompiledWithCompensation = false;
uint_fast8_t _power = 0;
T *_precompiledWindowingFactors = nullptr;
uint_fast16_t _samples;
T _samplingFrequency;
T *_vImag;
T *_vReal;
FFTWindow _windowFunction;
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 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
#endif
template
template
ArduinoFFT
T samplingFrequency, bool windowingFactors)
: _samples(samples), _samplingFrequency(samplingFrequency), _vImag(vImag),
_vReal(vReal) {
if (windowingFactors) {
_precompiledWindowingFactors = new T[samples / 2];
}
_power = exponent(samples);
#ifdef FFT_SPEED_OVER_PRECISION
_oneOverSamples = 1.0 / samples;
#endif
}
template
if (_precompiledWindowingFactors) {
delete[] _precompiledWindowingFactors;
}
}
template
complexToMagnitude(this->_vReal, this->_vImag, this->_samples);
}
template
void ArduinoFFT
uint_fast16_t samples) const {
for (uint_fast16_t i = 0; i < samples; i++) {
vReal[i] = sqrt_internal(sq(vReal[i]) + sq(vImag[i]));
}
}
template
compute(this->_vReal, this->_vImag, this->_samples, exponent(this->_samples),
dir);
}
template
void ArduinoFFT
FFTDirection dir) const {
compute(vReal, vImag, samples, exponent(samples), dir);
}
template
void ArduinoFFT
uint_fast8_t power, FFTDirection dir) const {
#ifdef FFT_SPEED_OVER_PRECISION
T oneOverSamples = this->_oneOverSamples;
if (!this->_oneOverSamples)
oneOverSamples = 1.0 / samples;
#endif
uint_fast16_t j = 0;
for (uint_fast16_t i = 0; i < (samples - 1); i++) {
if (i < j) {
swap(&vReal[i], &vReal[j]);
if (dir == FFTDirection::Reverse)
swap(&vImag[i], &vImag[j]);
}
uint_fast16_t k = (samples >> 1);
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
T c1 = -1.0;
T c2 = 0.0;
uint_fast16_t l2 = 1;
for (uint_fast8_t l = 0; (l < power); l++) {
uint_fast16_t l1 = l2;
l2 <<= 1;
T u1 = 1.0;
T u2 = 0.0;
for (j = 0; j < l1; j++) {
for (uint_fast16_t i = j; i < samples; i += l2) {
uint_fast16_t i1 = i + l1;
T t1 = u1 * vReal[i1] - u2 * vImag[i1];
T t2 = u1 * vImag[i1] + u2 * vReal[i1];
vReal[i1] = vReal[i] - t1;
vImag[i1] = vImag[i] - t2;
vReal[i] += t1;
vImag[i] += t2;
}
T z = ((u1 * c1) - (u2 * c2));
u2 = ((u1 * c2) + (u2 * c1));
u1 = z;
}
#if defined(__AVR__) && defined(USE_AVR_PROGMEM)
c2 = pgm_read_float_near(&(_c2[l]));
c1 = pgm_read_float_near(&(_c1[l]));
#else
T cTemp = 0.5 * c1;
c2 = sqrt_internal(0.5 - cTemp);
c1 = sqrt_internal(0.5 + cTemp);
#endif
if (dir == FFTDirection::Forward) {
c2 = -c2;
}
}
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
}
}
vReal[0] = 0;
}
template
dcRemoval(this->_vReal, this->_samples);
}
template
void ArduinoFFT
T mean = 0;
for (uint_fast16_t i = 0; i < samples; i++) {
mean += vData[i];
}
mean /= samples;
for (uint_fast16_t i = 0; i < samples; i++) {
vData[i] -= mean;
}
}
template
return majorPeak(this->_vReal, this->_samples, this->_samplingFrequency);
}
template
majorPeak(this->_vReal, this->_samples, this->_samplingFrequency, f, v);
}
template
T ArduinoFFT
T samplingFrequency) const {
T frequency;
majorPeak(vData, samples, samplingFrequency, &frequency, nullptr);
return frequency;
}
template
void ArduinoFFT
T samplingFrequency, T *frequency,
T *magnitude) const {
T maxY = 0;
uint_fast16_t IndexOfMaxY = 0;
findMaxY(vData, (samples >> 1) + 1, &maxY, &IndexOfMaxY);
T delta = 0.5 * ((vData[IndexOfMaxY – 1] – vData[IndexOfMaxY + 1]) /
(vData[IndexOfMaxY – 1] – (2.0 * vData[IndexOfMaxY]) +
vData[IndexOfMaxY + 1]));
T interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples – 1);
if (IndexOfMaxY == (samples >> 1)) // To improve calculation on edge values
interpolatedX = ((IndexOfMaxY + delta) * samplingFrequency) / (samples);
*frequency = interpolatedX;
if (magnitude != nullptr) {
#if defined(ESP8266) || defined(ESP32)
*magnitude = fabs(vData[IndexOfMaxY – 1] – (2.0 * vData[IndexOfMaxY]) +
vData[IndexOfMaxY + 1]);
#else
*magnitude = abs(vData[IndexOfMaxY – 1] – (2.0 * vData[IndexOfMaxY]) +
vData[IndexOfMaxY + 1]);
#endif
}
}
template
T freq = 0;
majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency,
&freq, nullptr);
return freq;
}
template
void ArduinoFFT
majorPeakParabola(this->_vReal, this->_samples, this->_samplingFrequency,
frequency, magnitude);
}
template
T ArduinoFFT
T samplingFrequency) const {
T freq = 0;
majorPeakParabola(vData, samples, samplingFrequency, &freq, nullptr);
return freq;
}
template
void ArduinoFFT
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) {
T a, b, c;
parabola(IndexOfMaxY – 1, vData[IndexOfMaxY – 1], IndexOfMaxY,
vData[IndexOfMaxY], IndexOfMaxY + 1, vData[IndexOfMaxY + 1], &a,
&b, &c);
T x = -b / (2 * a);
if (magnitude != nullptr) {
*magnitude = (a * x * x) + (b * x) + c;
}
*frequency = (x * samplingFrequency) / samples;
}
}
template
return (FFT_LIB_REV);
}
template
void ArduinoFFT
_vReal = vReal;
_vImag = vImag;
if (samples) {
_samples = samples;
#ifdef FFT_SPEED_OVER_PRECISION
_oneOverSamples = 1.0 / samples;
#endif
if (_precompiledWindowingFactors) {
delete[] _precompiledWindowingFactors;
}
_precompiledWindowingFactors = new T[samples / 2];
}
}
template
void ArduinoFFT
bool withCompensation) {
if (this->_precompiledWindowingFactors && this->_isPrecompiled &&
this->_windowFunction == windowType &&
this->_precompiledWithCompensation == withCompensation) {
windowing(this->_vReal, this->_samples, FFTWindow::Precompiled, dir,
this->_precompiledWindowingFactors, withCompensation);
} else if (this->_precompiledWindowingFactors) {
windowing(this->_vReal, this->_samples, windowType, dir,
this->_precompiledWindowingFactors, withCompensation);
this->_isPrecompiled = true;
this->_precompiledWithCompensation = withCompensation;
this->_windowFunction = windowType;
} else {
windowing(this->_vReal, this->_samples, windowType, dir, nullptr,
withCompensation);
}
}
template
void ArduinoFFT
FFTWindow windowType, FFTDirection dir,
T *windowingFactors, bool withCompensation) {
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
}
for (uint_fast16_t i = 0; i < (samples >> 1); i++) {
T indexMinusOne = T(i);
T ratio = (indexMinusOne / samplesMinusOne);
T weighingFactor = 1.0;
switch (windowType) {
case FFTWindow::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:
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
}
}
}
}
template
uint_fast8_t ArduinoFFT
uint_fast8_t result = 0;
while (value >>= 1)
result++;
return result;
}
template
void ArduinoFFT
uint_fast16_t *index) const {
*maxY = 0;
*index = 0;
for (uint_fast16_t i = 1; i < length; i++) {
if ((vData[i - 1] < vData[i]) && (vData[i] > vData[i + 1])) {
if (vData[i] > vData[*index]) {
*index = i;
}
}
}
*maxY = vData[*index];
}
template
void ArduinoFFT
T *c) const {
const T reversed_denom = -0.5;
*a = (x3 * (y2 – y1) + x2 * (y1 – y3) + x1 * (y3 – y2)) * reversed_denom;
*b = (x3 * x3 * (y1 – y2) + x2 * x2 * (y3 – y1) + x1 * x1 * (y2 – y3)) *
reversed_denom;
*c = (x2 * x3 * (x2 – x3) * y1 + x3 * x1 * (x3 – x1) * y2 +
x1 * x2 * (x1 – x2) * y3) *
reversed_denom;
}
template
T temp = *a;
*a = *b;
*b = temp;
}
#ifdef FFT_SQRT_APPROXIMATION
template
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;
u.x = (0.125 * 3.0) * xu * (5.0 – xu2 * ((10.0 / 3.0) – xu2));
return u.x;
}
template
#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;
u.x = (0.125 * 3.0) * xu * (5.0 – xu2 * ((10.0 / 3.0) – xu2));
return u.x;
#endif
}
#endif
template
const T ArduinoFFT
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
1.0 // Custom, precompiled value.
};
template class ArduinoFFT
template class ArduinoFFT
#define CHANNEL A1
const uint16_t samples = 128; //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
#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02
#define SCL_PLOT 0x03
void setup()
{
analogReference(INTERNAL);
sampling_period_us = round(1000000*(1.0/samplingFrequency));
Serial.begin(115200);
Serial.println(“Ready”);
}
void loop()
{
/*SAMPLING*/
microseconds = micros();
for(int i=0; i
float x = FFT.majorPeak();
Serial.println(x, 6); //Print out what frequency is the most dominant.
Serial.println(x*60, 6); //Print out what frequency is the most dominant.
// while(1); /* Run Once */
delay(1000); /* 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();
}
Deixe um comentário