From ba328ae9cc2316d7716419e5217c619afebf58c2 Mon Sep 17 00:00:00 2001 From: Andrew Belt Date: Tue, 6 Jan 2026 01:13:08 -0500 Subject: [PATCH] VCF: Optimize function approximations. Increase maximum frequency to 22,000 Hz. --- src/VCF.cpp | 67 ++++++++++++++++++++++++++++------------------------- 1 file changed, 36 insertions(+), 31 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index 5db0615..932a282 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -1,53 +1,56 @@ #include "plugin.hpp" -using simd::float_4; - - -/** Approximates tan(x) for x in [0, π*0.49]. -Max error: 1.8e-7 +/** Approximates tan(pi*x) for x in [0, 0.5). +Optimized coefficients for max rel error: 1.18e-05. */ template -inline T tan_5_4(T x) { +inline T tan_pi_1_2(T x) { T x2 = x * x; - return x * (T(1) + x2 * (T(-0.1122835153) + T(0.00114264086) * x2)) - / (T(1) + x2 * (T(-0.4456155807) + T(0.01634547741) * x2)); + T num = T(1) + T(-0.9622845476351338) * x2; + T den = T(1) + x2 * (T(-4.252711329946427) + T(1.0108965067534537) * x2); + return T(M_PI) * x * num / den; } /** Approximates tanh(x) for x in [-4, 4]. -Max error: 1.1e-5 +Optimized coefficients for max rel error: 5.56e-07. */ template -inline T tanh_5_6(T x) { +inline T tanh_2_3(T x) { x = simd::clamp(x, T(-4), T(4)); T x2 = x * x; - return x * (T(1450.73) + x2 * (T(152.494) + T(1.12357) * x2)) - / (T(1450.80) + x2 * (T(635.839) + x2 * (T(19.8913) + T(0.00188) * x2))); + T num = T(1) + x2 * (T(0.11823399425250458) + x2 * T(0.0017593473306865004)); + T den = T(1) + x2 * (T(0.4515649138214517) + x2 * (T(0.01895237471013944) + x2 * T(7.340776670083926e-05))); + return x * num / den; } -/** Approximates ln(cosh(x)) for x >= 0. -Max error: 6.2e-5 over [0, 10]. +/** Approximates ln(cosh(x)) for all x. +Optimized coefficients for max deriv rel error: 2.18e-04. +Max func abs error: 2.33e-04. */ template -inline T lncosh_5_3(T x) { +inline T ln_cosh_3_3(T x) { + x = simd::abs(x); T x2 = x * x; - return x2 * T(0.5) * (T(1) + x * (T(0.5825995649809) + x * (T(0.29698351914621) + x * T(-0.0001879624507617)))) - / (T(1) + x * (T(0.5938203513329) + x * (T(0.4266725633273) + x * T(0.1456088792904)))); + T num = T(1) + x * (T(0.5536821172234367) + x * (T(0.25181887093756017) + x * T(-0.00022483065845568806))); + T den = T(1) + x * (T(0.5571167780242312) + x * (T(0.4011383699340152) + x * T(0.12250932588907415))); + return x2 * T(0.5) * num / den; } -/** First antiderivative of tanh: F(x) = ln(cosh(x)). -For numerical stability, uses asymptotic form for |x| > 10. +/** Approximates ln(cosh(x)) for all x. +Optimized coefficients for max deriv rel error: 1.73e-05. +Max func abs error: 1.59e-05. */ template -inline T tanh_F1(T x) { +inline T ln_cosh_3_4(T x) { x = simd::abs(x); - // For |x| > 10, ln(cosh(x)) ≈ |x| - ln(2) - T asymptotic = x - T(0.6931471805599453); - T exact = lncosh_5_3(x); - return simd::ifelse(x < T(10), exact, asymptotic); + T x2 = x * x; + T num = T(1) + x * (T(0.6818733052030992) + x * (T(0.3310714761981197) + x * T(0.07328250587472884))); + T den = T(1) + x * (T(0.6818667996836973) + x * (T(0.4970931090131887) + x * (T(0.18905704720922234) + x * T(0.0366985697841935)))); + return x2 * T(0.5) * num / den; } @@ -68,13 +71,13 @@ struct TanhADAA1 { T diff = x - xPrev; // Normal case: finite difference of antiderivative - T Fx = tanh_F1(x); - T FxPrev = tanh_F1(xPrev); + T Fx = ln_cosh_3_4(x); + T FxPrev = ln_cosh_3_4(xPrev); T adaaResult = (Fx - FxPrev) / diff; // Fallback: use tanh at midpoint when diff is small T midpoint = (x + xPrev) * T(0.5); - T fallbackResult = tanh_5_6(midpoint); + T fallbackResult = tanh_2_3(midpoint); T y = simd::ifelse(simd::abs(diff) > eps, adaaResult, fallbackResult); @@ -120,7 +123,7 @@ struct LadderFilter { void process(Frame& frame) { // Pre-warp cutoff frequency using bilinear transform - T g = tan_5_4(T(M_PI) * frame.cutoff); + T g = tan_pi_1_2(frame.cutoff); // Compute the one-pole lowpass gain coefficient G = g/(1+g) T G = g / (T(1) + g); @@ -156,12 +159,14 @@ struct LadderFilter { frame.lowpass4 = y3; // Binomial expansion (1-H)^4 = 1 - 4H + 6H^2 - 4H^3 + H^4 - // Using saturated input and stage outputs frame.highpass4 = sat0 - T(4) * y0 + T(6) * y1 - T(4) * y2 + y3; } }; +using simd::float_4; + + struct VCF : Module { enum ParamIds { FREQ_PARAM, @@ -196,7 +201,7 @@ struct VCF : Module { // or // param = (log2(freq / C4) + 5) / 10 const float minFreq = (std::log2(8.f / dsp::FREQ_C4) + 5) / 10; - const float maxFreq = (std::log2(20000.f / dsp::FREQ_C4) + 5) / 10; + const float maxFreq = (std::log2(22000.f / dsp::FREQ_C4) + 5) / 10; const float defaultFreq = (minFreq + maxFreq) / 2.f; configParam(FREQ_PARAM, minFreq, maxFreq, defaultFreq, "Cutoff frequency", " Hz", std::pow(2, 10.f), dsp::FREQ_C4 / std::pow(2, 5.f)); configParam(RES_PARAM, 0.f, 1.f, 0.f, "Resonance", "%", 0.f, 100.f); @@ -267,7 +272,7 @@ struct VCF : Module { // Cutoff frequency float_4 pitch = freqParam + inputs[FREQ_INPUT].getPolyVoltageSimd(c) * freqCvParam; float_4 cutoff = dsp::FREQ_C4 * dsp::exp2_taylor5(pitch); - frame.cutoff = simd::clamp(cutoff * args.sampleTime, 0.f, 0.49f); + frame.cutoff = simd::clamp(cutoff * args.sampleTime, 0.f, 0.499f); // Process filters[c / 4].process(frame);