From 8c620f88e4626282f596b3544f5741c7f4c2ef06 Mon Sep 17 00:00:00 2001 From: daniel-c-harvey Date: Wed, 29 Jul 2026 23:36:00 -0400 Subject: [PATCH] temp-cortex filter dump --- src/temp-cortex/biquad.hpp | 73 ++++++++++++++ src/temp-cortex/biquad.tpp | 195 +++++++++++++++++++++++++++++++++++++ src/temp-cortex/filter.hpp | 62 ++++++++++++ src/temp-cortex/util.hpp | 46 +++++++++ 4 files changed, 376 insertions(+) create mode 100644 src/temp-cortex/biquad.hpp create mode 100644 src/temp-cortex/biquad.tpp create mode 100644 src/temp-cortex/filter.hpp create mode 100644 src/temp-cortex/util.hpp diff --git a/src/temp-cortex/biquad.hpp b/src/temp-cortex/biquad.hpp new file mode 100644 index 0000000..b8c916e --- /dev/null +++ b/src/temp-cortex/biquad.hpp @@ -0,0 +1,73 @@ +#pragma once + +#include "util.hpp" +#include "filter.hpp" +#include "filter_params.hpp" + +template +class Biqaud : public Filter +{ + public: + Biqaud(const uint32_t& sample_rate, FilterParameters *params); + + void prepare_parameters(const TUIParams& params) override; + + protected: + uint32_t sample_rate; + + void process_channel_frame(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) override; + + void filter(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) override; + + void update_feedback(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) override; +}; + +template +class BiquadHP : public Biqaud +{ + public: + BiquadHP(const uint32_t& sample_rate, FilterParameters *params); + + NormalCoefficients prepare_coefficients() override; + + protected: + void process_channel_frame(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) override; +}; + +template +class BiquadLP : public Biqaud +{ + public: + BiquadLP(const uint32_t& sample_rate, FilterParameters *params); + + NormalCoefficients prepare_coefficients() override; +}; + +template +class Biquad1PoleLP : public BiquadLP +{ + public: + Biquad1PoleLP(const uint32_t& sample_rate, FilterParameters *params); + + NormalCoefficients prepare_coefficients() override; + + protected: + void process_channel_frame(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) override; +}; + +#include "biquad.tpp" \ No newline at end of file diff --git a/src/temp-cortex/biquad.tpp b/src/temp-cortex/biquad.tpp new file mode 100644 index 0000000..e193be0 --- /dev/null +++ b/src/temp-cortex/biquad.tpp @@ -0,0 +1,195 @@ +#pragma once + +#include "basicmaths.h" +#include "biquad.hpp" + +template +Biqaud::Biqaud(const uint32_t& p_sample_rate, FilterParameters *p_params) +: Filter(p_params), sample_rate(p_sample_rate) +{ + // Initialize filter state to zero to prevent random behavior + for (int i = 0; i < k_channels; i++) { + this->state[i].x[0] = 0.0f; + this->state[i].x[1] = 0.0f; + this->state[i].y[0] = 0.0f; + this->state[i].y[1] = 0.0f; + this->state[i].fb = 0.0f; + } +} + +template +void Biqaud::prepare_parameters(const TUIParams& params) +{ + // Direct logarithmic interpolation for smooth frequency scaling using standard math + const float min_freq = 10.f; + const float max_freq = 23000.f; + float log_freq = logf(min_freq) + params.p_cutoff * (logf(max_freq) - logf(min_freq)); + float raw_cutoff = expf(log_freq); + this->params->cutoff = fminf(raw_cutoff, 0.48f * this->sample_rate); // Allow closer to Nyquist + + // Resonance response + this->params->res = params.p_resonance; + + // Base Q of 0.707 plus resonance + this->params->Q = M_SQRT1_2 + this->params->res; +} + +template +void Biqaud::process_channel_frame(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) +{ + this->filter(state, coeff, x, y); + this->update_feedback(state, coeff, x, y); +} + +template +void Biqaud::filter(FeedbackLine &state, const NormalCoefficients &coeff, const float &x, float &y) +{ + // debugMessage("Biqaud::filter"); + // Direct Form I biquad - matches Audio EQ Cookbook exactly + y = coeff.b0 * x + coeff.b1 * state.x[0] + coeff.b2 * state.x[1] + - coeff.a1 * state.y[0] - coeff.a2 * state.y[1]; +} + +template +void Biqaud::update_feedback(FeedbackLine& state, const NormalCoefficients& coeff, const float& x, float& y) +{ + // debugMessage("State x[0], x[1], y[0]: ", state.x[0], state.x[1], state.y[0]); + // Update feedback state + state.x[1] = state.x[0]; + state.x[0] = x; + state.y[1] = state.y[0]; + state.y[0] = y; + state.fb = y; +} + +template +BiquadHP::BiquadHP(const uint32_t& p_sample_rate, FilterParameters *p_params) +: Biqaud(p_sample_rate, p_params) {} + +template +void BiquadHP::process_channel_frame(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) +{ + // CRITICAL: Highpass filters require input feedback to work properly + // This compensates for coefficient collapse at low frequencies + const float fb_amount = this->params->res * 0.24f; + float input = x - fb_amount * feedback_saturate(state.fb * 0.9f); + + Biqaud::process_channel_frame(state, coeff, input, y); +} + +template +NormalCoefficients BiquadHP::prepare_coefficients() +{ + // Pre-warped bilinear transform - same topology as lowpass but for highpass + const float w = tanf(M_PI * this->params->cutoff / this->sample_rate); + const float w2 = w * w; + const float cosw = (1.0f - w2) / (1.0f + w2); + const float sinw = 2.0f * w / (1.0f + w2); + const float alpha = sinw / (2.0f * this->params->Q); + + // Standard RBJ highpass with pre-warped frequency + const float norm = 1.0f / (1.0f + alpha); + const float b0 = (1.0f + cosw) * 0.5f * norm; + const float b1 = -(1.0f + cosw) * norm; + const float b2 = (1.0f + cosw) * 0.5f * norm; + const float a1 = -2.0f * cosw * norm; + const float a2 = (1.0f - alpha) * norm; + + NormalCoefficients coeff = { + .a1 = a1, + .a2 = a2, + .b0 = b0, + .b1 = b1, + .b2 = b2 + }; + + return coeff; +} + +template +BiquadLP::BiquadLP(const uint32_t& p_sample_rate, FilterParameters *p_params) +: Biqaud(p_sample_rate, p_params) {} + +template +NormalCoefficients BiquadLP::prepare_coefficients() +{ + // Pre-warped bilinear transform - correct implementation + const float w = tanf(M_PI * this->params->cutoff / this->sample_rate); + const float w2 = w * w; + const float cosw = (1.0f - w2) / (1.0f + w2); + const float sinw = 2.0f * w / (1.0f + w2); + const float alpha = sinw / (2.0f * this->params->Q); + + // Standard RBJ lowpass with pre-warped frequency + const float norm = 1.0f / (1.0f + alpha); + const float b0 = (1.0f - cosw) * 0.5f * norm; + const float b1 = (1.0f - cosw) * norm; + const float b2 = (1.0f - cosw) * 0.5f * norm; + const float a1 = -2.0f * cosw * norm; + const float a2 = (1.0f - alpha) * norm; + + NormalCoefficients coeff = { + .a1 = a1, + .a2 = a2, + .b0 = b0, + .b1 = b1, + .b2 = b2 + }; + + return coeff; +} + +template +Biquad1PoleLP::Biquad1PoleLP(const uint32_t& p_sample_rate, FilterParameters *p_params) +: BiquadLP(p_sample_rate, p_params) {} + +template +void Biquad1PoleLP::process_channel_frame(FeedbackLine& state, + const NormalCoefficients& coeff, + const float& x, + float& y) +{ + // Stable Moog-style feedback with conservative limits + // Much more conservative k values for single-pole stability + const float k_max = 3.8f; // Much lower max for stability + const float k = fminf(k_max, fmaxf(0.0f, (this->params->Q - M_SQRT1_2))); // Conservative Q mapping + + // Conservative gain compensation + const float makeup_gain = 1.0f + k * 0.5f; // Gentler compensation + + // Stable global feedback with limiting + float resonant_input = (x - k * feedback_saturate(state.fb * 0.8f)) * makeup_gain; + + // Process with stable resonant input + Biqaud::process_channel_frame(state, coeff, resonant_input, y); +} + +template +NormalCoefficients Biquad1PoleLP::prepare_coefficients() +{ + // Correct 1-pole lowpass using bilinear transform + // H(s) = wc/(s + wc) -> H(z) = b0*(1+z^-1)/(1 + a1*z^-1) + const float w = tanf(M_PI * this->params->cutoff / this->sample_rate); + + // Bilinear transform gives both b0 and b1 coefficients + const float norm = 1.0f / (1.0f + w); + const float b0 = w * norm; // Coefficient for x[n] + const float b1 = w * norm; // Coefficient for x[n-1] (same as b0) + const float a1 = (w - 1.0f) * norm; // Pole coefficient + + NormalCoefficients coeff = { + .a1 = a1, // Pole coefficient + .a2 = 0.0f, // 1-pole has no second pole + .b0 = b0, // Current input coefficient + .b1 = b1, // Previous input coefficient + .b2 = 0.0f // 1-pole has no z^-2 numerator + }; + + return coeff; +} \ No newline at end of file diff --git a/src/temp-cortex/filter.hpp b/src/temp-cortex/filter.hpp new file mode 100644 index 0000000..ccc3d0a --- /dev/null +++ b/src/temp-cortex/filter.hpp @@ -0,0 +1,62 @@ +#pragma once + +#include "util.hpp" +// #include "fdecorator.hpp" + +template +class FilterBase +{ + public: + FilterBase(TFilterParams *p) : params(p) {} + + /// @brief Prepare the filter channels to process all frames in this block + virtual void prepare_parameters(const TUIParams& params) = 0; + + /// @brief Prepare the filter channels to process all frames in this block + virtual TCoefficients prepare_coefficients() = 0; + + /// @brief process the current frame samples for all channels + /// @param x inputs samples + /// @param y output samples + virtual void process_frame(const TCoefficients& coeff, const float x[k_channels], float y[k_channels]) + { + // Handle channel iteration in the base class to ensure virtual dispatch through decorator chain + for (uint16_t channel = 0; channel < k_channels; channel++) { + this->process_channel_frame(this->state[channel], coeff, x[channel], y[channel]); + } + } + + protected: + /// @brief process the current frame sample for given channel + /// @param x inputs sample + /// @param y output sample + virtual void process_channel_frame(TFeedbackLine& state, const TCoefficients& coeff, const float& x, float& y) = 0; + + /// @brief filter the current frame sample for given channel + /// @param state filter state + /// @param coeff filter coefficients + /// @param x input sample + /// @param y output sample + virtual void filter(TFeedbackLine& state, const TCoefficients& coeff, const float& x, float& y) = 0; + + /// @brief update the feedback line for the next frame + /// @param state filter state + /// @param coeff filter coefficients + /// @param x input sample + /// @param y output sample + virtual void update_feedback(TFeedbackLine& state, const TCoefficients& coeff, const float& x, float& y) = 0; + + TFilterParams* params; + TFeedbackLine state[k_channels]; + + template + friend class FilterDecorator; +}; + +template +class Filter : public FilterBase +{ + public: + Filter(TFilterParams *p) + : FilterBase(p) {} +}; \ No newline at end of file diff --git a/src/temp-cortex/util.hpp b/src/temp-cortex/util.hpp new file mode 100644 index 0000000..88e0b42 --- /dev/null +++ b/src/temp-cortex/util.hpp @@ -0,0 +1,46 @@ +#pragma once +#include "basicmaths.h" + +// Improved tanh approximation with proper continuity +static inline float tanh_saturate(float x, float threshold, float a, float b) +{ + if (x > threshold) { + float excess = x - threshold; + float sat_val = threshold * a / (a + b + threshold * threshold); // Value at threshold + return sat_val + excess * 0.1f; // Gentle slope beyond threshold + } + if (x < -threshold) { + float excess = x + threshold; + float sat_val = -threshold * a / (a + b + threshold * threshold); // Value at -threshold + return sat_val + excess * 0.1f; // Gentle slope beyond -threshold + } + const float x2 = x * x; + return x * a / (a + b + x2); +} + +// TB-303 style feedback saturation +// Hard saturation for filter feedback (handles large values) +static inline float feedback_saturate(float x) +{ + // More aggressive saturation for feedback control + return tanh_saturate(x, 2.0f, 27.f, 9.f); +} + +// Gentle saturation for audio signals (subtle, musical) +static inline float audio_saturate(float x) +{ + // Adjusted parameters to maintain more volume at threshold + // At x=0.92: output ≈ 0.85 (much better than previous 0.57) + return tanh_saturate(x, 0.92f, 15.0f, 1.0f); +} + +/// @brief Tunable logistic function (sigmoid) +/// @param a slope +/// @param b slope 2 +/// @param c offset +/// @param z portion scalar +/// @return H(x) +static inline float H(float x, float a, float b, float c, float z) +{ + return z * a / (a + expf(b * (c - x))) - 0.02f; +} \ No newline at end of file