Use stable Zero-Delay Feedback (ZDF) SVF topology for filters

Author: juzzlin

src/domain/devices/high_pass_filter_effect.cpp

@@ -32,25 +32,48 @@ void HighPassFilterEffect::process(float & left, float & right, uint32_t sampleR
         return;
     }
 
+    // Zero-Delay Feedback State Variable Filter (2nd order)
     const float freq = 20.0f * std::pow(std::min(20000.0f, sampleRate * 0.49f) / 20.0f, m_cutoff);
-    const float f = std::min(1.0f, 2.0f * std::sin(std::numbers::pi_v<float> * freq / static_cast<float>(sampleRate)));
-    const float q = 0.5f;
-
-    m_hpL = left - m_lpL - q * m_bpL;
-    m_bpL += f * m_hpL;
-    m_lpL += f * m_bpL;
-    left = m_hpL;
-
-    m_hpR = right - m_lpR - q * m_bpR;
-    m_bpR += f * m_hpR;
-    m_lpR += f * m_bpR;
-    right = m_hpR;
+    const double g = std::tan(std::numbers::pi * static_cast<double>(freq) / static_cast<double>(sampleRate));
+    const double k = 1.0; // Q = 1.0 / k
+    const double damping = 1.0 / (1.0 + g * (g + k));
+
+    // Left channel
+    {
+        const double hp = (static_cast<double>(left) - (g + k) * m_s1L - m_s2L) * damping;
+        const double v1 = g * hp;
+        const double v = v1 + m_s1L;
+        m_s1L = v1 + v;
+        const double v2 = g * v;
+        const double lp = v2 + m_s2L;
+        m_s2L = v2 + lp;
+        left = static_cast<float>(hp);
+    }
+
+    // Right channel
+    {
+        const double hp = (static_cast<double>(right) - (g + k) * m_s1R - m_s2R) * damping;
+        const double v1 = g * hp;
+        const double v = v1 + m_s1R;
+        m_s1R = v1 + v;
+        const double v2 = g * v;
+        const double lp = v2 + m_s2R;
+        m_s2R = v2 + lp;
+        right = static_cast<float>(hp);
+    }
+
+    // NaN protection
+    if (std::isnan(left) || std::isnan(right)) {
+        reset();
+        left = 0.0f;
+        right = 0.0f;
+    }
 }
 
 void HighPassFilterEffect::reset()
 {
-    m_lpL = m_hpL = m_bpL = 0.0f;
-    m_lpR = m_hpR = m_bpR = 0.0f;
+    m_s1L = m_s2L = 0.0;
+    m_s1R = m_s2R = 0.0;
 }
 
 } // namespace noteahead

src/domain/devices/high_pass_filter_effect.hpp

@@ -29,8 +29,8 @@ class HighPassFilterEffect : public Effect
 
 private:
     float m_cutoff { 0.0f };
-    float m_lpL { 0.0f }, m_hpL { 0.0f }, m_bpL { 0.0f };
-    float m_lpR { 0.0f }, m_hpR { 0.0f }, m_bpR { 0.0f };
+    double m_s1L { 0.0 }, m_s2L { 0.0 };
+    double m_s1R { 0.0 }, m_s2R { 0.0 };
 };
 
 } // namespace noteahead

src/domain/devices/low_pass_filter_effect.cpp

@@ -32,25 +32,48 @@ void LowPassFilterEffect::process(float & left, float & right, uint32_t sampleRa
         return;
     }
 
+    // Zero-Delay Feedback State Variable Filter (2nd order)
     const float freq = 20.0f * std::pow(std::min(20000.0f, sampleRate * 0.49f) / 20.0f, m_cutoff);
-    const float f = 2.0f * std::sin(std::numbers::pi_v<float> * freq / static_cast<float>(sampleRate));
-    const float q = 0.5f;
-
-    m_hpL = left - m_lpL - q * m_bpL;
-    m_bpL += f * m_hpL;
-    m_lpL += f * m_bpL;
-    left = m_lpL;
-
-    m_hpR = right - m_lpR - q * m_bpR;
-    m_bpR += f * m_hpR;
-    m_lpR += f * m_bpR;
-    right = m_lpR;
+    const double g = std::tan(std::numbers::pi * static_cast<double>(freq) / static_cast<double>(sampleRate));
+    const double k = 1.0; // Q = 1.0 / k
+    const double damping = 1.0 / (1.0 + g * (g + k));
+
+    // Left channel
+    {
+        const double hp = (static_cast<double>(left) - (g + k) * m_s1L - m_s2L) * damping;
+        const double v1 = g * hp;
+        const double v = v1 + m_s1L;
+        m_s1L = v1 + v;
+        const double v2 = g * v;
+        const double lp = v2 + m_s2L;
+        m_s2L = v2 + lp;
+        left = static_cast<float>(lp);
+    }
+
+    // Right channel
+    {
+        const double hp = (static_cast<double>(right) - (g + k) * m_s1R - m_s2R) * damping;
+        const double v1 = g * hp;
+        const double v = v1 + m_s1R;
+        m_s1R = v1 + v;
+        const double v2 = g * v;
+        const double lp = v2 + m_s2R;
+        m_s2R = v2 + lp;
+        right = static_cast<float>(lp);
+    }
+
+    // NaN protection
+    if (std::isnan(left) || std::isnan(right)) {
+        reset();
+        left = 0.0f;
+        right = 0.0f;
+    }
 }
 
 void LowPassFilterEffect::reset()
 {
-    m_lpL = m_hpL = m_bpL = 0.0f;
-    m_lpR = m_hpR = m_bpR = 0.0f;
+    m_s1L = m_s2L = 0.0;
+    m_s1R = m_s2R = 0.0;
 }
 
 } // namespace noteahead

src/domain/devices/low_pass_filter_effect.hpp

@@ -29,8 +29,8 @@ class LowPassFilterEffect : public Effect
 
 private:
     float m_cutoff { 1.0f };
-    float m_lpL { 0.0f }, m_hpL { 0.0f }, m_bpL { 0.0f };
-    float m_lpR { 0.0f }, m_hpR { 0.0f }, m_bpR { 0.0f };
+    double m_s1L { 0.0 }, m_s2L { 0.0 };
+    double m_s1R { 0.0 }, m_s2R { 0.0 };
 };
 
 } // namespace noteahead

src/unit_tests/effects_test/effects_test.cpp

@@ -139,6 +139,27 @@ void EffectsTest::test_highPassFilterEffect()
     }
 }
 
+void EffectsTest::test_filterStability()
+{
+    LowPassFilterEffect lp;
+    HighPassFilterEffect hp;
+
+    for (int i = 0; i < 1000; ++i) {
+        float left = 1.0f;
+        float right = 1.0f;
+        const float cutoff = 0.5f + 0.49f * std::sin(i * 0.1f);
+
+        lp.setCutoff(cutoff);
+        hp.setCutoff(cutoff);
+
+        lp.process(left, right, 44100);
+        hp.process(left, right, 44100);
+
+        QVERIFY(!std::isnan(left));
+        QVERIFY(!std::isnan(right));
+    }
+}
+
 } // namespace noteahead
 
 QTEST_GUILESS_MAIN(noteahead::EffectsTest)

src/unit_tests/effects_test/effects_test.hpp

@@ -29,6 +29,7 @@ private slots:
     void test_panningEffect();
     void test_lowPassFilterEffect();
     void test_highPassFilterEffect();
+    void test_filterStability();
 };
 
 } // namespace noteahead