temp-cortex filter dump

2026-07-29 23:36:00 -04:00
parent a689fb75eb
commit 8c620f88e4
4 changed files with 376 additions and 0 deletions
@@ -0,0 +1,73 @@
+#pragma once
+
+#include "util.hpp"
+#include "filter.hpp"
+#include "filter_params.hpp"
+
+template <int k_channels, typename TUIParams>
+class Biqaud : public Filter<k_channels, FeedbackLine, NormalCoefficients, TUIParams, FilterParameters>
+{
+    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 <int k_channels, typename TUIParams>
+class BiquadHP : public Biqaud<k_channels, TUIParams>
+{
+    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 <int k_channels, typename TUIParams>
+class BiquadLP : public Biqaud<k_channels, TUIParams>
+{
+    public:
+        BiquadLP(const uint32_t& sample_rate, FilterParameters *params);
+
+        NormalCoefficients prepare_coefficients() override;
+};
+
+template <int k_channels, typename TUIParams>
+class Biquad1PoleLP : public BiquadLP<k_channels, TUIParams>
+{
+    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"
@@ -0,0 +1,195 @@
+#pragma once
+
+#include "basicmaths.h"
+#include "biquad.hpp"
+
+template <int k_channels, typename TUIParams>
+Biqaud<k_channels, TUIParams>::Biqaud(const uint32_t& p_sample_rate, FilterParameters *p_params)
+: Filter<k_channels, FeedbackLine, NormalCoefficients, TUIParams, FilterParameters>(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 <int k_channels, typename TUIParams>
+void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
+void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
+void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
+void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
+BiquadHP<k_channels, TUIParams>::BiquadHP(const uint32_t& p_sample_rate, FilterParameters *p_params)
+: Biqaud<k_channels, TUIParams>(p_sample_rate, p_params) {}
+
+template <int k_channels, typename TUIParams>
+void BiquadHP<k_channels, TUIParams>::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<k_channels, TUIParams>::process_channel_frame(state, coeff, input, y);
+}
+
+template <int k_channels, typename TUIParams>
+NormalCoefficients BiquadHP<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
+BiquadLP<k_channels, TUIParams>::BiquadLP(const uint32_t& p_sample_rate, FilterParameters *p_params)
+: Biqaud<k_channels, TUIParams>(p_sample_rate, p_params) {}
+
+template <int k_channels, typename TUIParams>
+NormalCoefficients BiquadLP<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
+Biquad1PoleLP<k_channels, TUIParams>::Biquad1PoleLP(const uint32_t& p_sample_rate, FilterParameters *p_params)
+: BiquadLP<k_channels, TUIParams>(p_sample_rate, p_params) {}
+
+template <int k_channels, typename TUIParams>
+void Biquad1PoleLP<k_channels, TUIParams>::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<k_channels, TUIParams>::process_channel_frame(state, coeff, resonant_input, y);
+}
+
+template <int k_channels, typename TUIParams>
+NormalCoefficients Biquad1PoleLP<k_channels, TUIParams>::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;
+}
@@ -0,0 +1,62 @@
+#pragma once
+
+#include "util.hpp"
+// #include "fdecorator.hpp"
+
+template <int k_channels, typename TFeedbackLine, typename TCoefficients, typename TUIParams, typename TFilterParams>
+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 <int, typename, typename, typename, typename, typename, typename> 
+        friend class FilterDecorator;
+};
+
+template <int k_channels, typename TFeedbackLine, typename TCoefficients, typename TUIParams, typename TFilterParams>
+class Filter : public FilterBase<k_channels, TFeedbackLine, TCoefficients, TUIParams, TFilterParams>
+{
+    public:  
+        Filter(TFilterParams *p) 
+        : FilterBase<k_channels, TFeedbackLine, TCoefficients, TUIParams, TFilterParams>(p) {}
+};
@@ -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;
+}