Add explicit computation of third-order bilinear transform

Chowdhury-DSP · Jan 13, 2024 · c7a38b2 · c7a38b2
1 parent 4355c19
commit c7a38b2
Show file tree

Hide file tree

Showing 2 changed files with 51 additions and 8 deletions.
diff --git a/modules/dsp/chowdsp_filters/Utils/chowdsp_ConformalMaps.h b/modules/dsp/chowdsp_filters/Utils/chowdsp_ConformalMaps.h
@@ -132,14 +132,20 @@ namespace ConformalMaps
         {
             // the mobius transform adds a little computational overhead, so here's the optimized bilinear transform.
             const auto KSq = K * K;
-            const auto a0_inv = (T) 1 / (as[0] * KSq + as[1] * K + as[2]);
+
+            const auto as0_KSq = as[0] * KSq;
+            const auto as1_K = as[1] * K;
+            const auto bs0_KSq = bs[0] * KSq;
+            const auto bs1_K = bs[1] * K;
+
+            const auto a0_inv = (T) 1 / (as0_KSq + as1_K + as[2]);
 
             a[0] = (T) 1;
-            a[1] = (T) 2 * (as[2] - as[0] * KSq) * a0_inv;
-            a[2] = (as[0] * KSq - as[1] * K + as[2]) * a0_inv;
-            b[0] = (bs[0] * KSq + bs[1] * K + bs[2]) * a0_inv;
-            b[1] = (T) 2 * (bs[2] - bs[0] * KSq) * a0_inv;
-            b[2] = (bs[0] * KSq - bs[1] * K + bs[2]) * a0_inv;
+            a[1] = (T) 2 * (as[2] - as0_KSq) * a0_inv;
+            a[2] = (as0_KSq - as1_K + as[2]) * a0_inv;
+            b[0] = (bs0_KSq + bs1_K + bs[2]) * a0_inv;
+            b[1] = (T) 2 * (bs[2] - bs0_KSq) * a0_inv;
+            b[2] = (bs0_KSq - bs1_K + bs[2]) * a0_inv;
         }
 
         /** Second-order alpha transform */
@@ -151,6 +157,37 @@ namespace ConformalMaps
         }
     };
 
+    /** Transforms for a third-order filter */
+    template <typename T>
+    struct Transform<T, 3>
+    {
+        /** Third-order bilinear transform */
+        static inline void bilinear (T (&b)[4], T (&a)[4], const T (&bs)[4], const T (&as)[4], T K)
+        {
+            // the mobius transform adds a little computational overhead, so here's the optimized bilinear transform.
+            const auto KSq = K * K;
+            const auto KCb = K * KSq;
+
+            const auto as0_KCb = as[0] * KCb;
+            const auto as1_KSq = as[1] * KSq;
+            const auto as2_K = as[2] * K;
+            const auto bs0_KCb = bs[0] * KCb;
+            const auto bs1_KSq = bs[1] * KSq;
+            const auto bs2_K = bs[2] * K;
+
+            const auto a0_inv = (T) 1 / (as0_KCb + as1_KSq + as2_K + as[3]);
+
+            a[0] = (T) 1;
+            a[1] = ((T) -3 * as0_KCb - as1_KSq + as2_K + (T) 3 * as[3]) * a0_inv;
+            a[2] = ((T) 3 * as0_KCb - as1_KSq - as2_K + (T) 3 * as[3]) * a0_inv;
+            a[3] = (-as0_KCb + as1_KSq - as2_K + as[3]) * a0_inv;
+            b[0] = (bs0_KCb + bs1_KSq + bs[2] * K + bs[3]) * a0_inv;
+            b[1] = ((T) -3 * bs0_KCb - bs1_KSq + bs2_K + (T) 3 * bs[3]) * a0_inv;
+            b[2] = ((T) 3 * bs0_KCb - bs1_KSq - bs2_K + (T) 3 * bs[3]) * a0_inv;
+            b[3] = (-bs0_KCb + bs1_KSq - bs2_K + bs[3]) * a0_inv;
+        }
+    };
+
     /** Computes the warping factor "K" so that the angular frequency wc is matched at sample rate fs */
     template <typename T, typename NumericType>
     inline T computeKValueAngular (T wc, NumericType fs)

diff --git a/tests/dsp_tests/chowdsp_filters_test/ConformalMapsTest.cpp b/tests/dsp_tests/chowdsp_filters_test/ConformalMapsTest.cpp
@@ -13,8 +13,8 @@ void testBilinearTransform (const float (&bs)[N], const float (&as)[N], const fl
 
     for (int i = 0; i < N; ++i)
     {
-        REQUIRE_MESSAGE (b[i] == Catch::Approx (b_exp[i]).margin (1.0e-6f), "B coefficient " << std::to_string (i) << " is incorrect!");
-        REQUIRE_MESSAGE (a[i] == Catch::Approx (a_exp[i]).margin (1.0e-6f), "A coefficient " << std::to_string (i) << " is incorrect!");
+        REQUIRE_MESSAGE (b[i] == Catch::Approx (b_exp[i]).margin (1.0e-3f), "B coefficient " << std::to_string (i) << " is incorrect!");
+        REQUIRE_MESSAGE (a[i] == Catch::Approx (a_exp[i]).margin (1.0e-3f), "A coefficient " << std::to_string (i) << " is incorrect!");
     }
 }
 
@@ -35,6 +35,12 @@ TEST_CASE ("Conformal Maps Test", "[dsp][filters]")
         testBilinearTransform ({ 1.0f, 2.0f, 3.0 }, { 2.0f, 3.0f, 4.0f }, { 0.50000304f, -0.999981701f, 0.499978632f }, { 1.f, -1.9999634f, 0.999963343f }, kVal);
     }
 
+    SECTION ("Third-Order Bilinear Transform")
+    {
+        testBilinearTransform ({ 1.0f, 2.0f, 3.0f, 4.0f }, { 2.0f, 3.0f, 4.0f, 5.0f }, { 0.499997396f, -1.50001302f, 1.50003385f, -0.500018230f }, { 1.0f, -3.00003125f, 3.00006250f, -1.00003125f }, 2 * fs);
+        testBilinearTransform ({ 1.0f, 2.0f, 3.0f, 4.0f }, { 2.0f, 3.0f, 4.0f, 5.0f }, { 0.499996947f, -1.50001527f, 1.50003969f, -0.500021372f }, { 1.0f, -3.00003664f, 3.00007327f, -1.00003664f }, kVal);
+    }
+
     SECTION ("Fourth-Order Bilinear Transform")
     {
         testBilinearTransform ({ 1.0f, 2.0f, 3.0f, 4.0f, 5.0f }, { 2.0f, 3.0f, 4.0f, 5.0f, 6.0f }, { 0.500002623f, -1.99998963f, 2.99995327f, -1.99994779f, 0.499981761f }, { 1.f, -3.99996853f, 5.99990654f, -3.9999063f, 0.999968707f }, 2 * fs);