Skip to content

Commit

Permalink
update the Matrix filter, turn it into a dynamic notch filter (emufli…
Browse files Browse the repository at this point in the history
…ght#500)

* update the Matrix filter, turn it into a dynamic notch filter -- This filter is far more precise than the previous filter was. I've currently set this up with 36 bins (so accuracy is (max-min)/36 which for defaults is, 12.5hz, the old filter had a bin size of 16 which would have given an accuracy of 28.125hz, so this is over double as accurate.) This can also track up to 5 peaks and place 5 notch filters. It is also nicer on the cpu than the older version was, so all around really good. Using this filter + abg i was easily able to raise my gyro lpf filters to 300hz with dterm filters at around 150hz. I'm really excited to see what you all think about this filter :), P.S. all credit goes to @KarateBrot for coding this into betaflight, I may want to update this later to have dynamic widths and dynamic on the number of notches it uses, but so far, hey its awesome. Probably going to stay like this for 0.4.0. Only the defaults need testing right now.
* fix emugravity
* fix potential bug
* the approximation is slower and less accurate
* fixes
* update q correctly
* dterm dyn notch
* dterm dyn notch can be tuned
* Update gyro_filter_impl.h
* OSD and BB headers for dterm_dyn_notch_enable & dterm_dyn_notch_q (emuflight#638)
* OSD and BB headers for dterm_dyn_notch_enable & dterm_dyn_notch_q
* fix OSD toggle
* dynamicFilter/dynamicDtermNotch (emuflight#657)
* fix Codacy; set dynDtermNotch disabled by default
* update the Matrix filter, turn it into a dynamic notch filter
* fix Codacy again
* make codacy happy
Co-authored-by: nerdCopter <56646290+nerdCopter@users.noreply.github.com>
  • Loading branch information
Quick-Flash authored Aug 2, 2021
1 parent 9030ccc commit 75a10b7
Show file tree
Hide file tree
Showing 16 changed files with 561 additions and 321 deletions.
9 changes: 1 addition & 8 deletions make/source.mk
Original file line number Diff line number Diff line change
Expand Up @@ -190,6 +190,7 @@ SPEED_OPTIMISED_SRC := $(SPEED_OPTIMISED_SRC) \
common/encoding.c \
common/filter.c \
common/maths.c \
common/sdft.c \
common/typeconversion.c \
drivers/accgyro/accgyro_fake.c \
drivers/accgyro/accgyro_mpu.c \
Expand Down Expand Up @@ -335,14 +336,6 @@ ifneq ($(DSP_LIB),)

INCLUDE_DIRS += $(DSP_LIB)/Include

SRC += $(DSP_LIB)/Source/BasicMathFunctions/arm_mult_f32.c
SRC += $(DSP_LIB)/Source/TransformFunctions/arm_rfft_fast_f32.c
SRC += $(DSP_LIB)/Source/TransformFunctions/arm_cfft_f32.c
SRC += $(DSP_LIB)/Source/TransformFunctions/arm_rfft_fast_init_f32.c
SRC += $(DSP_LIB)/Source/TransformFunctions/arm_cfft_radix8_f32.c
SRC += $(DSP_LIB)/Source/CommonTables/arm_common_tables.c

SRC += $(DSP_LIB)/Source/ComplexMathFunctions/arm_cmplx_mag_f32.c
SRC += $(DSP_LIB)/Source/StatisticsFunctions/arm_max_f32.c

SRC += $(wildcard $(DSP_LIB)/Source/*/*.S)
Expand Down
5 changes: 4 additions & 1 deletion src/main/blackbox/blackbox.c
Original file line number Diff line number Diff line change
Expand Up @@ -1176,6 +1176,8 @@ static bool blackboxWriteSysinfo(void) {
BLACKBOX_PRINT_HEADER_LINE("dterm_lowpass2_hz_roll", "%d", currentPidProfile->dFilter[ROLL].dLpf2);
BLACKBOX_PRINT_HEADER_LINE("dterm_lowpass2_hz_pitch", "%d", currentPidProfile->dFilter[PITCH].dLpf2);
BLACKBOX_PRINT_HEADER_LINE("dterm_lowpass2_hz_yaw", "%d", currentPidProfile->dFilter[YAW].dLpf2);
BLACKBOX_PRINT_HEADER_LINE("dterm_dyn_notch_enable", "%d", currentPidProfile->dtermDynNotch);
BLACKBOX_PRINT_HEADER_LINE("dterm_dyn_notch_q", "%d", currentPidProfile->dterm_dyn_notch_q);
BLACKBOX_PRINT_HEADER_LINE("dterm_ABG_alpha", "%d", currentPidProfile->dterm_ABG_alpha);
BLACKBOX_PRINT_HEADER_LINE("dterm_ABG_boost", "%d", currentPidProfile->dterm_ABG_boost);
BLACKBOX_PRINT_HEADER_LINE("dterm_ABG_half_life", "%d", currentPidProfile->dterm_ABG_half_life);
Expand Down Expand Up @@ -1230,7 +1232,8 @@ static bool blackboxWriteSysinfo(void) {
BLACKBOX_PRINT_HEADER_LINE("gyro_notch_cutoff", "%d,%d", gyroConfig()->gyro_soft_notch_cutoff_1,
gyroConfig()->gyro_soft_notch_cutoff_2);
#if defined(USE_GYRO_DATA_ANALYSE)
BLACKBOX_PRINT_HEADER_LINE("dynamic_gyro_notch_q", "%d", gyroConfig()->dyn_notch_q_factor);
BLACKBOX_PRINT_HEADER_LINE("dynamic_gyro_width", "%d", gyroConfig()->dyn_notch_q);
BLACKBOX_PRINT_HEADER_LINE("dynamic_gyro_notch_count", "%d", gyroConfig()->dyn_notch_count);
BLACKBOX_PRINT_HEADER_LINE("dynamic_gyro_notch_min_hz", "%d", gyroConfig()->dyn_notch_min_hz);
BLACKBOX_PRINT_HEADER_LINE("dynamic_gyro_notch_max_hz", "%d", gyroConfig()->dyn_notch_max_hz);
#endif
Expand Down
52 changes: 24 additions & 28 deletions src/main/cms/cms_menu_imu.c
Original file line number Diff line number Diff line change
Expand Up @@ -490,13 +490,10 @@ static uint16_t gyroConfig_gyro_lowpass_hz_yaw;
static uint16_t gyroConfig_gyro_lowpass2_hz_roll;
static uint16_t gyroConfig_gyro_lowpass2_hz_pitch;
static uint16_t gyroConfig_gyro_lowpass2_hz_yaw;
static uint16_t gyroConfig_gyro_soft_notch_hz_1;
static uint16_t gyroConfig_gyro_soft_notch_cutoff_1;
static uint16_t gyroConfig_gyro_soft_notch_hz_2;
static uint16_t gyroConfig_gyro_soft_notch_cutoff_2;
static uint16_t gyroConfig_gyro_matrix_q;
static uint16_t gyroConfig_gyro_matrix_min_hz;
static uint16_t gyroConfig_gyro_matrix_max_hz;
static uint16_t gyroConfig_gyro_q;
static uint8_t gyroConfig_gyro_notch_count;
static uint16_t gyroConfig_gyro_notch_min_hz;
static uint16_t gyroConfig_gyro_notch_max_hz;
static uint16_t gyroConfig_gyro_abg_alpha;
static uint16_t gyroConfig_gyro_abg_boost;
static uint8_t gyroConfig_gyro_abg_half_life;
Expand All @@ -521,13 +518,10 @@ static long cmsx_menuGyro_onEnter(void) {
gyroConfig_gyro_lowpass2_hz_roll = gyroConfig()->gyro_lowpass2_hz[ROLL];
gyroConfig_gyro_lowpass2_hz_pitch = gyroConfig()->gyro_lowpass2_hz[PITCH];
gyroConfig_gyro_lowpass2_hz_yaw = gyroConfig()->gyro_lowpass2_hz[YAW];
gyroConfig_gyro_soft_notch_hz_1 = gyroConfig()->gyro_soft_notch_hz_1;
gyroConfig_gyro_soft_notch_cutoff_1 = gyroConfig()->gyro_soft_notch_cutoff_1;
gyroConfig_gyro_soft_notch_hz_2 = gyroConfig()->gyro_soft_notch_hz_2;
gyroConfig_gyro_soft_notch_cutoff_2 = gyroConfig()->gyro_soft_notch_cutoff_2;
gyroConfig_gyro_matrix_q = gyroConfig()->dyn_notch_q_factor;
gyroConfig_gyro_matrix_min_hz = gyroConfig()->dyn_notch_min_hz;
gyroConfig_gyro_matrix_max_hz = gyroConfig()->dyn_notch_max_hz;
gyroConfig_gyro_q = gyroConfig()->dyn_notch_q;
gyroConfig_gyro_notch_count = gyroConfig()->dyn_notch_count;
gyroConfig_gyro_notch_min_hz = gyroConfig()->dyn_notch_min_hz;
gyroConfig_gyro_notch_max_hz = gyroConfig()->dyn_notch_max_hz;
gyroConfig_gyro_abg_alpha = gyroConfig()->gyro_ABG_alpha;
gyroConfig_gyro_abg_boost = gyroConfig()->gyro_ABG_boost;
gyroConfig_gyro_abg_half_life = gyroConfig()->gyro_ABG_half_life;
Expand Down Expand Up @@ -556,13 +550,10 @@ static long cmsx_menuGyro_onExit(const OSD_Entry *self) {
gyroConfigMutable()->gyro_lowpass2_hz[ROLL] = gyroConfig_gyro_lowpass2_hz_roll;
gyroConfigMutable()->gyro_lowpass2_hz[PITCH] = gyroConfig_gyro_lowpass2_hz_pitch;
gyroConfigMutable()->gyro_lowpass2_hz[YAW] = gyroConfig_gyro_lowpass2_hz_yaw;
gyroConfigMutable()->gyro_soft_notch_hz_1 = gyroConfig_gyro_soft_notch_hz_1;
gyroConfigMutable()->gyro_soft_notch_cutoff_1 = gyroConfig_gyro_soft_notch_cutoff_1;
gyroConfigMutable()->gyro_soft_notch_hz_2 = gyroConfig_gyro_soft_notch_hz_2;
gyroConfigMutable()->gyro_soft_notch_cutoff_2 = gyroConfig_gyro_soft_notch_cutoff_2;
gyroConfigMutable()->dyn_notch_q_factor = gyroConfig_gyro_matrix_q;
gyroConfigMutable()->dyn_notch_min_hz = gyroConfig_gyro_matrix_min_hz;
gyroConfigMutable()->dyn_notch_max_hz = gyroConfig_gyro_matrix_max_hz;
gyroConfigMutable()->dyn_notch_q = gyroConfig_gyro_q;
gyroConfigMutable()->dyn_notch_count = gyroConfig_gyro_notch_count;
gyroConfigMutable()->dyn_notch_min_hz = gyroConfig_gyro_notch_min_hz;
gyroConfigMutable()->dyn_notch_max_hz = gyroConfig_gyro_notch_max_hz;
gyroConfigMutable()->gyro_ABG_alpha = gyroConfig_gyro_abg_alpha;
gyroConfigMutable()->gyro_ABG_boost = gyroConfig_gyro_abg_boost;
gyroConfigMutable()->gyro_ABG_half_life = gyroConfig_gyro_abg_half_life;
Expand Down Expand Up @@ -593,20 +584,17 @@ static OSD_Entry cmsx_menuFilterGlobalEntries[] = {
{ "GYRO LPF2 PITCH", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_lowpass2_hz_pitch, 0, 16000, 1 }, 0 },
{ "GYRO LPF2 YAW", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_lowpass2_hz_yaw, 0, 16000, 1 }, 0 },
#endif
{ "MATRIX Q", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_matrix_q, 0, 1000, 1 }, 0 },
{ "MATRIX MIN HZ", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_matrix_min_hz, 30, 1000, 1 }, 0 },
{ "MATRIX MAX HZ", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_matrix_max_hz, 200, 1000, 1 }, 0 },
{ "DYN NOTCH Q", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_q, 0, 1000, 1 }, 0 },
{ "DYN NOTCH COUNT", OME_UINT8, NULL, &(OSD_UINT8_t) { &gyroConfig_gyro_notch_count, 1, 5, 1 }, 0 },
{ "DYN NOTCH MIN HZ", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_notch_min_hz, 30, 1000, 1 }, 0 },
{ "DYN NOTCH MAX HZ", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_notch_max_hz, 200, 1000, 1 }, 0 },
#ifndef USE_GYRO_IMUF9001
{ "IMUF W", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_imuf_w, 0, 512, 1 }, 0 },
{ "ROLL Q", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_imuf_roll_q, 100, 16000, 100 }, 0 },
{ "PITCH Q", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_imuf_pitch_q, 100, 16000, 100 }, 0 },
{ "YAW Q", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_imuf_yaw_q, 100, 16000, 100 }, 0 },

#endif
{ "GYRO NF1", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_soft_notch_hz_1, 0, 500, 1 }, 0 },
{ "GYRO NF1C", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_soft_notch_cutoff_1, 0, 500, 1 }, 0 },
{ "GYRO NF2", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_soft_notch_hz_2, 0, 500, 1 }, 0 },
{ "GYRO NF2C", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_soft_notch_cutoff_2, 0, 500, 1 }, 0 },

{ "GYRO ABG ALPHA", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_abg_alpha, 0, 1000, 1 }, 0 },
{ "GYRO ABG BOOST", OME_UINT16, NULL, &(OSD_UINT16_t) { &gyroConfig_gyro_abg_boost, 0, 2000, 1 }, 0 },
Expand Down Expand Up @@ -721,6 +709,8 @@ static uint16_t cmsx_dterm_lowpass2_type;
static uint16_t cmsx_dterm_lowpass2_hz_roll;
static uint16_t cmsx_dterm_lowpass2_hz_pitch;
static uint16_t cmsx_dterm_lowpass2_hz_yaw;
static uint8_t cmsx_dterm_dyn_notch_enable ;
static uint16_t cmsx_dterm_dyn_notch_q;
static uint16_t cmsx_dterm_abg_alpha;
static uint16_t cmsx_dterm_abg_boost;
static uint8_t cmsx_dterm_abg_half_life;
Expand All @@ -735,6 +725,8 @@ static long cmsx_FilterPerProfileRead(void) {
cmsx_dterm_lowpass2_hz_roll = pidProfile->dFilter[ROLL].dLpf2;
cmsx_dterm_lowpass2_hz_pitch = pidProfile->dFilter[PITCH].dLpf2;
cmsx_dterm_lowpass2_hz_yaw = pidProfile->dFilter[YAW].dLpf2;
cmsx_dterm_dyn_notch_enable = pidProfile->dtermDynNotch;
cmsx_dterm_dyn_notch_q = pidProfile->dterm_dyn_notch_q;
cmsx_dterm_abg_alpha = pidProfile->dterm_ABG_alpha;
cmsx_dterm_abg_boost = pidProfile->dterm_ABG_boost;
cmsx_dterm_abg_half_life = pidProfile->dterm_ABG_half_life;
Expand All @@ -752,6 +744,8 @@ static long cmsx_FilterPerProfileWriteback(const OSD_Entry *self) {
pidProfile->dFilter[ROLL].dLpf2 = cmsx_dterm_lowpass2_hz_roll;
pidProfile->dFilter[PITCH].dLpf2 = cmsx_dterm_lowpass2_hz_pitch;
pidProfile->dFilter[YAW].dLpf2 = cmsx_dterm_lowpass2_hz_yaw;
pidProfile->dtermDynNotch = cmsx_dterm_dyn_notch_enable;
pidProfile->dterm_dyn_notch_q = cmsx_dterm_dyn_notch_q;
pidProfile->dterm_ABG_alpha = cmsx_dterm_abg_alpha;
pidProfile->dterm_ABG_boost = cmsx_dterm_abg_boost;
pidProfile->dterm_ABG_half_life = cmsx_dterm_abg_half_life;
Expand All @@ -769,6 +763,8 @@ static OSD_Entry cmsx_menuFilterPerProfileEntries[] = {
{ "DTERM LPF2 ROLL", OME_UINT16, NULL, &(OSD_UINT16_t){ &cmsx_dterm_lowpass2_hz_roll, 0, 500, 1 }, 0 },
{ "DTERM LPF2 PITCH", OME_UINT16, NULL, &(OSD_UINT16_t){ &cmsx_dterm_lowpass2_hz_pitch, 0, 500, 1 }, 0 },
{ "DTERM LPF2 YAW", OME_UINT16, NULL, &(OSD_UINT16_t){ &cmsx_dterm_lowpass2_hz_yaw, 0, 500, 1 }, 0 },
{ "DTERM DYN ENABLE", OME_TAB, NULL, &(OSD_TAB_t) { (uint8_t *) &cmsx_dterm_dyn_notch_enable, 1, cms_offOnLabels }, 0 },
{ "DTERM DYN NOT Q", OME_UINT16, NULL, &(OSD_UINT16_t){ &cmsx_dterm_dyn_notch_q, 0, 2000, 1 }, 0 },
{ "DTERM ABG ALPHA", OME_UINT16, NULL, &(OSD_UINT16_t){ &cmsx_dterm_abg_alpha, 0, 1000, 1 }, 0 },
{ "DTERM ABG BOOST", OME_UINT16, NULL, &(OSD_UINT16_t){ &cmsx_dterm_abg_boost, 0, 2000, 1 }, 0 },
{ "DTERM ABG HL", OME_UINT8, NULL, &(OSD_UINT8_t) { &cmsx_dterm_abg_half_life, 0, 250, 1 }, 0 },
Expand Down
10 changes: 4 additions & 6 deletions src/main/common/filter.c
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,6 @@

#include "fc/fc_rc.h"

#define M_LN2_FLOAT 0.69314718055994530942f
#define M_PI_FLOAT 3.14159265358979323846f
#define BIQUAD_Q (1.0f / sqrtf(2.0f)) /* quality factor - 2nd order butterworth*/

// NULL filter
Expand All @@ -46,7 +44,7 @@ FAST_CODE float nullFilterApply(filter_t *filter, float input) {
// PT1 Low Pass filter

float pt1FilterGain(uint16_t f_cut, float dT) {
const float RC = 0.5f / (M_PI_FLOAT * f_cut);
const float RC = 0.5f / (M_PIf * f_cut);
return dT / (RC + dT);
}

Expand Down Expand Up @@ -108,7 +106,7 @@ void biquadFilterInit(biquadFilter_t *filter, float filterFreq, uint32_t refresh

FAST_CODE void biquadFilterUpdate(biquadFilter_t *filter, float filterFreq, uint32_t refreshRate, float Q, biquadFilterType_e filterType) {
// setup variables
const float omega = 2.0f * M_PI_FLOAT * filterFreq * refreshRate * 0.000001f;
const float omega = 2.0f * M_PIf * filterFreq * refreshRate * 0.000001f;
const float sn = sin_approx(omega);
const float cs = cos_approx(omega);
const float alpha = sn / (2.0f * Q);
Expand Down Expand Up @@ -257,13 +255,13 @@ FAST_CODE void ptnFilterInit(ptnFilter_t *filter, uint8_t order, uint16_t f_cut,

Adj_f_cut = (float)f_cut * ScaleF[filter->order - 1];

filter->k = dT / ((1.0f / (2.0f * M_PI_FLOAT * Adj_f_cut)) + dT);
filter->k = dT / ((1.0f / (2.0f * M_PIf * Adj_f_cut)) + dT);
} // ptnFilterInit

FAST_CODE void ptnFilterUpdate(ptnFilter_t *filter, float f_cut, float ScaleF, float dT) {
float Adj_f_cut;
Adj_f_cut = (float)f_cut * ScaleF;
filter->k = dT / ((1.0f / (2.0f * M_PI_FLOAT * Adj_f_cut)) + dT);
filter->k = dT / ((1.0f / (2.0f * M_PIf * Adj_f_cut)) + dT);
}

FAST_CODE float ptnFilterApply(ptnFilter_t *filter, float input) {
Expand Down
3 changes: 3 additions & 0 deletions src/main/common/maths.h
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@
#define M_PIf 3.14159265358979323846f
#define M_PI_HALFf 3.14159265358979323846f/2
#define M_EULERf 2.71828182845904523536f
#define M_SQRT2f 1.41421356237309504880f
#define M_LN2f 0.69314718055994530942f


#define RAD (M_PIf / 180.0f)
#define DEGREES_TO_DECIDEGREES(angle) ((angle) * 10)
Expand Down
153 changes: 153 additions & 0 deletions src/main/common/sdft.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
/*
* This file is part of Cleanflight and Betaflight.
*
* Cleanflight and Betaflight are free software. You can redistribute
* this software and/or modify this software 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.
*
* Cleanflight and Betaflight are distributed in the hope that they
* 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 software.
*
* If not, see <http://www.gnu.org/licenses/>.
*/

#include <stdbool.h>
#include <math.h>

#include "platform.h"

#include "common/maths.h"
#include "common/sdft.h"

#define SDFT_R 0.9999f // damping factor for guaranteed SDFT stability (r < 1.0f)

static FAST_RAM_ZERO_INIT float rPowerN; // SDFT_R to the power of SDFT_SAMPLE_SIZE
static FAST_RAM_ZERO_INIT bool isInitialized;
static FAST_RAM_ZERO_INIT complex_t twiddle[SDFT_BIN_COUNT];

static void applySqrt(const sdft_t *sdft, float *data);


void sdftInit(sdft_t *sdft, const uint8_t startBin, const uint8_t endBin, const uint8_t numBatches)
{
if (!isInitialized) {
rPowerN = powf(SDFT_R, SDFT_SAMPLE_SIZE);
const float c = 2.0f * M_PIf / (float)SDFT_SAMPLE_SIZE;
for (uint8_t i = 0; i < SDFT_BIN_COUNT; i++) {
float phi = 0.0f;
phi = c * i;
twiddle[i] = SDFT_R * (cos_approx(phi) + _Complex_I * sin_approx(phi));
}
isInitialized = true;
}

sdft->idx = 0;
sdft->startBin = startBin;
sdft->endBin = endBin;
sdft->numBatches = numBatches;
sdft->batchSize = (sdft->endBin - sdft->startBin + 1) / sdft->numBatches + 1;

for (uint8_t i = 0; i < SDFT_SAMPLE_SIZE; i++) {
sdft->samples[i] = 0.0f;
}

for (uint8_t i = 0; i < SDFT_BIN_COUNT; i++) {
sdft->data[i] = 0.0f;
}
}


// Add new sample to frequency spectrum
FAST_CODE void sdftPush(sdft_t *sdft, const float *sample)
{
const float delta = *sample - rPowerN * sdft->samples[sdft->idx];

sdft->samples[sdft->idx] = *sample;
sdft->idx = (sdft->idx + 1) % SDFT_SAMPLE_SIZE;

for (uint8_t i = sdft->startBin; i <= sdft->endBin; i++) {
sdft->data[i] = twiddle[i] * (sdft->data[i] + delta);
}
}


// Add new sample to frequency spectrum in parts
FAST_CODE void sdftPushBatch(sdft_t* sdft, const float *sample, const uint8_t *batchIdx)
{
const uint8_t batchStart = sdft->batchSize * *batchIdx;
uint8_t batchEnd = batchStart;

const float delta = *sample - rPowerN * sdft->samples[sdft->idx];

if (*batchIdx == sdft->numBatches - 1) {
sdft->samples[sdft->idx] = *sample;
sdft->idx = (sdft->idx + 1) % SDFT_SAMPLE_SIZE;
batchEnd += sdft->endBin - batchStart + 1;
} else {
batchEnd += sdft->batchSize;
}

for (uint8_t i = batchStart; i < batchEnd; i++) {
sdft->data[i] = twiddle[i] * (sdft->data[i] + delta);
}
}


// Get squared magnitude of frequency spectrum
FAST_CODE void sdftMagSq(const sdft_t *sdft, float *output)
{
for (uint8_t i = sdft->startBin; i <= sdft->endBin; i++) {
float re = crealf(sdft->data[i]);
float im = cimagf(sdft->data[i]);
output[i] = re * re + im * im;
}
}


// Get magnitude of frequency spectrum (slower)
FAST_CODE void sdftMagnitude(const sdft_t *sdft, float *output)
{
sdftMagSq(sdft, output);
applySqrt(sdft, output);
}


// Get squared magnitude of frequency spectrum with Hann window applied
// Hann window in frequency domain: X[k] = -0.25 * X[k-1] +0.5 * X[k] -0.25 * X[k+1]
FAST_CODE void sdftWinSq(const sdft_t *sdft, float *output)
{
complex_t val;

for (uint8_t i = (sdft->startBin + 1); i < sdft->endBin; i++) {
float re;
float im;
val = sdft->data[i] - 0.5f * (sdft->data[i - 1] + sdft->data[i + 1]); // multiply by 2 to save one multiplication
re = crealf(val);
im = cimagf(val);
output[i] = re * re + im * im;
}
}


// Get magnitude of frequency spectrum with Hann window applied (slower)
FAST_CODE void sdftWindow(const sdft_t *sdft, float *output)
{
sdftWinSq(sdft, output);
applySqrt(sdft, output);
}


// Apply fast square root approximation to the whole sdft range
static FAST_CODE void applySqrt(const sdft_t *sdft, float *data)
{
for (uint8_t i = sdft->startBin; i <= sdft->endBin; i++) {
data[i] = sqrtf(data[i]);
}
}
Loading

0 comments on commit 75a10b7

Please sign in to comment.