diff options
26 files changed, 2669 insertions, 204 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 574a0ad5..6f022e46 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -488,10 +488,8 @@ list(APPEND SOURCE_FILES drivers/TwiMaster.cpp heartratetask/HeartRateTask.cpp - components/heartrate/Ppg.cpp - components/heartrate/Biquad.cpp - components/heartrate/Ptagc.cpp components/heartrate/HeartRateController.cpp + components/heartrate/Ppg.cpp buttonhandler/ButtonHandler.cpp touchhandler/TouchHandler.cpp @@ -552,8 +550,7 @@ list(APPEND RECOVERY_SOURCE_FILES components/heartrate/HeartRateController.cpp heartratetask/HeartRateTask.cpp components/heartrate/Ppg.cpp - components/heartrate/Biquad.cpp - components/heartrate/Ptagc.cpp + components/motor/MotorController.cpp components/fs/FS.cpp buttonhandler/ButtonHandler.cpp @@ -666,9 +663,10 @@ set(INCLUDE_FILES drivers/TwiMaster.h heartratetask/HeartRateTask.h components/heartrate/Ppg.h - components/heartrate/Biquad.h - components/heartrate/Ptagc.h components/heartrate/HeartRateController.h + libs/arduinoFFT-develop/src/arduinoFFT.h + libs/arduinoFFT-develop/src/defs.h + libs/arduinoFFT-develop/src/types.h components/motor/MotorController.h buttonhandler/ButtonHandler.h touchhandler/TouchHandler.h diff --git a/src/components/heartrate/Biquad.cpp b/src/components/heartrate/Biquad.cpp deleted file mode 100644 index b7edd403..00000000 --- a/src/components/heartrate/Biquad.cpp +++ /dev/null @@ -1,26 +0,0 @@ -/* - SPDX-License-Identifier: LGPL-3.0-or-later - Original work Copyright (C) 2020 Daniel Thompson - C++ port Copyright (C) 2021 Jean-François Milants -*/ - -#include "components/heartrate/Biquad.h" - -using namespace Pinetime::Controllers; - -/** Original implementation from wasp-os : https://github.com/daniel-thompson/wasp-os/blob/master/wasp/ppg.py */ -Biquad::Biquad(float b0, float b1, float b2, float a1, float a2) : b0 {b0}, b1 {b1}, b2 {b2}, a1 {a1}, a2 {a2} { -} - -float Biquad::Step(float x) { - auto v1 = this->v1; - auto v2 = this->v2; - - auto v = x - (a1 * v1) - (a2 * v2); - auto y = (b0 * v) + (b1 * v1) + (b2 * v2); - - this->v2 = v1; - this->v1 = v; - - return y; -} diff --git a/src/components/heartrate/Biquad.h b/src/components/heartrate/Biquad.h deleted file mode 100644 index 7c8ca58f..00000000 --- a/src/components/heartrate/Biquad.h +++ /dev/null @@ -1,22 +0,0 @@ -#pragma once - -namespace Pinetime { - namespace Controllers { - /// Direct Form II Biquad Filter - class Biquad { - public: - Biquad(float b0, float b1, float b2, float a1, float a2); - float Step(float x); - - private: - float b0; - float b1; - float b2; - float a1; - float a2; - - float v1 = 0.0f; - float v2 = 0.0f; - }; - } -} diff --git a/src/components/heartrate/Ppg.cpp b/src/components/heartrate/Ppg.cpp index 900b1c22..3a6988ae 100644 --- a/src/components/heartrate/Ppg.cpp +++ b/src/components/heartrate/Ppg.cpp @@ -1,107 +1,292 @@ -/* - SPDX-License-Identifier: LGPL-3.0-or-later - Original work Copyright (C) 2020 Daniel Thompson - C++ port Copyright (C) 2021 Jean-François Milants -*/ - #include "components/heartrate/Ppg.h" -#include <vector> #include <nrf_log.h> +#include <vector> + using namespace Pinetime::Controllers; -/** Original implementation from wasp-os : https://github.com/daniel-thompson/wasp-os/blob/master/wasp/ppg.py */ namespace { - int Compare(int8_t* d1, int8_t* d2, size_t count) { - int e = 0; - for (size_t i = 0; i < count; i++) { - auto d = d1[i] - d2[i]; - e += d * d; + float LinearInterpolation(const float* xValues, const float* yValues, int length, float pointX) { + if (pointX > xValues[length - 1]) { + return yValues[length - 1]; + } else if (pointX <= xValues[0]) { + return yValues[0]; + } + int index = 0; + while (pointX > xValues[index] && index < length - 1) { + index++; + } + float pointX0 = xValues[index - 1]; + float pointX1 = xValues[index]; + float pointY0 = yValues[index - 1]; + float pointY1 = yValues[index]; + float mu = (pointX - pointX0) / (pointX1 - pointX0); + + return (pointY0 * (1 - mu) + pointY1 * mu); + } + + float PeakSearch(float* xVals, float* yVals, float threshold, float& width, float start, float end, int length) { + int peaks = 0; + bool enabled = false; + float minBin = 0.0f; + float maxBin = 0.0f; + float peakCenter = 0.0f; + float prevValue = LinearInterpolation(xVals, yVals, length, start - 0.01f); + float currValue = LinearInterpolation(xVals, yVals, length, start); + float idx = start; + while (idx < end) { + float nextValue = LinearInterpolation(xVals, yVals, length, idx + 0.01f); + if (currValue < threshold) { + enabled = true; + } + if (currValue >= threshold and enabled) { + if (prevValue < threshold) { + minBin = idx; + } else if (nextValue <= threshold) { + maxBin = idx; + peaks++; + width = maxBin - minBin; + peakCenter = width / 2.0f + minBin; + } + } + prevValue = currValue; + currValue = nextValue; + idx += 0.01f; } - return e; + if (peaks != 1) { + width = 0.0f; + peakCenter = 0.0f; + } + return peakCenter; } - int CompareShift(int8_t* d, int shift, size_t count) { - return Compare(d + shift, d, count - shift); + float SpectrumMean(const std::array<float, Ppg::spectrumLength>& signal, int start, int end) { + int total = 0; + float mean = 0.0f; + for (int idx = start; idx < end; idx++) { + mean += signal.at(idx); + total++; + } + if (total > 0) { + mean /= static_cast<float>(total); + } + return mean; } - int Trough(int8_t* d, size_t size, uint8_t mn, uint8_t mx) { - auto z2 = CompareShift(d, mn - 2, size); - auto z1 = CompareShift(d, mn - 1, size); - for (int i = mn; i < mx + 1; i++) { - auto z = CompareShift(d, i, size); - if (z2 > z1 && z1 < z) { - return i; + float SignalToNoise(const std::array<float, Ppg::spectrumLength>& signal, int start, int end, float max) { + float mean = SpectrumMean(signal, start, end); + return max / mean; + } + + // Simple bandpass filter using exponential moving average + void Filter30to240(std::array<float, Ppg::dataLength>& signal) { + // From: + // https://www.norwegiancreations.com/2016/03/arduino-tutorial-simple-high-pass-band-pass-and-band-stop-filtering/ + + int length = signal.size(); + // 0.268 is ~0.5Hz and 0.816 is ~4Hz cutoff at 10Hz sampling + float expAlpha = 0.816f; + float expAvg = 0.0f; + for (int loop = 0; loop < 4; loop++) { + expAvg = signal.front(); + for (int idx = 0; idx < length; idx++) { + expAvg = (expAlpha * signal.at(idx)) + ((1 - expAlpha) * expAvg); + signal[idx] = expAvg; + } + } + expAlpha = 0.268f; + for (int loop = 0; loop < 4; loop++) { + expAvg = signal.front(); + for (int idx = 0; idx < length; idx++) { + expAvg = (expAlpha * signal.at(idx)) + ((1 - expAlpha) * expAvg); + signal[idx] -= expAvg; } - z2 = z1; - z1 = z; } - return -1; } -} -Ppg::Ppg() - : hpf {0.87033078, -1.74066156, 0.87033078, -1.72377617, 0.75754694}, - agc {20, 0.971, 2}, - lpf {0.11595249, 0.23190498, 0.11595249, -0.72168143, 0.18549138} { -} + float SpectrumMax(const std::array<float, Ppg::spectrumLength>& data, int start, int end) { + float max = 0.0f; + for (int idx = start; idx < end; idx++) { + if (data.at(idx) > max) { + max = data.at(idx); + } + } + return max; + } -int8_t Ppg::Preprocess(float spl) { - spl -= offset; - spl = hpf.Step(spl); - spl = agc.Step(spl); - spl = lpf.Step(spl); + void Detrend(std::array<float, Ppg::dataLength>& signal) { + int size = signal.size(); + float offset = signal.front(); + float slope = (signal.at(size - 1) - offset) / static_cast<float>(size - 1); - auto spl_int = static_cast<int8_t>(spl); + for (int idx = 0; idx < size; idx++) { + signal[idx] -= (slope * static_cast<float>(idx) + offset); + } + for (int idx = 0; idx < size - 1; idx++) { + signal[idx] = signal[idx + 1] - signal[idx]; + } + } - if (dataIndex < 200) { - data[dataIndex++] = spl_int; + // Hanning Coefficients from numpy: python -c 'import numpy;print(numpy.hanning(64))' + // Note: Harcoded and must be updated if constexpr dataLength is changed. Prevents the need to + // use cosf() which results in an extra ~5KB in storage. + // This data is symetrical so just using the first half (saves 128B when dataLength is 64). + static constexpr float hanning[Ppg::dataLength >> 1] { + 0.0f, 0.00248461f, 0.00991376f, 0.0222136f, 0.03926189f, 0.06088921f, 0.08688061f, 0.11697778f, + 0.15088159f, 0.1882551f, 0.22872687f, 0.27189467f, 0.31732949f, 0.36457977f, 0.41317591f, 0.46263495f, + 0.51246535f, 0.56217185f, 0.61126047f, 0.65924333f, 0.70564355f, 0.75f, 0.79187184f, 0.83084292f, + 0.86652594f, 0.89856625f, 0.92664544f, 0.95048443f, 0.96984631f, 0.98453864f, 0.99441541f, 0.99937846f}; +} + +Ppg::Ppg() { + dataAverage.fill(0.0f); + spectrum.fill(0.0f); +} + +int8_t Ppg::Preprocess(uint32_t hrs, uint32_t als) { + if (dataIndex < dataLength) { + dataHRS[dataIndex++] = hrs; } - return spl_int; + alsValue = als; + if (alsValue > alsThreshold) { + return 1; + } + return 0; } int Ppg::HeartRate() { - if (dataIndex < 200) { + if (dataIndex < dataLength) { return 0; } - - NRF_LOG_INFO("PREPROCESS, offset = %d", offset); - auto hr = ProcessHeartRate(); - dataIndex = 0; + int hr = 0; + hr = ProcessHeartRate(resetSpectralAvg); + resetSpectralAvg = false; + // Make room for overlapWindow number of new samples + for (int idx = 0; idx < dataLength - overlapWindow; idx++) { + dataHRS[idx] = dataHRS[idx + overlapWindow]; + } + dataIndex = dataLength - overlapWindow; return hr; } -int Ppg::ProcessHeartRate() { - int t0 = Trough(data.data(), dataIndex, 7, 48); - if (t0 < 0) { - return 0; +void Ppg::Reset(bool resetDaqBuffer) { + if (resetDaqBuffer) { + dataIndex = 0; } + avgIndex = 0; + dataAverage.fill(0.0f); + lastPeakLocation = 0.0f; + alsThreshold = UINT16_MAX; + alsValue = 0; + resetSpectralAvg = true; + spectrum.fill(0.0f); +} - int t1 = t0 * 2; - t1 = Trough(data.data(), dataIndex, t1 - 5, t1 + 5); - if (t1 < 0) { - return 0; +// Pass init == true to reset spectral averaging. +// Returns -1 (Reset Acquisition), 0 (Unable to obtain HR) or HR (BPM). +int Ppg::ProcessHeartRate(bool init) { + std::copy(dataHRS.begin(), dataHRS.end(), vReal.begin()); + Detrend(vReal); + Filter30to240(vReal); + vImag.fill(0.0f); + // Apply Hanning Window + int hannIdx = 0; + for (int idx = 0; idx < dataLength; idx++) { + if (idx >= dataLength >> 1) { + hannIdx--; + } + vReal[idx] *= hanning[hannIdx]; + if (idx < dataLength >> 1) { + hannIdx++; + } } - - int t2 = (t1 * 3) / 2; - t2 = Trough(data.data(), dataIndex, t2 - 5, t2 + 5); - if (t2 < 0) { - return 0; + // Compute in place power spectrum + ArduinoFFT<float> FFT = ArduinoFFT<float>(vReal.data(), vImag.data(), dataLength, sampleFreq); + FFT.compute(FFTDirection::Forward); + FFT.complexToMagnitude(); + FFT.~ArduinoFFT(); + SpectrumAverage(vReal.data(), spectrum.data(), spectrum.size(), init); + peakLocation = 0.0f; + float threshold = peakDetectionThreshold; + float peakWidth = 0.0f; + int specLen = spectrum.size(); + float max = SpectrumMax(spectrum, hrROIbegin, hrROIend); + float signalToNoiseRatio = SignalToNoise(spectrum, hrROIbegin, hrROIend, max); + if (signalToNoiseRatio > signalToNoiseThreshold && spectrum.at(0) < dcThreshold) { + threshold *= max; + // Reuse VImag for interpolation x values passed to PeakSearch + for (int idx = 0; idx < dataLength; idx++) { + vImag[idx] = idx; + } + peakLocation = PeakSearch(vImag.data(), + spectrum.data(), + threshold, + peakWidth, + static_cast<float>(hrROIbegin), + static_cast<float>(hrROIend), + specLen); + peakLocation *= freqResolution; } - - int t3 = (t2 * 4) / 3; - t3 = Trough(data.data(), dataIndex, t3 - 4, t3 + 4); - if (t3 < 0) { - return (60 * 24 * 3) / t2; + // Peak too wide? (broad spectrum noise or large, rapid HR change) + if (peakWidth > maxPeakWidth) { + peakLocation = 0.0f; } - - return (60 * 24 * 4) / t3; + // Check HR limits + if (peakLocation < minHR || peakLocation > maxHR) { + peakLocation = 0.0f; + } + // Reset spectral averaging if bad reading + if (peakLocation == 0.0f) { + resetSpectralAvg = true; + } + // Set the ambient light threshold and return HR in BPM + alsThreshold = static_cast<uint16_t>(alsValue * alsFactor); + // Get current average HR. If HR reduced to zero, return -1 (reset) else HR + peakLocation = HeartRateAverage(peakLocation); + int rtn = -1; + if (peakLocation == 0.0f && lastPeakLocation > 0.0f) { + lastPeakLocation = 0.0f; + } else { + lastPeakLocation = peakLocation; + rtn = static_cast<int>((peakLocation * 60.0f) + 0.5f); + } + return rtn; } -void Ppg::SetOffset(uint16_t offset) { - this->offset = offset; - dataIndex = 0; +void Ppg::SpectrumAverage(const float* data, float* spectrum, int length, bool reset) { + if (reset) { + spectralAvgCount = 0; + } + float count = static_cast<float>(spectralAvgCount); + for (int idx = 0; idx < length; idx++) { + spectrum[idx] = (spectrum[idx] * count + data[idx]) / (count + 1); + } + if (spectralAvgCount < spectralAvgMax) { + spectralAvgCount++; + } } -void Ppg::Reset() { - dataIndex = 0; +float Ppg::HeartRateAverage(float hr) { + avgIndex++; + avgIndex %= dataAverage.size(); + dataAverage[avgIndex] = hr; + float avg = 0.0f; + float total = 0.0f; + float min = 300.0f; + float max = 0.0f; + for (const float& value : dataAverage) { + if (value > 0.0f) { + avg += value; + if (value < min) + min = value; + if (value > max) + max = value; + total++; + } + } + if (total > 0) { + avg /= total; + } else { + avg = 0.0f; + } + return avg; } diff --git a/src/components/heartrate/Ppg.h b/src/components/heartrate/Ppg.h index 1f709bab..2f8a1faa 100644 --- a/src/components/heartrate/Ppg.h +++ b/src/components/heartrate/Ppg.h @@ -3,29 +3,77 @@ #include <array> #include <cstddef> #include <cstdint> -#include "components/heartrate/Biquad.h" -#include "components/heartrate/Ptagc.h" +// Note: Change internal define 'sqrt_internal sqrt' to +// 'sqrt_internal sqrtf' to save ~3KB of flash. +#define FFT_SPEED_OVER_PRECISION +#include "libs/arduinoFFT-develop/src/arduinoFFT.h" namespace Pinetime { namespace Controllers { class Ppg { public: Ppg(); - int8_t Preprocess(float spl); + int8_t Preprocess(uint32_t hrs, uint32_t als); int HeartRate(); - - void SetOffset(uint16_t offset); - void Reset(); + void Reset(bool resetDaqBuffer); + static constexpr int deltaTms = 100; + // Daq dataLength: Must be power of 2 + static constexpr uint16_t dataLength = 64; + static constexpr uint16_t spectrumLength = dataLength >> 1; private: - std::array<int8_t, 200> data; - size_t dataIndex = 0; - float offset; - Biquad hpf; - Ptagc agc; - Biquad lpf; + // The sampling frequency (Hz) based on sampling time in milliseconds (DeltaTms) + static constexpr float sampleFreq = 1000.0f / static_cast<float>(deltaTms); + // The frequency resolution (Hz) + static constexpr float freqResolution = sampleFreq / dataLength; + // Number of samples before each analysis + // 0.5 second update rate at 10Hz + static constexpr uint16_t overlapWindow = 5; + // Maximum number of spectrum running averages + // Note: actual number of spectra averaged = spectralAvgMax + 1 + static constexpr uint16_t spectralAvgMax = 2; + // Multiple Peaks above this threshold (% of max) are rejected + static constexpr float peakDetectionThreshold = 0.6f; + // Maximum peak width (bins) at threshold for valid peak. + static constexpr float maxPeakWidth = 2.5f; + // Metric for spectrum noise level. + static constexpr float signalToNoiseThreshold = 3.0f; + // Heart rate Region Of Interest begin (bins) + static constexpr uint16_t hrROIbegin = static_cast<uint16_t>((30.0f / 60.0f) / freqResolution + 0.5f); + // Heart rate Region Of Interest end (bins) + static constexpr uint16_t hrROIend = static_cast<uint16_t>((240.0f / 60.0f) / freqResolution + 0.5f); + // Minimum HR (Hz) + static constexpr float minHR = 40.0f / 60.0f; + // Maximum HR (Hz) + static constexpr float maxHR = 230.0f / 60.0f; + // Threshold for high DC level after filtering + static constexpr float dcThreshold = 0.5f; + // ALS detection factor + static constexpr float alsFactor = 2.0f; + + // Raw ADC data + std::array<uint16_t, dataLength> dataHRS; + // Stores Real numbers from FFT + std::array<float, dataLength> vReal; + // Stores Imaginary numbers from FFT + std::array<float, dataLength> vImag; + // Stores power spectrum calculated from FFT real and imag values + std::array<float, (spectrumLength)> spectrum; + // Stores each new HR value (Hz). Non zero values are averaged for HR output + std::array<float, 20> dataAverage; + + uint16_t avgIndex = 0; + uint16_t spectralAvgCount = 0; + float lastPeakLocation = 0.0f; + uint16_t alsThreshold = UINT16_MAX; + uint16_t alsValue = 0; + uint16_t dataIndex = 0; + float peakLocation; + bool resetSpectralAvg = true; - int ProcessHeartRate(); + int ProcessHeartRate(bool init); + float HeartRateAverage(float hr); + void SpectrumAverage(const float* data, float* spectrum, int length, bool reset); }; } } diff --git a/src/components/heartrate/Ptagc.cpp b/src/components/heartrate/Ptagc.cpp deleted file mode 100644 index 221be460..00000000 --- a/src/components/heartrate/Ptagc.cpp +++ /dev/null @@ -1,29 +0,0 @@ -/* - SPDX-License-Identifier: LGPL-3.0-or-later - Original work Copyright (C) 2020 Daniel Thompson - C++ port Copyright (C) 2021 Jean-François Milants -*/ - -#include "components/heartrate/Ptagc.h" -#include <cmath> - -using namespace Pinetime::Controllers; - -/** Original implementation from wasp-os : https://github.com/daniel-thompson/wasp-os/blob/master/wasp/ppg.py */ -Ptagc::Ptagc(float start, float decay, float threshold) : peak {start}, decay {decay}, boost {1.0f / decay}, threshold {threshold} { -} - -float Ptagc::Step(float spl) { - if (std::abs(spl) > peak) { - peak *= boost; - } else { - peak *= decay; - } - - if ((spl > (peak * threshold)) || (spl < (peak * -threshold))) { - return 0.0f; - } - - spl = 100.0f * spl / (2.0f * peak); - return spl; -} diff --git a/src/components/heartrate/Ptagc.h b/src/components/heartrate/Ptagc.h deleted file mode 100644 index 3476636b..00000000 --- a/src/components/heartrate/Ptagc.h +++ /dev/null @@ -1,17 +0,0 @@ -#pragma once - -namespace Pinetime { - namespace Controllers { - class Ptagc { - public: - Ptagc(float start, float decay, float threshold); - float Step(float spl); - - private: - float peak; - float decay; - float boost; - float threshold; - }; - } -} diff --git a/src/drivers/Hrs3300.cpp b/src/drivers/Hrs3300.cpp index 6c47ae28..004c5362 100644 --- a/src/drivers/Hrs3300.cpp +++ b/src/drivers/Hrs3300.cpp @@ -16,6 +16,8 @@ using namespace Pinetime::Drivers; /** Driver for the HRS3300 heart rate sensor. * Original implementation from wasp-os : https://github.com/daniel-thompson/wasp-os/blob/master/wasp/drivers/hrs3300.py + * + * Experimentaly derived changes to improve signal/noise (see comments below) - Ceimour */ Hrs3300::Hrs3300(TwiMaster& twiMaster, uint8_t twiAddress) : twiMaster {twiMaster}, twiAddress {twiAddress} { } @@ -26,19 +28,21 @@ void Hrs3300::Init() { Disable(); vTaskDelay(100); - // HRS disabled, 12.5 ms wait time between cycles, (partly) 20mA drive - WriteRegister(static_cast<uint8_t>(Registers::Enable), 0x60); + // HRS disabled, 50ms wait time between ADC conversion period, current 12.5mA + WriteRegister(static_cast<uint8_t>(Registers::Enable), 0x50); - // (partly) 20mA drive, power on, "magic" (datasheet says both - // "reserved" and "set low nibble to 8" but 0xe gives better results - // and is used by at least two other HRS3300 drivers - WriteRegister(static_cast<uint8_t>(Registers::PDriver), 0x6E); + // Current 12.5mA and low nibble 0xF. + // Note: Setting low nibble to 0x8 per the datasheet results in + // modulated LED driver output. Setting to 0xF results in clean, + // steady output during the ADC conversion period. + WriteRegister(static_cast<uint8_t>(Registers::PDriver), 0x2f); - // HRS and ALS both in 16-bit mode - WriteRegister(static_cast<uint8_t>(Registers::Res), 0x88); + // HRS and ALS both in 15-bit mode results in ~50ms LED drive period + // and presumably ~50ms ADC conversion period. + WriteRegister(static_cast<uint8_t>(Registers::Res), 0x77); - // 8x gain, non default, reduced value for better readings - WriteRegister(static_cast<uint8_t>(Registers::Hgain), 0xc); + // Gain set to 1x + WriteRegister(static_cast<uint8_t>(Registers::Hgain), 0x00); } void Hrs3300::Enable() { diff --git a/src/heartratetask/HeartRateTask.cpp b/src/heartratetask/HeartRateTask.cpp index 50833ab2..414cdf2c 100644 --- a/src/heartratetask/HeartRateTask.cpp +++ b/src/heartratetask/HeartRateTask.cpp @@ -26,10 +26,11 @@ void HeartRateTask::Process(void* instance) { void HeartRateTask::Work() { int lastBpm = 0; while (true) { - auto delay = portMAX_DELAY; + Messages msg; + uint32_t delay; if (state == States::Running) { if (measurementStarted) { - delay = 40; + delay = ppg.deltaTms; } else { delay = 100; } @@ -37,8 +38,7 @@ void HeartRateTask::Work() { delay = portMAX_DELAY; } - Messages msg; - if (xQueueReceive(messageQueue, &msg, delay) == pdTRUE) { + if (xQueueReceive(messageQueue, &msg, delay)) { switch (msg) { case Messages::GoToSleep: StopMeasurement(); @@ -70,12 +70,28 @@ void HeartRateTask::Work() { } if (measurementStarted) { - ppg.Preprocess(static_cast<float>(heartRateSensor.ReadHrs())); - auto bpm = ppg.HeartRate(); + int8_t ambient = ppg.Preprocess(heartRateSensor.ReadHrs(), heartRateSensor.ReadAls()); + int bpm = ppg.HeartRate(); + + // If ambient light detected or a reset requested (bpm < 0) + if (ambient > 0) { + // Reset all DAQ buffers + ppg.Reset(true); + // Force state to NotEnoughData (below) + lastBpm = 0; + bpm = 0; + } else if (bpm < 0) { + // Reset all DAQ buffers except HRS buffer + ppg.Reset(false); + // Set HR to zero and update + bpm = 0; + controller.Update(Controllers::HeartRateController::States::Running, bpm); + } if (lastBpm == 0 && bpm == 0) { - controller.Update(Controllers::HeartRateController::States::NotEnoughData, 0); + controller.Update(Controllers::HeartRateController::States::NotEnoughData, bpm); } + if (bpm != 0) { lastBpm = bpm; controller.Update(Controllers::HeartRateController::States::Running, lastBpm); @@ -87,7 +103,7 @@ void HeartRateTask::Work() { void HeartRateTask::PushMessage(HeartRateTask::Messages msg) { BaseType_t xHigherPriorityTaskWoken = pdFALSE; xQueueSendFromISR(messageQueue, &msg, &xHigherPriorityTaskWoken); - if (xHigherPriorityTaskWoken == pdTRUE) { + if (xHigherPriorityTaskWoken) { /* Actual macro used here is port specific. */ // TODO : should I do something here? } @@ -95,11 +111,12 @@ void HeartRateTask::PushMessage(HeartRateTask::Messages msg) { void HeartRateTask::StartMeasurement() { heartRateSensor.Enable(); + ppg.Reset(true); vTaskDelay(100); - ppg.SetOffset(heartRateSensor.ReadHrs()); } void HeartRateTask::StopMeasurement() { heartRateSensor.Disable(); + ppg.Reset(true); vTaskDelay(100); } diff --git a/src/libs/arduinoFFT-develop/.gitignore b/src/libs/arduinoFFT-develop/.gitignore new file mode 100644 index 00000000..669c7704 --- /dev/null +++ b/src/libs/arduinoFFT-develop/.gitignore @@ -0,0 +1,3 @@ +/.project +/sync.ffs_db +*.*bak diff --git a/src/libs/arduinoFFT-develop/Examples/FFT_01/FFT_01.ino b/src/libs/arduinoFFT-develop/Examples/FFT_01/FFT_01.ino new file mode 100644 index 00000000..22b5024a --- /dev/null +++ b/src/libs/arduinoFFT-develop/Examples/FFT_01/FFT_01.ino @@ -0,0 +1,119 @@ +/* + + Example of use of the FFT libray + + Copyright (C) 2014 Enrique Condes + Copyright (C) 2020 Bim Overbohm (header-only, 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/>. + +*/ + +/* + In this example, the Arduino simulates the sampling of a sinusoidal 1000 Hz + signal with an amplitude of 100, sampled at 5000 Hz. Samples are stored + inside the vReal array. The samples are windowed according to Hamming + function. The FFT is computed using the windowed samples. Then the magnitudes + of each of the frequencies that compose the signal are calculated. Finally, + the frequency with the highest peak is obtained, being that the main frequency + present in the signal. +*/ + +#include "arduinoFFT.h" + +/* +These values can be changed in order to evaluate the functions +*/ +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 +*/ +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 +#define SCL_PLOT 0x03 + +void setup() +{ + Serial.begin(115200); + Serial.println("Ready"); +} + +void loop() +{ + /* Build raw data */ + double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read + for (uint16_t i = 0; i < samples; i++) + { + vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/ + //vReal[i] = uint8_t((amplitude * (sin((i * (twoPi * cycles)) / samples) + 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 + } + /* Print the results of the simulated 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); + double x = FFT.majorPeak(); + Serial.println(x, 6); + while(1); /* Run Once */ + // delay(2000); /* Repeat after delay */ +} + +void PrintVector(double *vData, uint16_t bufferSize, uint8_t scaleType) +{ + for (uint16_t i = 0; i < bufferSize; i++) + { + double 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(); +} diff --git a/src/libs/arduinoFFT-develop/Examples/FFT_02/FFT_02.ino b/src/libs/arduinoFFT-develop/Examples/FFT_02/FFT_02.ino new file mode 100644 index 00000000..7164dab1 --- /dev/null +++ b/src/libs/arduinoFFT-develop/Examples/FFT_02/FFT_02.ino @@ -0,0 +1,125 @@ +/* + + 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 (header-only, template, speed improvements) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + +*/ + +#include "arduinoFFT.h" + +/* +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; +const double startFrequency = 2; +const double stopFrequency = 16.4; +const double step_size = 0.1; + +/* +These are the input and output vectors +Input vectors receive computed results from FFT +*/ +double vReal[samples]; +double vImag[samples]; + +/* Create FFT object */ +ArduinoFFT<double> FFT = ArduinoFFT<double>(vReal, vImag, samples, sampling); + +unsigned long startTime; + +#define SCL_INDEX 0x00 +#define SCL_TIME 0x01 +#define SCL_FREQUENCY 0x02 +#define SCL_PLOT 0x03 + +void setup() +{ + Serial.begin(115200); + Serial.println("Ready"); +} + +void loop() +{ + Serial.println("Frequency\tDetected\ttakes (ms)"); + Serial.println("=======================================\n"); + for(double frequency = startFrequency; frequency<=stopFrequency; frequency+=step_size) + { + /* Build raw data */ + double cycles = (((samples-1) * frequency) / sampling); + for (uint16_t i = 0; i < samples; i++) + { + vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0); + vImag[i] = 0; //Reset the imaginary values vector for each new frequency + } + /*Serial.println("Data:"); + PrintVector(vReal, samples, SCL_TIME);*/ + startTime=millis(); + 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);*/ + double x = FFT.majorPeak(); + Serial.print(frequency); + Serial.print(": \t\t"); + Serial.print(x, 4); + Serial.print("\t\t"); + Serial.print(millis()-startTime); + Serial.println(" ms"); + // delay(2000); /* Repeat after delay */ + } + while(1); /* Run Once */ +} + +void PrintVector(double *vData, uint16_t bufferSize, uint8_t scaleType) +{ + for (uint16_t i = 0; i < bufferSize; i++) + { + double abscissa; + /* Print abscissa value */ + switch (scaleType) + { + case SCL_INDEX: + abscissa = (i * 1.0); + break; + case SCL_TIME: + abscissa = ((i * 1.0) / sampling); + break; + case SCL_FREQUENCY: + abscissa = ((i * 1.0 * sampling) / samples); + break; + } + Serial.print(abscissa, 6); + if(scaleType==SCL_FREQUENCY) + Serial.print("Hz"); + Serial.print(" "); + Serial.println(vData[i], 4); + } + Serial.println(); +} diff --git a/src/libs/arduinoFFT-develop/Examples/FFT_03/FFT_03.ino b/src/libs/arduinoFFT-develop/Examples/FFT_03/FFT_03.ino new file mode 100644 index 00000000..ee2b2946 --- /dev/null +++ b/src/libs/arduinoFFT-develop/Examples/FFT_03/FFT_03.ino @@ -0,0 +1,114 @@ +/* + + 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 (header-only, template, speed improvements) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + +*/ + +#include "arduinoFFT.h" + +/* +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; + +/* +These are the input and output vectors +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 +#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); + 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 */ +} + +void PrintVector(double *vData, uint16_t bufferSize, uint8_t scaleType) +{ + for (uint16_t i = 0; i < bufferSize; i++) + { + double 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(); +} diff --git a/src/libs/arduinoFFT-develop/Examples/FFT_04/FFT_04.ino b/src/libs/arduinoFFT-develop/Examples/FFT_04/FFT_04.ino new file mode 100644 index 00000000..b125991d --- /dev/null +++ b/src/libs/arduinoFFT-develop/Examples/FFT_04/FFT_04.ino @@ -0,0 +1,110 @@ +/* + + Example of use of the FFT libray + + Copyright (C) 2018 Enrique Condes + Copyright (C) 2020 Bim Overbohm (header-only, 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/>. + +*/ + +/* + In this example, the Arduino simulates the sampling of a sinusoidal 1000 Hz + signal with an amplitude of 100, sampled at 5000 Hz. Samples are stored + inside the vReal array. The samples are windowed according to Hamming + function. The FFT is computed using the windowed samples. Then the magnitudes + of each of the frequencies that compose the signal are calculated. Finally, + the frequency spectrum magnitudes are printed. If you use the Arduino IDE + serial plotter, you will see a single spike corresponding to the 1000 Hz + frecuency. +*/ + +#include "arduinoFFT.h" + +/* +These values can be changed in order to evaluate the functions +*/ +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 +*/ +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 +#define SCL_PLOT 0x03 + +void setup() +{ + Serial.begin(115200); +} + +void loop() +{ + /* Build raw data */ + double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read + for (uint16_t i = 0; i < samples; i++) + { + vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/ + //vReal[i] = uint8_t((amplitude * (sin((i * (twoPi * cycles)) / samples) + 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.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(); + while(1); /* Run Once */ + // delay(2000); /* Repeat after delay */ +} + +void PrintVector(double *vData, uint16_t bufferSize, uint8_t scaleType) +{ + for (uint16_t i = 0; i < bufferSize; i++) + { + double 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; + } + 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/libs/arduinoFFT-develop/Examples/FFT_05/FFT_05.ino b/src/libs/arduinoFFT-develop/Examples/FFT_05/FFT_05.ino new file mode 100644 index 00000000..a6f4df7a --- /dev/null +++ b/src/libs/arduinoFFT-develop/Examples/FFT_05/FFT_05.ino @@ -0,0 +1,124 @@ +/* + + Example of use of the FFT libray + + Copyright (C) 2014 Enrique Condes + Copyright (C) 2020 Bim Overbohm (header-only, 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/>. + +*/ + +/* + In this example, the Arduino simulates the sampling of a sinusoidal 1000 Hz + signal with an amplitude of 100, sampled at 5000 Hz. Samples are stored + inside the vReal array. The samples are windowed according to Hamming + function. The FFT is computed using the windowed samples. Then the magnitudes + of each of the frequencies that compose the signal are calculated. Finally, + the frequency with the highest peak is obtained, being that the main frequency + present in the signal. This frequency is printed, along with the magnitude of + the peak. +*/ + +#include "arduinoFFT.h" + +/* +These values can be changed in order to evaluate the functions +*/ +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 +*/ +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 +#define SCL_PLOT 0x03 + +void setup() +{ + Serial.begin(115200); + Serial.println("Ready"); +} + +void loop() +{ + /* Build raw data */ + double cycles = (((samples-1) * signalFrequency) / samplingFrequency); //Number of signal cycles that the sampling will read + for (uint16_t i = 0; i < samples; i++) + { + vReal[i] = int8_t((amplitude * (sin((i * (TWO_PI * cycles)) / samples))) / 2.0);/* Build data with positive and negative values*/ + //vReal[i] = uint8_t((amplitude * (sin((i * (twoPi * cycles)) / samples) + 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 + } + /* Print the results of the simulated 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); + double x; + double v; + FFT.majorPeak(x, v); + Serial.print(x, 6); + Serial.print(", "); + Serial.println(v, 6); + while(1); /* Run Once */ + // delay(2000); /* Repeat after delay */ +} + +void PrintVector(double *vData, uint16_t bufferSize, uint8_t scaleType) +{ + for (uint16_t i = 0; i < bufferSize; i++) + { + double 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(); +} diff --git a/src/libs/arduinoFFT-develop/Examples/FFT_speedup/FFT_speedup.ino b/src/libs/arduinoFFT-develop/Examples/FFT_speedup/FFT_speedup.ino new file mode 100644 index 00000000..a059a170 --- /dev/null +++ b/src/libs/arduinoFFT-develop/Examples/FFT_speedup/FFT_speedup.ino @@ -0,0 +1,129 @@ +/* + + 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 (header-only, 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]; + +/* +Allocate space for FFT window weighing factors, so they are calculated only the first time windowing() is called. +If you don't do this, a lot of calculations are necessary, depending on the window function. +*/ +float weighingFactors[samples]; + +/* Create FFT object with weighing factor storage */ +ArduinoFFT<float> FFT = ArduinoFFT<float>(vReal, vImag, samples, samplingFrequency, weighingFactors); + +#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(); +} diff --git a/src/libs/arduinoFFT-develop/LICENSE b/src/libs/arduinoFFT-develop/LICENSE new file mode 100644 index 00000000..70566f2d --- /dev/null +++ b/src/libs/arduinoFFT-develop/LICENSE @@ -0,0 +1,674 @@ +GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/> + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + {one line to give the program's name and a brief idea of what it does.} + Copyright (C) {year} {name of author} + + 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/>. + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + {project} Copyright (C) {year} {fullname} + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +<http://www.gnu.org/licenses/>. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +<http://www.gnu.org/philosophy/why-not-lgpl.html>.
\ No newline at end of file diff --git a/src/libs/arduinoFFT-develop/README.md b/src/libs/arduinoFFT-develop/README.md new file mode 100644 index 00000000..f9229ef9 --- /dev/null +++ b/src/libs/arduinoFFT-develop/README.md @@ -0,0 +1,129 @@ +arduinoFFT +========== + +# Fast Fourier Transform for Arduino + +This is a fork from https://code.google.com/p/makefurt/ which has been abandoned since 2011. +~~This is a C++ library for Arduino for computing FFT.~~ Now it works both on Arduino and C projects. This is version 2.0 of the library, which has a different [API](#api). See here [how to migrate from 1.x to 2.x](#migrating-from-1x-to-2x). +Tested on Arduino 1.6.11 and 1.8.10. + +## 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: + +`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) + +## 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>` + +## API + +* ```ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, T * weighingFactors = nullptr);``` +Constructor. +The type `T` can be `float` or `double`. `vReal` and `vImag` are pointers to arrays of real and imaginary data and have to be allocated outside of ArduinoFFT. `samples` is the number of samples in `vReal` and `vImag` and `weighingFactors` (if specified). `samplingFrequency` is the sample frequency of the data. `weighingFactors` can optionally be specified to cache weighing factors for the windowing function. This speeds up repeated calls to **windowing()** significantly. You can deallocate `vReal` and `vImag` after you are done using the library, or only use specific library functions that only need one of those arrays. + +```C++ +const uint32_t nrOfSamples = 1024; +auto real = new float[nrOfSamples]; +auto imag = new float[nrOfSamples]; +auto fft = ArduinoFFT<float>(real, imag, nrOfSamples, 10000); +// ... fill real + imag and use it ... +fft.compute(); +fft.complexToMagnitude(); +delete [] imag; +// ... continue using real and only functions that use real ... +auto peak = fft.majorPeak(); +``` +* ```~ArduinoFFT()``` +Destructor. +* ```void complexToMagnitude() const;``` +Convert complex values to their magnitude and store in vReal. Uses vReal and vImag. +* ```void compute(FFTDirection dir) const;``` +Calcuates the Fast Fourier Transform. Uses vReal and vImag. +* ```void dcRemoval() const;``` +Removes the DC component from the sample data. Uses vReal. +* ```T majorPeak() const;``` +Returns the frequency of the biggest spike in the analyzed signal. Uses vReal. +* ```void majorPeak(T &frequency, T &value) const;``` +Returns the frequency and the value of the biggest spike in the analyzed signal. Uses vReal. +* ```uint8_t revision() const;``` +Returns the library revision. +* ```void setArrays(T *vReal, T *vImag);``` +Replace the data array pointers. +* ```void windowing(FFTWindow windowType, FFTDirection dir, bool withCompensation = false);``` +Performs a windowing function on the values array. Uses vReal. The possible windowing options are: + * Rectangle + * Hamming + * Hann + * Triangle + * Nuttall + * Blackman + * Blackman_Nuttall + * Blackman_Harris + * Flat_top + * Welch + + If `withCompensation` == true, the following compensation factors are used: + * Rectangle: 1.0 * 2.0 + * Hamming: 1.8549343278 * 2.0 + * Hann: 1.8554726898 * 2.0 + * Triangle: 2.0039186079 * 2.0 + * Nuttall: 2.8163172034 * 2.0 + * Blackman: 2.3673474360 * 2.0 + * Blackman Nuttall: 2.7557840395 * 2.0 + * Blackman Harris: 2.7929062517 * 2.0 + * Flat top: 3.5659039231 * 2.0 + * Welch: 1.5029392863 * 2.0 + +## Special flags + +You can define these before including arduinoFFT.h: + +* #define FFT_SPEED_OVER_PRECISION +Define this to use reciprocal multiplication for division and some more speedups that might decrease precision. + +* #define FFT_SQRT_APPROXIMATION +Define this to use a low-precision square root approximation instead of the regular sqrt() call. This might only work for specific use cases, but is significantly faster. Only works if `T == float`. + +See the `FFT_speedup.ino` example in `Examples/FFT_speedup/FFT_speedup.ino`. + +# Migrating from 1.x to 2.x + +* The function signatures where you could pass in pointers were deprecated and have been removed. Pass in pointers to your real / imaginary array in the ArduinoFFT() constructor. If you have the need to replace those pointers during usage of the library (e.g. to free memory) you can do the following: + +```C++ +const uint32_t nrOfSamples = 1024; +auto real = new float[nrOfSamples]; +auto imag = new float[nrOfSamples]; +auto fft = ArduinoFFT<float>(real, imag, nrOfSamples, 10000); +// ... fill real + imag and use it ... +fft.compute(); +fft.complexToMagnitude(); +delete [] real; +// ... replace vReal in library with imag ... +fft.setArrays(imag, nullptr); +// ... keep doing whatever ... +``` +* All function names are camelCase case now (start with lower-case character), e.g. "windowing()" instead of "Windowing()". + +## TODO +* Ratio table for windowing function. +* Document windowing functions advantages and disadvantages. +* Optimize usage and arguments. +* Add new windowing functions. +* ~~Spectrum table?~~ diff --git a/src/libs/arduinoFFT-develop/changeLog.txt b/src/libs/arduinoFFT-develop/changeLog.txt new file mode 100644 index 00000000..d49b8548 --- /dev/null +++ b/src/libs/arduinoFFT-develop/changeLog.txt @@ -0,0 +1,40 @@ +02/22/20 v1.9.2 +Fix compilation on AVR systems. + +02/22/20 v1.9.1 +Add setArrays() function because of issue #32. +Add API migration info to README and improve README. +Use better sqrtf() approximation. + +02/19/20 v1.9.0 +Remove deprecated API. Consistent renaming of functions to lowercase. +Make template to be able to use float or double type (float brings a ~70% speed increase on ESP32). +Add option to provide cache for window function weighing factors (~50% speed increase on ESP32). +Add some #defines to enable math approximisations to further speed up code (~40% speed increase on ESP32). + +01/27/20 v1.5.5 +Lookup table for constants c1 and c2 used during FFT comupting. This increases the FFT computing speed in around 5%. + +02/10/18 v1.4 +Transition version. Minor optimization to functions. New API. Deprecation of old functions. + +12/06/18 v1.3 +Add support for mbed development boards. + +09/04/17 v1.2.3 +Finally solves the issue of Arduino IDE not correctly detecting and highlighting the keywords. + +09/03/17 v1.2.2 +Solves a format issue in keywords.txt that prevented keywords from being detected. + +08/28/17 v1.2.1 +Fix to issues 6 and 7. Not cleaning the imaginary vector after each cycle leaded to erroneous calculations and could cause buffer overflows. + +08/04/17 v1.2 +Fix to bug preventing the number of samples to be greater than 128. New logical limit is 32768 samples but it is bound to the RAM on the chip. + +05/12/17 v1.1 +Fix issue that prevented installation through the Arduino Library Manager interface. + +05/11/17 v1.0 +Initial commit to Arduino Library Manager. diff --git a/src/libs/arduinoFFT-develop/keywords.txt b/src/libs/arduinoFFT-develop/keywords.txt new file mode 100644 index 00000000..3807cdbd --- /dev/null +++ b/src/libs/arduinoFFT-develop/keywords.txt @@ -0,0 +1,41 @@ +####################################### +# Syntax Coloring Map For arduinoFFT +####################################### + +####################################### +# Datatypes (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 +setArrays KEYWORD2 + +####################################### +# Constants (LITERAL1) +####################################### + +Forward LITERAL1 +Reverse LITERAL1 +Rectangle LITERAL1 +Hamming LITERAL1 +Hann LITERAL1 +Triangle LITERAL1 +Nuttall LITERAL1 +Blackman LITERAL1 +Blackman_Nuttall LITERAL1 +Blackman_Harris LITERAL1 +Flat_top LITERAL1 +Welch LITERAL1 diff --git a/src/libs/arduinoFFT-develop/library.json b/src/libs/arduinoFFT-develop/library.json new file mode 100644 index 00000000..6c354193 --- /dev/null +++ b/src/libs/arduinoFFT-develop/library.json @@ -0,0 +1,31 @@ +{ + "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" + }, + { + "name": "Bim Overbohm", + "url": "https://github.com/HorstBaerbel", + "email": "bim.overbohm@googlemail.com" + } + ], + "version": "1.9.2", + "frameworks": ["arduino","mbed","espidf"], + "platforms": "*" +} diff --git a/src/libs/arduinoFFT-develop/library.properties b/src/libs/arduinoFFT-develop/library.properties new file mode 100644 index 00000000..0a909477 --- /dev/null +++ b/src/libs/arduinoFFT-develop/library.properties @@ -0,0 +1,10 @@ +name=arduinoFFT +version=1.9.2 +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. +paragraph=With this library you can calculate the frequency of a sampled signal. +category=Data Processing +url=https://github.com/kosme/arduinoFFT +architectures=* +includes=arduinoFFT.h diff --git a/src/libs/arduinoFFT-develop/src/.gitignore b/src/libs/arduinoFFT-develop/src/.gitignore new file mode 100644 index 00000000..00e95bf6 --- /dev/null +++ b/src/libs/arduinoFFT-develop/src/.gitignore @@ -0,0 +1 @@ +/arduinoFFT.h.gch diff --git a/src/libs/arduinoFFT-develop/src/arduinoFFT.h b/src/libs/arduinoFFT-develop/src/arduinoFFT.h new file mode 100644 index 00000000..fe8f9d91 --- /dev/null +++ b/src/libs/arduinoFFT-develop/src/arduinoFFT.h @@ -0,0 +1,498 @@ +/* + + FFT library + Copyright (C) 2010 Didier Longueville + Copyright (C) 2014 Enrique Condes + Copyright (C) 2020 Bim Overbohm (header-only, 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/>. + +*/ + +#ifndef ArduinoFFT_h /* Prevent loading library twice */ +#define ArduinoFFT_h +#ifdef ARDUINO +#if ARDUINO >= 100 +#include "Arduino.h" +#else +#include "WProgram.h" /* This is where the standard Arduino code lies */ +#endif +#else +#include <stdlib.h> +#include <stdio.h> +#ifdef __AVR__ +#include <avr/io.h> +#include <avr/pgmspace.h> +#endif +#include <math.h> +#include "defs.h" +#include "types.h" +#endif + +// 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 + +#ifdef FFT_SQRT_APPROXIMATION + #include <type_traits> +#else + #define sqrt_internal sqrtf +#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 +}; + +template <typename T> +class ArduinoFFT +{ +public: + // Constructor + ArduinoFFT(T *vReal, T *vImag, uint_fast16_t samples, T samplingFrequency, T *windowWeighingFactors = nullptr) + : _vReal(vReal) + , _vImag(vImag) + , _samples(samples) +#ifdef FFT_SPEED_OVER_PRECISION + , _oneOverSamples(1.0 / samples) +#endif + , _samplingFrequency(samplingFrequency) + , _windowWeighingFactors(windowWeighingFactors) + { + // Calculates the base 2 logarithm of sample count + _power = 0; + while (((samples >> _power) & 1) != 1) + { + _power++; + } + } + + // Destructor + ~ArduinoFFT() + { + } + + // Get library revision + static uint8_t revision() + { + return 0x19; + } + + // Replace the data array pointers + void setArrays(T *vReal, T *vImag) + { + _vReal = vReal; + _vImag = vImag; + } + + // Computes in-place complex-to-complex FFT + void compute(FFTDirection dir) const + { + // Reverse bits / + uint_fast16_t j = 0; + for (uint_fast16_t i = 0; i < (this->_samples - 1); i++) + { + if (i < j) + { + Swap(this->_vReal[i], this->_vReal[j]); + if (dir == FFTDirection::Reverse) + { + Swap(this->_vImag[i], this->_vImag[j]); + } + } + uint_fast16_t k = (this->_samples >> 1); + while (k <= j) + { + j -= k; + k >>= 1; + } + j += k; + } + // Compute the FFT +#ifdef __AVR__ + uint_fast8_t index = 0; +#endif + T c1 = -1.0; + T c2 = 0.0; + uint_fast16_t l2 = 1; + for (uint_fast8_t l = 0; (l < this->_power); l++) + { + uint_fast16_t l1 = l2; + l2 <<= 1; + T u1 = 1.0; + T u2 = 0.0; + for (j = 0; j < l1; j++) + { + for (uint_fast16_t i = j; i < this->_samples; i += l2) + { + uint_fast16_t i1 = i + l1; + T t1 = u1 * this->_vReal[i1] - u2 * this->_vImag[i1]; + T t2 = u1 * this->_vImag[i1] + u2 * this->_vReal[i1]; + this->_vReal[i1] = this->_vReal[i] - t1; + this->_vImag[i1] = this->_vImag[i] - t2; + this->_vReal[i] += t1; + this->_vImag[i] += t2; + } + T z = ((u1 * c1) - (u2 * c2)); + u2 = ((u1 * c2) + (u2 * c1)); + u1 = z; + } +#ifdef __AVR__ + c2 = pgm_read_float_near(&(_c2[index])); + c1 = pgm_read_float_near(&(_c1[index])); + index++; +#else + T cTemp = 0.5 * c1; + c2 = sqrt_internal(0.5 - cTemp); + c1 = sqrt_internal(0.5 + cTemp); +#endif + c2 = dir == FFTDirection::Forward ? -c2 : c2; + } + // Scaling for reverse transform + if (dir != FFTDirection::Forward) + { + for (uint_fast16_t i = 0; i < this->_samples; i++) + { +#ifdef FFT_SPEED_OVER_PRECISION + this->_vReal[i] *= _oneOverSamples; + this->_vImag[i] *= _oneOverSamples; +#else + this->_vReal[i] /= this->_samples; + this->_vImag[i] /= this->_samples; +#endif + } + } + } + + void complexToMagnitude() const + { + // vM is half the size of vReal and vImag + for (uint_fast16_t i = 0; i < this->_samples; i++) + { + this->_vReal[i] = sqrt_internal(sq(this->_vReal[i]) + sq(this->_vImag[i])); + } + } + + void dcRemoval() const + { + // calculate the mean of vData + T mean = 0; + for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++) + { + mean += this->_vReal[i]; + } + mean /= this->_samples; + // Subtract the mean from vData + for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++) + { + this->_vReal[i] -= mean; + } + } + + void windowing(FFTWindow windowType, FFTDirection dir, bool withCompensation = false) + { + // check if values are already pre-computed for the correct window type and compensation + if (_windowWeighingFactors && _weighingFactorsComputed && + _weighingFactorsFFTWindow == windowType && + _weighingFactorsWithCompensation == withCompensation) + { + // yes. values are precomputed + if (dir == FFTDirection::Forward) + { + for (uint_fast16_t i = 0; i < (this->_samples >> 1); i++) + { + this->_vReal[i] *= _windowWeighingFactors[i]; + this->_vReal[this->_samples - (i + 1)] *= _windowWeighingFactors[i]; + } + } + else + { + for (uint_fast16_t i = 0; i < (this->_samples >> 1); i++) + { +#ifdef FFT_SPEED_OVER_PRECISION + // on many architectures reciprocals and multiplying are much faster than division + T oneOverFactor = 1.0 / _windowWeighingFactors[i]; + this->_vReal[i] *= oneOverFactor; + this->_vReal[this->_samples - (i + 1)] *= oneOverFactor; +#else + this->_vReal[i] /= _windowWeighingFactors[i]; + this->_vReal[this->_samples - (i + 1)] /= _windowWeighingFactors[i]; +#endif + } + } + } + else + { + // no. values need to be pre-computed or applied + T samplesMinusOne = (T(this->_samples) - 1.0); + T compensationFactor = _WindowCompensationFactors[static_cast<uint_fast8_t>(windowType)]; + for (uint_fast16_t i = 0; i < (this->_samples >> 1); i++) + { + T indexMinusOne = T(i); + T ratio = (indexMinusOne / samplesMinusOne); + T weighingFactor = 1.0; + // Compute and record weighting factor + switch (windowType) + { + case FFTWindow::Rectangle: // rectangle (box car) + weighingFactor = 1.0; + break; + case FFTWindow::Hamming: // hamming + weighingFactor = 0.54 - (0.46 * cos(TWO_PI * ratio)); + break; + case FFTWindow::Hann: // hann + weighingFactor = 0.54 * (1.0 - cos(TWO_PI * ratio)); + break; + case FFTWindow::Triangle: // triangle (Bartlett) + weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne); + break; + case FFTWindow::Nuttall: // nuttall + weighingFactor = 0.355768 - (0.487396 * (cos(TWO_PI * ratio))) + (0.144232 * (cos(FOUR_PI * ratio))) - (0.012604 * (cos(SIX_PI * ratio))); + break; + case FFTWindow::Blackman: // blackman + weighingFactor = 0.42323 - (0.49755 * (cos(TWO_PI * ratio))) + (0.07922 * (cos(FOUR_PI * ratio))); + break; + case FFTWindow::Blackman_Nuttall: // blackman nuttall + weighingFactor = 0.3635819 - (0.4891775 * (cos(TWO_PI * ratio))) + (0.1365995 * (cos(FOUR_PI * ratio))) - (0.0106411 * (cos(SIX_PI * ratio))); + break; + case FFTWindow::Blackman_Harris: // blackman harris + weighingFactor = 0.35875 - (0.48829 * (cos(TWO_PI * ratio))) + (0.14128 * (cos(FOUR_PI * ratio))) - (0.01168 * (cos(SIX_PI * ratio))); + break; + case FFTWindow::Flat_top: // flat top + weighingFactor = 0.2810639 - (0.5208972 * cos(TWO_PI * ratio)) + (0.1980399 * cos(FOUR_PI * ratio)); + break; + case FFTWindow::Welch: // welch + weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0)); + break; + } + if (withCompensation) + { + weighingFactor *= compensationFactor; + } + if (_windowWeighingFactors) + { + _windowWeighingFactors[i] = weighingFactor; + } + if (dir == FFTDirection::Forward) + { + this->_vReal[i] *= weighingFactor; + this->_vReal[this->_samples - (i + 1)] *= weighingFactor; + } + else + { +#ifdef FFT_SPEED_OVER_PRECISION + // on many architectures reciprocals and multiplying are much faster than division + T oneOverFactor = 1.0 / weighingFactor; + this->_vReal[i] *= oneOverFactor; + this->_vReal[this->_samples - (i + 1)] *= oneOverFactor; +#else + this->_vReal[i] /= weighingFactor; + this->_vReal[this->_samples - (i + 1)] /= weighingFactor; +#endif + } + } + // mark cached values as pre-computed + _weighingFactorsFFTWindow = windowType; + _weighingFactorsWithCompensation = withCompensation; + _weighingFactorsComputed = true; + } + } + + T majorPeak() const + { + T maxY = 0; + uint_fast16_t IndexOfMaxY = 0; + //If sampling_frequency = 2 * max_frequency in signal, + //value would be stored at position samples/2 + for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++) + { + if ((this->_vReal[i - 1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i + 1])) + { + if (this->_vReal[i] > maxY) + { + maxY = this->_vReal[i]; + IndexOfMaxY = i; + } + } + } + T delta = 0.5 * ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) / (this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1])); + T interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1); + if (IndexOfMaxY == (this->_samples >> 1)) + { + //To improve calculation on edge values + interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples); + } + // returned value: interpolated frequency peak apex + return interpolatedX; + } + + void majorPeak(T &frequency, T &value) const + { + T maxY = 0; + uint_fast16_t IndexOfMaxY = 0; + //If sampling_frequency = 2 * max_frequency in signal, + //value would be stored at position samples/2 + for (uint_fast16_t i = 1; i < ((this->_samples >> 1) + 1); i++) + { + if ((this->_vReal[i - 1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i + 1])) + { + if (this->_vReal[i] > maxY) + { + maxY = this->_vReal[i]; + IndexOfMaxY = i; + } + } + } + T delta = 0.5 * ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) / (this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1])); + T interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples - 1); + if (IndexOfMaxY == (this->_samples >> 1)) + { + //To improve calculation on edge values + interpolatedX = ((IndexOfMaxY + delta) * this->_samplingFrequency) / (this->_samples); + } + // returned value: interpolated frequency peak apex + frequency = interpolatedX; + value = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]); + } + +private: +#ifdef __AVR__ + static const float _c1[] PROGMEM; + static const float _c2[] PROGMEM; +#endif + static const T _WindowCompensationFactors[10]; + + // Mathematial constants +#ifndef TWO_PI + static constexpr T TWO_PI = 6.28318531; // might already be defined in Arduino.h +#endif + static constexpr T FOUR_PI = 12.56637061; + static constexpr T SIX_PI = 18.84955593; + + static inline void Swap(T &x, T &y) + { + T 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 V = T> + static inline V sqrt_internal(typename std::enable_if<std::is_same<V, float>::value, V>::type x) + { + union // get bits for float value + { + float x; + int32_t i; + } u; + u.x = x; + u.i = 0x5f375a86 - (u.i >> 1); // gives initial guess y0. + float xu = x * u.x; + float xu2 = xu * u.x; + u.x = (0.125 * 3.0) * xu * (5.0 - xu2 * ((10.0 / 3.0) - xu2)); // Halley's method, repeating increases accuracy + return u.x; + } + + template <typename V = T> + static inline V sqrt_internal(typename std::enable_if<std::is_same<V, double>::value, V>::type x) + { + // According to HosrtBaerbel, on the ESP32 the approximation is not faster, so we use the standard function + #ifdef ESP32 + return sqrt(x); + #else + union // get bits for float value + { + double x; + int64_t i; + } u; + u.x = x; + u.i = 0x5fe6ec85e7de30da - (u.i >> 1); // gives initial guess y0. + double xu = x * u.x; + double xu2 = xu * u.x; + u.x = (0.125 * 3.0) * xu * (5.0 - xu2 * ((10.0 / 3.0) - xu2)); // Halley's method, repeating increases accuracy + return u.x; + #endif + } +#endif + + /* Variables */ + T *_vReal = nullptr; + T *_vImag = nullptr; + uint_fast16_t _samples = 0; +#ifdef FFT_SPEED_OVER_PRECISION + T _oneOverSamples = 0.0; +#endif + T _samplingFrequency = 0; + T *_windowWeighingFactors = nullptr; + FFTWindow _weighingFactorsFFTWindow; + bool _weighingFactorsWithCompensation = false; + bool _weighingFactorsComputed = false; + uint_fast8_t _power = 0; +}; + +#ifdef __AVR__ +template <typename T> +const float ArduinoFFT<T>::_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}; + +template <typename T> +const float ArduinoFFT<T>::_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 + +template <typename T> +const T ArduinoFFT<T>::_WindowCompensationFactors[10] = { + 1.0000000000 * 2.0, // rectangle (Box car) + 1.8549343278 * 2.0, // hamming + 1.8554726898 * 2.0, // hann + 2.0039186079 * 2.0, // triangle (Bartlett) + 2.8163172034 * 2.0, // nuttall + 2.3673474360 * 2.0, // blackman + 2.7557840395 * 2.0, // blackman nuttall + 2.7929062517 * 2.0, // blackman harris + 3.5659039231 * 2.0, // flat top + 1.5029392863 * 2.0 // welch +}; + +#endif diff --git a/src/libs/arduinoFFT-develop/src/defs.h b/src/libs/arduinoFFT-develop/src/defs.h new file mode 100644 index 00000000..2422b243 --- /dev/null +++ b/src/libs/arduinoFFT-develop/src/defs.h @@ -0,0 +1,90 @@ +/*! \file avrlibdefs.h \brief AVRlib global defines and macros. */ +//***************************************************************************** +// +// File Name : 'avrlibdefs.h' +// Title : AVRlib global defines and macros include file +// Author : Pascal Stang +// Created : 7/12/2001 +// Revised : 9/30/2002 +// Version : 1.1 +// Target MCU : Atmel AVR series +// Editor Tabs : 4 +// +// Description : This include file is designed to contain items useful to all +// code files and projects, regardless of specific implementation. +// +// This code is distributed under the GNU Public License +// which can be found at http://www.gnu.org/licenses/gpl.txt +// +//***************************************************************************** + + +#ifndef AVRLIBDEFS_H +#define AVRLIBDEFS_H + +//#define F_CPU 4000000 +#define MEM_TYPE 1 + +// Code compatibility to new AVR-libc +// outb(), inb(), inw(), outw(), BV(), sbi(), cbi(), sei(), cli() +#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 cbi +// #define cbi(reg,bit) reg &= ~(BV(bit)) +//#endif +//#ifndef sbi +// #define sbi(reg,bit) reg |= (BV(bit)) +//#endif +#ifndef cli + #define cli() __asm__ __volatile__ ("cli" ::) +#endif +#ifndef sei + #define sei() __asm__ __volatile__ ("sei" ::) +#endif + +// support for individual port pin naming in the mega128 +// see port128.h for details +#ifdef __AVR_ATmega128__ +// not currently necessary due to inclusion +// of these defines in newest AVR-GCC +// do a quick test to see if include is needed +#ifndef PD0 + //#include "port128.h" +#endif +#endif + +// use this for packed structures +// (this is seldom necessary on an 8-bit architecture like AVR, +// but can assist in code portability to AVR) +#define GNUC_PACKED __attribute__((packed)) + +// port address helpers +#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 + +// MIN/MAX/ABS macros +#define MIN(a,b) ((a<b)?(a):(b)) +#define MAX(a,b) ((a>b)?(a):(b)) +#define ABS(x) ((x>0)?(x):(-x)) + +// constants +#define PI 3.14159265359 + +//Math +#define sq(x) ((x)*(x)) +#define constrain(amt,low,high) ((amt)<(low)?(low):((amt)>(high)?(high):(amt))) + +#endif diff --git a/src/libs/arduinoFFT-develop/src/types.h b/src/libs/arduinoFFT-develop/src/types.h new file mode 100644 index 00000000..6cd7d854 --- /dev/null +++ b/src/libs/arduinoFFT-develop/src/types.h @@ -0,0 +1,69 @@ +//useful things to include in code + +#ifndef TYPES_H +#define TYPES_H + +#ifndef WIN32 + // true/false defines + #define FALSE 0 + #define TRUE -1 +#endif + +// datatype definitions macros +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; + +// #ifndef __AVR__ +#ifdef __MBED__ + // use inttypes.h instead + // C99 standard integer type definitions + typedef unsigned char uint8_t; + typedef signed char int8_t; + typedef unsigned short uint16_t; + typedef signed short int16_t; + /*typedef unsigned long uint32_t; + typedef signed long int32_t; + typedef unsigned long uint64_t; + typedef signed long int64_t; + */ +#endif + +// maximum value that can be held +// by unsigned data types (8,16,32bits) +#define MAX_U08 255 +#define MAX_U16 65535 +#define MAX_U32 4294967295 + +// maximum values that can be held +// by signed data types (8,16,32bits) +#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 + // more type redefinitions + 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 |