diff --git a/CMakeLists.txt b/CMakeLists.txt index 2f45f00..ecfed7e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -995,6 +995,18 @@ target_link_libraries(curve_popup PUBLIC editor_geometry) add_library(master_gain STATIC src/core/instrument/engine/master_gain.cpp) target_include_directories(master_gain PUBLIC src) +# filter — the per-voice 2-pole resonant low/high-pass, ported from Daniel's Cortex-M4 filter +# with its virtual FilterBase/Filter/Biquad hierarchy flattened away (that hierarchy dispatched +# virtually per channel per sample, which the per-voice per-sample path forbids). Control +# mapping, RBJ coefficient math, feedback saturation, and the filter type each get their own +# file; VoiceFilter::process is header-inline so the biquad kernel still inlines at the call +# site. Standard library only. NEITHER SDK. +add_library(filter STATIC + src/core/instrument/engine/filter/filter_params.cpp + src/core/instrument/engine/filter/filter_coeffs.cpp + src/core/instrument/engine/filter/voice_filter.cpp) +target_include_directories(filter PUBLIC src) + # sample_bands: the band-stack allocator's vertical inventory, asserted as pure geometry # (chrome / two-lane waveform / deck row) independent of any paint call — the contract the # band owners downstream read. @@ -1095,6 +1107,13 @@ add_executable(master_gain_tests tests/test_master_gain.cpp) target_link_libraries(master_gain_tests PRIVATE master_gain) add_test(NAME master_gain_tests COMMAND master_gain_tests) +# filter: the per-voice resonant filter. Pins the RBJ coefficients against an independent +# textbook cos/sin derivation, asserts the cutoff/Q control mappings at their anchors, and +# measures the resonant peak both analytically and by driving real sines. NEITHER SDK. +add_executable(filter_tests tests/test_filter.cpp) +target_link_libraries(filter_tests PRIVATE filter) +add_test(NAME filter_tests COMMAND filter_tests) + # --------------------------------------------------------------------------- # 4) The REAPER extension — a loadable module (dlopen'd by REAPER, not linked). # --------------------------------------------------------------------------- diff --git a/src/core/instrument/engine/filter/CLAUDE.md b/src/core/instrument/engine/filter/CLAUDE.md new file mode 100644 index 0000000..ecec0f7 --- /dev/null +++ b/src/core/instrument/engine/filter/CLAUDE.md @@ -0,0 +1,92 @@ +# src/core/instrument/engine/filter — the per-voice resonant filter + +## Scope + +The pure 2-pole resonant low/high-pass a sounding voice runs. No REAPER, no VST3, no +allocation, no I/O. Four files, one responsibility each: + +- `filter_params` — the control domain: `FilterMode`, normalized [0,1] knob position → + cutoff Hz and Q, and the exact inverses. +- `filter_coeffs` — the DSP domain: `BiquadCoeffs` and the RBJ coefficient computation + from (mode, cutoff Hz, Q, sample rate). +- `filter_saturate` — the high-pass feedback saturator (`tanhSaturate` / + `feedbackSaturate`). Header-only inline; it sits on the per-sample path. +- `voice_filter` — `FilterSettings` and `VoiceFilter`, the concrete per-voice type. + `process()` is defined in the header. + +## Invariants + +### No vtable on the per-sample path + +This is a **port, not a relocation**. The Cortex-M4 source was a virtual hierarchy +(`FilterBase` → `Filter` → `Biquad` → `{BiquadHP, BiquadLP}`) whose base class routed the +channel loop through pure-virtual `process_channel_frame` / `filter` / `update_feedback` +so a `FilterDecorator` chain could wrap it. **None of that came across, and none of it may +come back.** `VoiceFilter` is concrete: mode is a member branch inside an inlined +`process()`, predicted perfectly because it cannot change within a note. There is no +`IFilter`, no decorator seam, no virtual `tick()`, and no allocation in `process()` — root +`CLAUDE.md`'s structural heuristic 3 names this class of dispatch blowout directly. + +A non-type template parameter for the mode was considered and rejected: mode is a +runtime-settable user parameter, so templating would only relocate the same branch to the +call site and force the voice to hold two instances or switch over them. + +### Two modes, and only two + +2-pole high-pass and 2-pole low-pass. The source's `Biquad1PoleLP` is struck and was not +ported. Further modes are deferred — **do not build a mode-extension framework** for them. + +### The cutoff control is sample-rate-free; the clamp is not + +`filterCutoffHzFromNorm` sweeps a fixed 20 Hz – 20 kHz (three exact decades, so norm 1/3 +is 200 Hz and 2/3 is 2 kHz) and takes no sample rate. The persisted value is the +normalized knob position, so a rate-derived endpoint would make one preset sound different +at 44.1k and 96k. The Nyquist clamp (`kFilterNyquistFraction`, 0.48) is a property of the +bilinear transform — `tan(pi*fc/sr)` diverges at Nyquist — so it lives in `biquadCoeffs` +where the rate is already a parameter. 20 kHz is under 0.48·sr at every supported rate, so +the clamp never eats live knob travel; the source's hardcoded 23 kHz endpoint did exactly +that at 44.1k. + +`biquadCoeffs` with a non-positive sample rate returns pass-through coefficients. It does +**not** fall back to 44100 — that would breach the standing no-hardcoded-sample-rates +ruling. + +### Q spans 0.1 → 10 with √2 at the center + +Settled by Daniel. The source's `Q = M_SQRT1_2 + resonance` mapping (floored at 0.707, no +center anchor) was **rewritten, not ported**. The curve is quadratic in log Q through the +three anchors rather than two spliced log segments — same anchors either way, but no slope +kink at the center detent. The quadratic term is nonzero only because √2 is not the +geometric mean of 0.1 and 10; `filterNormFromQ` divides by it. + +### The high-pass input feedback is load-bearing + +`kHighPassFeedbackShare` (0.24) times the raw **normalized** resonance, not Q — Q reaches +10 and scaling the feedback by it would push loop gain past unity. The high-pass numerator +collapses toward zero as cutoff falls, taking the resonance with it; the saturated +feedback restores the character down there. Ported behavior; the constant is the tuning +knob if the feel needs adjusting. `audio_saturate` and `H()` from the source were unused +by the biquads and were not ported. + +### Denormal flushing + +`process()` flushes the **y** history to exact zero below `kFilterDenormalFloor` (1e-30). +Only the recursive half needs it: a denormal in `y` self-sustains and stalls the FPU for +thousands of samples on a ringing-out voice, while the `x` history is an FIR tail that +shifts out within two samples. `isSilent()` reports the flushed state and is the honest +signal that a voice's filter can no longer contribute output. + +## Gotchas + +- **The tan pre-warp is not a different filter.** By the half-angle identity + `cos(w0) = (1-w²)/(1+w²)` and `sin(w0) = 2w/(1+w²)` with `w = tan(pi*fc/sr)`, these are + the textbook RBJ cos/sin coefficients exactly — just computed in a form that stays + conditioned at low cutoff where `cos(w0) → 1`. `tests/test_filter.cpp` asserts the + equivalence against an independent derivation. Don't "simplify" it back to `std::cos`. +- **`prepare()` deliberately does not clear history** — a live parameter move must glide, + not click. Call `reset()` at note-on. +- **`a1`/`a2` are stored for a subtracting difference equation** (`y = ... - a1*y1 - + a2*y2`), so the transfer denominator is `1 + a1*z^-1 + a2*z^-2`. A sign convention slip + here inverts the poles. +- **No call site yet.** Wiring the filter into the voice path is a separate track; nothing + in `sampler_core` references this module today. diff --git a/src/core/instrument/engine/filter/filter_coeffs.cpp b/src/core/instrument/engine/filter/filter_coeffs.cpp new file mode 100644 index 0000000..2cd973d --- /dev/null +++ b/src/core/instrument/engine/filter/filter_coeffs.cpp @@ -0,0 +1,43 @@ +#include "core/instrument/engine/filter/filter_coeffs.h" + +#include + +namespace reasampler::instrument::engine { +namespace { + +// M_PI is not standard C++ and is absent on MSVC without _USE_MATH_DEFINES. +constexpr double kPi = 3.14159265358979323846; + +double clampd(double v, double lo, double hi) { return v < lo ? lo : (v > hi ? hi : v); } + +} // namespace + +BiquadCoeffs biquadCoeffs(FilterMode mode, float cutoffHz, float q, double sampleRate) { + if (!(sampleRate > 0.0)) return BiquadCoeffs{}; + + const double nyquistCeiling = kFilterNyquistFraction * sampleRate; + const double fc = clampd(cutoffHz, kFilterCutoffMinHz, nyquistCeiling); + const double qq = clampd(q, kFilterQMin, kFilterQMax); + + const double w = std::tan(kPi * fc / sampleRate); + const double w2 = w * w; + const double cosw = (1.0 - w2) / (1.0 + w2); + const double sinw = 2.0 * w / (1.0 + w2); + const double alpha = sinw / (2.0 * qq); + const double norm = 1.0 / (1.0 + alpha); + + // Both modes share the denominator; only the numerator's sign on cosw differs, and b1 is + // always +/-2*b0 — folding that in keeps the two branches from drifting apart. + const double b0 = (mode == FilterMode::HighPass ? (1.0 + cosw) : (1.0 - cosw)) * 0.5 * norm; + const double b1 = (mode == FilterMode::HighPass ? -2.0 : 2.0) * b0; + + BiquadCoeffs c; + c.b0 = static_cast(b0); + c.b1 = static_cast(b1); + c.b2 = static_cast(b0); + c.a1 = static_cast(-2.0 * cosw * norm); + c.a2 = static_cast((1.0 - alpha) * norm); + return c; +} + +} // namespace reasampler::instrument::engine diff --git a/src/core/instrument/engine/filter/filter_coeffs.h b/src/core/instrument/engine/filter/filter_coeffs.h new file mode 100644 index 0000000..67699e2 --- /dev/null +++ b/src/core/instrument/engine/filter/filter_coeffs.h @@ -0,0 +1,31 @@ +// filter_coeffs.h — RBJ Audio EQ Cookbook Direct Form I biquad coefficients for the 2-pole +// low/high-pass. Computed via the tan half-angle substitution w = tan(pi*fc/sr): by the +// identity cos(w0) = (1-w^2)/(1+w^2), sin(w0) = 2w/(1+w^2) these ARE the textbook cos/sin +// coefficients, in a form that stays conditioned at low cutoff where cos(w0) -> 1. + +#pragma once + +#include "core/instrument/engine/filter/filter_params.h" + +namespace reasampler::instrument::engine { + +// Already normalized by a0. The denominator is 1 + a1*z^-1 + a2*z^-2, so the difference +// equation SUBTRACTS the a terms: y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2. +struct BiquadCoeffs { + float b0 = 1.0f; + float b1 = 0.0f; + float b2 = 0.0f; + float a1 = 0.0f; + float a2 = 0.0f; +}; + +// Highest fraction of the sample rate the pre-warp stays well-conditioned at: tan() diverges +// as fc approaches sr/2. Ported unchanged from the firmware, where it was already the ceiling. +inline constexpr double kFilterNyquistFraction = 0.48; + +// cutoffHz is clamped into [kFilterCutoffMinHz, kFilterNyquistFraction*sampleRate] and q into +// [kFilterQMin, kFilterQMax]. A non-positive sampleRate yields pass-through coefficients — the +// no-hardcoded-sample-rates ruling means we refuse to invent a rate rather than assume 44.1k. +BiquadCoeffs biquadCoeffs(FilterMode mode, float cutoffHz, float q, double sampleRate); + +} // namespace reasampler::instrument::engine diff --git a/src/core/instrument/engine/filter/filter_params.cpp b/src/core/instrument/engine/filter/filter_params.cpp new file mode 100644 index 0000000..7ac5bbf --- /dev/null +++ b/src/core/instrument/engine/filter/filter_params.cpp @@ -0,0 +1,59 @@ +#include "core/instrument/engine/filter/filter_params.h" + +#include + +namespace reasampler::instrument::engine { +namespace { + +double clamp01(double v) { return v < 0.0 ? 0.0 : (v > 1.0 ? 1.0 : v); } + +// log Q = A + B*n + C*n^2, solved from the three anchor points. C is nonzero precisely +// because the center anchor sqrt(2) is not the geometric mean of the endpoints (which is 1); +// were they equal the curve would degenerate to a plain log sweep and the inverse below +// would divide by zero. +struct QCurve { + double a, b, c; +}; + +QCurve qCurve() { + const double lo = std::log(static_cast(kFilterQMin)); + const double mid = std::log(static_cast(kFilterQCenter)); + const double hi = std::log(static_cast(kFilterQMax)); + return {lo, 4.0 * mid - 3.0 * lo - hi, 2.0 * lo + 2.0 * hi - 4.0 * mid}; +} + +} // namespace + +float filterCutoffHzFromNorm(float norm) { + const double lo = std::log(static_cast(kFilterCutoffMinHz)); + const double hi = std::log(static_cast(kFilterCutoffMaxHz)); + return static_cast(std::exp(lo + clamp01(norm) * (hi - lo))); +} + +float filterNormFromCutoffHz(float hz) { + if (!(hz > 0.0f)) return 0.0f; + const double lo = std::log(static_cast(kFilterCutoffMinHz)); + const double hi = std::log(static_cast(kFilterCutoffMaxHz)); + return static_cast(clamp01((std::log(static_cast(hz)) - lo) / (hi - lo))); +} + +float filterQFromNorm(float norm) { + const QCurve k = qCurve(); + const double n = clamp01(norm); + return static_cast(std::exp(k.a + n * (k.b + k.c * n))); +} + +float filterNormFromQ(float q) { + if (!(q > kFilterQMin)) return 0.0f; + if (q >= kFilterQMax) return 1.0f; + // Clamping first is load-bearing, not just tidy: the parabola peaks at log Q well below + // an arbitrarily large q, so an unclamped out-of-range value has no real root at all. + const QCurve k = qCurve(); + const double d = k.b * k.b - 4.0 * k.c * (k.a - std::log(static_cast(q))); + if (!(d >= 0.0)) return 0.0f; + // Of the two roots only this one lies on the rising branch inside [0,1]; the parabola's + // vertex sits well above 1 for the settled anchors. + return static_cast(clamp01((-k.b + std::sqrt(d)) / (2.0 * k.c))); +} + +} // namespace reasampler::instrument::engine diff --git a/src/core/instrument/engine/filter/filter_params.h b/src/core/instrument/engine/filter/filter_params.h new file mode 100644 index 0000000..9b642c5 --- /dev/null +++ b/src/core/instrument/engine/filter/filter_params.h @@ -0,0 +1,39 @@ +// filter_params.h — control-domain mapping for the voice filter: normalized [0,1] knob +// positions to cutoff Hz and Q, plus the two-mode enum. Deliberately sample-rate-free — +// the Nyquist clamp is a property of the bilinear transform and lives in filter_coeffs, +// so the persisted normalized cutoff means the same frequency at every project rate. + +#pragma once + +namespace reasampler::instrument::engine { + +enum class FilterMode { LowPass, HighPass }; + +// The audio band the cutoff control sweeps: three exact decades, so norm 1/3 is 200 Hz and +// norm 2/3 is 2 kHz. NOT derived from the sample rate — a rate-dependent endpoint would make +// one saved preset sound different at 44.1k and 96k, and at 44.1k the top of the travel would +// be dead against the Nyquist clamp (the ported firmware's 23 kHz endpoint had exactly that +// defect). 20 kHz sits under 0.48*sr at every rate we support, so the whole knob stays live. +inline constexpr float kFilterCutoffMinHz = 20.0f; +inline constexpr float kFilterCutoffMaxHz = 20000.0f; + +// Q spans the full range with Butterworth (sqrt(2)) at the control's center detent. +inline constexpr float kFilterQMin = 0.1f; +inline constexpr float kFilterQMax = 10.0f; +inline constexpr float kFilterQCenter = 1.41421356f; + +// Out-of-range norm clamps to the endpoints. +float filterCutoffHzFromNorm(float norm); + +// Exact inverse of filterCutoffHzFromNorm over the band; out-of-band Hz clamps to 0 or 1. +float filterNormFromCutoffHz(float hz); + +// A single smooth curve — quadratic in log Q — through (0, kFilterQMin), +// (0.5, kFilterQCenter), (1, kFilterQMax), rather than two spliced log segments. Same three +// anchors either way, but the single curve has no slope kink at the center detent. +float filterQFromNorm(float norm); + +// Exact inverse of filterQFromNorm; out-of-range Q clamps to 0 or 1. +float filterNormFromQ(float q); + +} // namespace reasampler::instrument::engine diff --git a/src/core/instrument/engine/filter/filter_saturate.h b/src/core/instrument/engine/filter/filter_saturate.h new file mode 100644 index 0000000..1a20da2 --- /dev/null +++ b/src/core/instrument/engine/filter/filter_saturate.h @@ -0,0 +1,27 @@ +// filter_saturate.h — the high-pass feedback-path saturator, ported from the Cortex-M4 +// filter. Header-inline: it sits on the per-voice per-sample path, and a rational +// approximation is here precisely to avoid a transcendental tanh() call there. + +#pragma once + +namespace reasampler::instrument::engine { + +// Rational tanh approximation inside +/-threshold, continued past it with a gentle 0.1 slope +// anchored at the threshold value so the curve stays continuous rather than hard-clipping. +inline float tanhSaturate(float x, float threshold, float a, float b) { + if (x > threshold) { + const float satAtThreshold = threshold * a / (a + b + threshold * threshold); + return satAtThreshold + (x - threshold) * 0.1f; + } + if (x < -threshold) { + const float satAtThreshold = -threshold * a / (a + b + threshold * threshold); + return satAtThreshold + (x + threshold) * 0.1f; + } + return x * a / (a + b + x * x); +} + +// TB-303-style hard feedback saturation. Tuned for the large excursions a resonant feedback +// path produces, not for audio-level signals — do not reuse it as a general waveshaper. +inline float feedbackSaturate(float x) { return tanhSaturate(x, 2.0f, 27.0f, 9.0f); } + +} // namespace reasampler::instrument::engine diff --git a/src/core/instrument/engine/filter/voice_filter.cpp b/src/core/instrument/engine/filter/voice_filter.cpp new file mode 100644 index 0000000..d926973 --- /dev/null +++ b/src/core/instrument/engine/filter/voice_filter.cpp @@ -0,0 +1,29 @@ +#include "core/instrument/engine/filter/voice_filter.h" + +namespace reasampler::instrument::engine { + +void VoiceFilter::prepare(const FilterSettings& settings, double sampleRate) { + mode_ = settings.mode; + const float cutoffHz = filterCutoffHzFromNorm(settings.cutoffNorm); + coeffs_ = biquadCoeffs(settings.mode, cutoffHz, filterQFromNorm(settings.resonanceNorm), + sampleRate); + const float res = settings.resonanceNorm < 0.0f + ? 0.0f + : (settings.resonanceNorm > 1.0f ? 1.0f : settings.resonanceNorm); + fbAmount_ = res * kHighPassFeedbackShare; +} + +void VoiceFilter::reset() { + for (State& s : state_) s = State{}; +} + +bool VoiceFilter::isSilent() const { + for (const State& s : state_) { + if (s.x1 != 0.0f || s.x2 != 0.0f || s.y1 != 0.0f || s.y2 != 0.0f || s.fb != 0.0f) { + return false; + } + } + return true; +} + +} // namespace reasampler::instrument::engine diff --git a/src/core/instrument/engine/filter/voice_filter.h b/src/core/instrument/engine/filter/voice_filter.h new file mode 100644 index 0000000..5e200ef --- /dev/null +++ b/src/core/instrument/engine/filter/voice_filter.h @@ -0,0 +1,113 @@ +// voice_filter.h — per-voice 2-pole resonant low/high-pass. Concrete type, no vtable: this +// sits on the per-voice per-sample path, so process() is header-inline and mode is a member +// branch. No allocation, no virtual dispatch, no I/O anywhere in process(). + +#pragma once + +#include +#include + +#include "core/instrument/engine/filter/filter_coeffs.h" +#include "core/instrument/engine/filter/filter_params.h" +#include "core/instrument/engine/filter/filter_saturate.h" + +namespace reasampler::instrument::engine { + +// Normalized control positions, as the editor moves them and the persisted state carries them. +struct FilterSettings { + FilterMode mode = FilterMode::LowPass; + float cutoffNorm = 1.0f; + float resonanceNorm = 0.0f; +}; + +// Share of the last output fed back into the high-pass input at full resonance. Driven by the +// raw control position rather than by Q: Q reaches 10, and scaling the feedback by it would +// push the loop gain past unity at the top of the range. +inline constexpr float kHighPassFeedbackShare = 0.24f; + +// Below this the recursion has decayed past -600 dB. Flushing keeps the history out of the +// subnormal range, where a ringing-out voice would otherwise stall the FPU for thousands of +// samples. Chosen well above FLT_MIN so a flushed state can never re-enter that range. +inline constexpr float kFilterDenormalFloor = 1e-30f; + +class VoiceFilter { +public: + // The instrument's output bus is permanently stereo; one history line per channel. + static constexpr int kMaxChannels = 2; + + struct State { + float x1 = 0.0f; + float x2 = 0.0f; + float y1 = 0.0f; + float y2 = 0.0f; + float fb = 0.0f; // last output; the high-pass input-feedback tap + }; + + // Recomputes coefficients from the control positions. History is deliberately preserved so + // a live parameter move glides instead of clicking; call reset() at note-on. + void prepare(const FilterSettings& settings, double sampleRate); + + void reset(); + + // Hot path. `channel` must be in [0, kMaxChannels). + float process(int channel, float x) { + assert(channel >= 0 && channel < kMaxChannels); + State& s = state_[channel]; + + // The high-pass numerator collapses toward zero as cutoff falls, taking the resonance + // with it; feeding a saturated share of the last output back into the input restores + // the character the coefficients alone stop producing down there. + const float in = (mode_ == FilterMode::HighPass) + ? x - fbAmount_ * feedbackSaturate(s.fb * 0.9f) + : x; + + const float y = coeffs_.b0 * in + coeffs_.b1 * s.x1 + coeffs_.b2 * s.x2 + - coeffs_.a1 * s.y1 - coeffs_.a2 * s.y2; + + s.x2 = s.x1; + s.x1 = in; + s.y2 = s.y1; + s.y1 = y; + + // Snap the WHOLE state once the recursion as a whole has decayed past -600 dB. + // Zeroing individual samples instead does not work: a resonator swings through zero + // twice a cycle, so a per-sample flush injects a step in phase with the resonance, + // which the resonance then amplifies — the filter limit-cycles at the floor forever + // rather than going quiet. Testing y1 AND y2 tests the envelope, not one sample. + if (s.y1 > -kFilterDenormalFloor && s.y1 < kFilterDenormalFloor && + s.y2 > -kFilterDenormalFloor && s.y2 < kFilterDenormalFloor) { + s = State{}; + } + s.fb = s.y1; + return y; + } + + void processFrame(float* samples, int channelCount) { + assert(channelCount >= 0 && channelCount <= kMaxChannels); + for (int c = 0; c < channelCount; ++c) samples[c] = process(c, samples[c]); + } + + // True once every history line has flushed to exact zero — the voice's filter has stopped + // ringing and cannot contribute further output. + bool isSilent() const; + + const State& state(int channel) const { + assert(channel >= 0 && channel < kMaxChannels); + return state_[channel]; + } + const BiquadCoeffs& coeffs() const { return coeffs_; } + +private: + BiquadCoeffs coeffs_{}; + FilterMode mode_ = FilterMode::LowPass; + float fbAmount_ = 0.0f; + State state_[kMaxChannels]{}; +}; + +// The port's whole point, enforced by the compiler rather than by review: the source was a +// virtual hierarchy dispatching per channel per sample, and this type must never grow one +// back. Trivially copyable also means nothing here is heap-owned. +static_assert(!std::is_polymorphic_v, "no vtable on the per-sample path"); +static_assert(std::is_trivially_copyable_v, "state is plain values, never owned"); + +} // namespace reasampler::instrument::engine diff --git a/src/temp-cortex/biquad.hpp b/src/temp-cortex/biquad.hpp deleted file mode 100644 index b8c916e..0000000 --- a/src/temp-cortex/biquad.hpp +++ /dev/null @@ -1,73 +0,0 @@ -#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 deleted file mode 100644 index e193be0..0000000 --- a/src/temp-cortex/biquad.tpp +++ /dev/null @@ -1,195 +0,0 @@ -#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 deleted file mode 100644 index ccc3d0a..0000000 --- a/src/temp-cortex/filter.hpp +++ /dev/null @@ -1,62 +0,0 @@ -#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/filter_params.hpp b/src/temp-cortex/filter_params.hpp deleted file mode 100644 index b9f51c3..0000000 --- a/src/temp-cortex/filter_params.hpp +++ /dev/null @@ -1,22 +0,0 @@ -#pragma once - -typedef struct -{ - float cutoff; - float res; - float Q; -} FilterParameters; - -typedef struct { - float a1; - float a2; - float b0; - float b1; - float b2; -} NormalCoefficients; - -typedef struct { - float x[2]; // Previous inputs - float y[2]; // Previous outputs - float fb; // Feedback value for resonance -} FeedbackLine; \ No newline at end of file diff --git a/src/temp-cortex/util.hpp b/src/temp-cortex/util.hpp deleted file mode 100644 index 88e0b42..0000000 --- a/src/temp-cortex/util.hpp +++ /dev/null @@ -1,46 +0,0 @@ -#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 diff --git a/tests/test_filter.cpp b/tests/test_filter.cpp new file mode 100644 index 0000000..22e06c7 --- /dev/null +++ b/tests/test_filter.cpp @@ -0,0 +1,450 @@ +// Standalone tests for the per-voice filter — no VST3, no REAPER, no framework. Same fast +// assert loop as the sibling pure tests. The coefficient pins are literals so a refactor that +// changes the DSP fails loudly; they are cross-checked in-test against a textbook RBJ +// derivation (std::cos/std::sin) that shares no code with the implementation. + +#include "../src/core/instrument/engine/filter/filter_coeffs.h" +#include "../src/core/instrument/engine/filter/filter_params.h" +#include "../src/core/instrument/engine/filter/filter_saturate.h" +#include "../src/core/instrument/engine/filter/voice_filter.h" + +#include +#include +#include +#include + +using namespace reasampler::instrument::engine; + +static int g_fail = 0; +#define CHECK(cond) do { if(!(cond)) { \ + std::printf("FAIL line %d: %s\n", __LINE__, #cond); ++g_fail; } } while(0) +#define CHECK_NEAR(a, b, eps) do { const double a_ = (a), b_ = (b); \ + if (!(std::fabs(a_ - b_) <= (eps))) { \ + std::printf("FAIL line %d: %s (%.10f) != %s (%.10f), delta %.3e\n", \ + __LINE__, #a, a_, #b, b_, std::fabs(a_ - b_)); ++g_fail; } } while(0) + +static constexpr double kPi = 3.14159265358979323846; + +// --------------------------------------------------------------------------- +// Cutoff mapping +// --------------------------------------------------------------------------- + +static void testCutoffMapsThreeDecadesLogarithmically() { + CHECK_NEAR(filterCutoffHzFromNorm(0.0f), 20.0, 1e-3); + CHECK_NEAR(filterCutoffHzFromNorm(1.0f), 20000.0, 1e-2); + + // Exactly three decades, so the decade midpoints land on round numbers. + CHECK_NEAR(filterCutoffHzFromNorm(1.0f / 3.0f), 200.0, 1e-3); + CHECK_NEAR(filterCutoffHzFromNorm(2.0f / 3.0f), 2000.0, 1e-2); + + // Half-decade steps confirm the sweep is log, not linear. + CHECK_NEAR(filterCutoffHzFromNorm(1.0f / 6.0f), 20.0 * std::sqrt(10.0), 1e-3); + CHECK_NEAR(filterCutoffHzFromNorm(0.5f), 20.0 * std::sqrt(1000.0), 1e-2); + + // A linear sweep would put the midpoint at 10010 Hz; the log sweep is nowhere near it. + CHECK(filterCutoffHzFromNorm(0.5f) < 1000.0f); + + CHECK_NEAR(filterCutoffHzFromNorm(-1.0f), 20.0, 1e-3); + CHECK_NEAR(filterCutoffHzFromNorm(2.0f), 20000.0, 1e-2); +} + +static void testCutoffNormRoundTrips() { + for (int i = 0; i <= 20; ++i) { + const float n = static_cast(i) / 20.0f; + CHECK_NEAR(filterNormFromCutoffHz(filterCutoffHzFromNorm(n)), n, 1e-6); + } + CHECK_NEAR(filterNormFromCutoffHz(200.0f), 1.0 / 3.0, 1e-6); + CHECK_NEAR(filterNormFromCutoffHz(2000.0f), 2.0 / 3.0, 1e-6); + CHECK(filterNormFromCutoffHz(1.0f) == 0.0f); + CHECK(filterNormFromCutoffHz(0.0f) == 0.0f); + CHECK(filterNormFromCutoffHz(48000.0f) == 1.0f); +} + +// --------------------------------------------------------------------------- +// Q mapping +// --------------------------------------------------------------------------- + +static void testQSpansPointOneToTenWithRootTwoAtCenter() { + CHECK_NEAR(filterQFromNorm(0.0f), 0.1, 1e-6); + CHECK_NEAR(filterQFromNorm(0.5f), std::sqrt(2.0), 1e-5); + CHECK_NEAR(filterQFromNorm(1.0f), 10.0, 1e-4); + + CHECK_NEAR(filterQFromNorm(-1.0f), 0.1, 1e-6); + CHECK_NEAR(filterQFromNorm(2.0f), 10.0, 1e-4); + + // Strictly monotonic across the whole travel — no fold-back from the quadratic term. + float prev = -1.0f; + for (int i = 0; i <= 1000; ++i) { + const float q = filterQFromNorm(static_cast(i) / 1000.0f); + CHECK(q > prev); + prev = q; + } +} + +static void testQNormRoundTrips() { + for (int i = 0; i <= 20; ++i) { + const float n = static_cast(i) / 20.0f; + CHECK_NEAR(filterNormFromQ(filterQFromNorm(n)), n, 1e-5); + } + CHECK_NEAR(filterNormFromQ(static_cast(std::sqrt(2.0))), 0.5, 1e-5); + CHECK(filterNormFromQ(0.0f) == 0.0f); + CHECK(filterNormFromQ(1000.0f) == 1.0f); +} + +// --------------------------------------------------------------------------- +// Coefficients — pinned literals plus an independent textbook derivation +// --------------------------------------------------------------------------- + +// Textbook RBJ Audio EQ Cookbook, computed straight from cos(w0)/sin(w0). Shares no code with +// filter_coeffs, which reaches the same numbers through the tan half-angle substitution. +static void rbjReference(bool highPass, double fc, double q, double sr, double out[5]) { + const double w0 = 2.0 * kPi * fc / sr; + const double c = std::cos(w0); + const double s = std::sin(w0); + const double alpha = s / (2.0 * q); + const double a0 = 1.0 + alpha; + const double n = highPass ? (1.0 + c) : (1.0 - c); + out[0] = n / 2.0 / a0; // b0 + out[1] = (highPass ? -n : n) / a0; // b1 + out[2] = n / 2.0 / a0; // b2 + out[3] = -2.0 * c / a0; // a1 + out[4] = (1.0 - alpha) / a0; // a2 +} + +static void testCoefficientsMatchPinnedRbjValues() { + const double sr = 48000.0, fc = 1000.0, q = std::sqrt(2.0); + + const BiquadCoeffs lp = biquadCoeffs(FilterMode::LowPass, static_cast(fc), + static_cast(q), sr); + const BiquadCoeffs hp = biquadCoeffs(FilterMode::HighPass, static_cast(fc), + static_cast(q), sr); + + // Pinned literals: change the math and these fail. + CHECK_NEAR(lp.b0, 0.0040888771, 2e-6); + CHECK_NEAR(lp.b1, 0.0081777542, 2e-6); + CHECK_NEAR(lp.b2, 0.0040888771, 2e-6); + CHECK_NEAR(lp.a1, -1.8954199076, 2e-6); + CHECK_NEAR(lp.a2, 0.9117754318, 2e-6); + + CHECK_NEAR(hp.b0, 0.9517988338, 2e-6); + CHECK_NEAR(hp.b1, -1.9035976676, 2e-6); + CHECK_NEAR(hp.b2, 0.9517988338, 2e-6); + CHECK_NEAR(hp.a1, -1.8954199076, 2e-6); + CHECK_NEAR(hp.a2, 0.9117754318, 2e-6); + + // Independent derivation — proves the pinned literals are RBJ and not just "what we emit". + double ref[5]; + rbjReference(false, fc, q, sr, ref); + CHECK_NEAR(lp.b0, ref[0], 1e-6); + CHECK_NEAR(lp.b1, ref[1], 1e-6); + CHECK_NEAR(lp.b2, ref[2], 1e-6); + CHECK_NEAR(lp.a1, ref[3], 1e-6); + CHECK_NEAR(lp.a2, ref[4], 1e-6); + + rbjReference(true, fc, q, sr, ref); + CHECK_NEAR(hp.b0, ref[0], 1e-6); + CHECK_NEAR(hp.b1, ref[1], 1e-6); + CHECK_NEAR(hp.b2, ref[2], 1e-6); + CHECK_NEAR(hp.a1, ref[3], 1e-6); + CHECK_NEAR(hp.a2, ref[4], 1e-6); +} + +static void testCoefficientsTrackSampleRateAndClampBelowNyquist() { + // Same fc at a different rate must give the RBJ answer for THAT rate, not a cached one. + double ref[5]; + rbjReference(false, 1000.0, 2.0, 44100.0, ref); + const BiquadCoeffs at441 = biquadCoeffs(FilterMode::LowPass, 1000.0f, 2.0f, 44100.0); + CHECK_NEAR(at441.a1, ref[3], 1e-6); + CHECK_NEAR(at441.a2, ref[4], 1e-6); + + // Requesting above 0.48*sr clamps rather than diverging through tan(). + const BiquadCoeffs clamped = biquadCoeffs(FilterMode::LowPass, 20000.0f, 1.0f, 32000.0); + rbjReference(false, 0.48 * 32000.0, 1.0, 32000.0, ref); + CHECK_NEAR(clamped.b0, ref[0], 1e-6); + CHECK(std::isfinite(clamped.a1) && std::isfinite(clamped.a2)); + + // A non-positive rate passes through instead of inventing 44.1k. + const BiquadCoeffs bypass = biquadCoeffs(FilterMode::LowPass, 1000.0f, 1.0f, 0.0); + CHECK(bypass.b0 == 1.0f && bypass.b1 == 0.0f && bypass.b2 == 0.0f); + CHECK(bypass.a1 == 0.0f && bypass.a2 == 0.0f); +} + +// DC gain of a lowpass and Nyquist gain of a highpass are both exactly unity — an independent +// structural check on the coefficient set that a sign slip would break. +static void testPassbandGainIsUnity() { + for (double q : {0.1, std::sqrt(2.0), 10.0}) { + const BiquadCoeffs lp = + biquadCoeffs(FilterMode::LowPass, 1000.0f, static_cast(q), 48000.0); + CHECK_NEAR((lp.b0 + lp.b1 + lp.b2) / (1.0 + lp.a1 + lp.a2), 1.0, 1e-4); + + const BiquadCoeffs hp = + biquadCoeffs(FilterMode::HighPass, 1000.0f, static_cast(q), 48000.0); + CHECK_NEAR((hp.b0 - hp.b1 + hp.b2) / (1.0 - hp.a1 + hp.a2), 1.0, 1e-4); + } +} + +// --------------------------------------------------------------------------- +// Resonance +// --------------------------------------------------------------------------- + +// |H(e^jw)| for y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2. +static double magnitudeAt(const BiquadCoeffs& c, double freqHz, double sr) { + const double w = 2.0 * kPi * freqHz / sr; + const double nRe = c.b0 + c.b1 * std::cos(w) + c.b2 * std::cos(2 * w); + const double nIm = -(c.b1 * std::sin(w) + c.b2 * std::sin(2 * w)); + const double dRe = 1.0 + c.a1 * std::cos(w) + c.a2 * std::cos(2 * w); + const double dIm = -(c.a1 * std::sin(w) + c.a2 * std::sin(2 * w)); + return std::sqrt(nRe * nRe + nIm * nIm) / std::sqrt(dRe * dRe + dIm * dIm); +} + +static void testHighQPeaksAtCutoffInBothModes() { + const double sr = 48000.0, fc = 1000.0; + const float qHigh = filterQFromNorm(1.0f); // 10 + const float qLow = filterQFromNorm(0.0f); // 0.1 + + for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) { + const BiquadCoeffs hi = biquadCoeffs(mode, static_cast(fc), qHigh, sr); + + // Scan a log grid and locate the maximum. + double peakMag = 0.0, peakFreq = 0.0; + for (int i = 0; i <= 600; ++i) { + const double f = 20.0 * std::pow(1000.0, static_cast(i) / 600.0); + const double m = magnitudeAt(hi, f, sr); + if (m > peakMag) { peakMag = m; peakFreq = f; } + } + // The peak is at the cutoff, not at a band edge — within a quarter octave. + CHECK(peakFreq > fc / 1.19 && peakFreq < fc * 1.19); + // An RBJ 2-pole peaks at Q; assert most of that emphasis is really there. + CHECK(peakMag > 8.0); + + // The emphasis is relative to the passband, not just a loud filter. + const double passband = magnitudeAt(hi, mode == FilterMode::LowPass ? 20.0 : 20000.0, sr); + CHECK_NEAR(passband, 1.0, 0.05); + CHECK(peakMag / passband > 8.0); + + // At the bottom of the Q control there is no peak at all: the response is monotone + // over the band, so high Q is genuinely doing the work. + const BiquadCoeffs lo = biquadCoeffs(mode, static_cast(fc), qLow, sr); + double prev = magnitudeAt(lo, 20.0, sr); + bool monotone = true; + for (int i = 1; i <= 600; ++i) { + const double f = 20.0 * std::pow(1000.0, static_cast(i) / 600.0); + const double m = magnitudeAt(lo, f, sr); + if (mode == FilterMode::LowPass ? (m > prev + 1e-9) : (m < prev - 1e-9)) { + monotone = false; + } + prev = m; + } + CHECK(monotone); + } +} + +// Drive real sines through VoiceFilter and measure steady-state RMS. Unlike the analytic +// check above this also exercises the high-pass input-feedback path, which is outside the +// coefficient transfer function. +static double measuredRms(FilterMode mode, float cutoffNorm, float resNorm, double freqHz, + double sr) { + VoiceFilter f; + f.prepare({mode, cutoffNorm, resNorm}, sr); + f.reset(); + + const int settle = 24000, measure = 24000; + double sumSq = 0.0; + for (int i = 0; i < settle + measure; ++i) { + const float x = static_cast(std::sin(2.0 * kPi * freqHz * i / sr)); + const float y = f.process(0, x); + if (i >= settle) sumSq += static_cast(y) * y; + } + return std::sqrt(sumSq / measure); +} + +static void testMeasuredResponsePeaksAtCutoffInBothModes() { + const double sr = 48000.0; + const float cutoffNorm = filterNormFromCutoffHz(1000.0f); + + for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) { + double peakRms = 0.0, peakFreq = 0.0; + for (int i = 0; i <= 40; ++i) { + const double f = 100.0 * std::pow(100.0, static_cast(i) / 40.0); + const double r = measuredRms(mode, cutoffNorm, 1.0f, f, sr); + if (r > peakRms) { peakRms = r; peakFreq = f; } + } + CHECK(peakFreq > 1000.0 / 1.3 && peakFreq < 1000.0 * 1.3); + + const double passband = + measuredRms(mode, cutoffNorm, 1.0f, mode == FilterMode::LowPass ? 100.0 : 10000.0, sr); + CHECK(peakRms / passband > 3.0); + + // Same measurement at the bottom of the resonance control shows no such emphasis. + const double flatAtCutoff = measuredRms(mode, cutoffNorm, 0.0f, 1000.0, sr); + const double flatPassband = + measuredRms(mode, cutoffNorm, 0.0f, mode == FilterMode::LowPass ? 100.0 : 10000.0, sr); + CHECK(flatAtCutoff / flatPassband < 1.0); + } +} + +// --------------------------------------------------------------------------- +// Stability +// --------------------------------------------------------------------------- + +static void testFullRangeCutoffSweepAtAudioRateStaysBounded() { + // Deterministic pseudo-noise; a fixed sine would miss the resonant frequency on most steps. + unsigned rng = 0x13579bdfu; + auto noise = [&rng]() { + rng = rng * 1664525u + 1013904223u; + return static_cast(static_cast(rng >> 9) - (1 << 22)) / static_cast(1 << 22); + }; + + for (double sr : {44100.0, 48000.0, 96000.0}) { + for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) { + for (float res : {0.0f, 0.5f, 1.0f}) { + for (int direction = 0; direction < 2; ++direction) { + VoiceFilter f; + f.reset(); + const int n = 48000; + for (int i = 0; i < n; ++i) { + const float t = static_cast(i) / static_cast(n - 1); + // Per-sample coefficient update across the whole cutoff travel. + f.prepare({mode, direction == 0 ? t : 1.0f - t, res}, sr); + const float y = f.process(0, noise()); + CHECK(std::isfinite(y)); + CHECK(std::fabs(y) < 100.0f); + if (!std::isfinite(y)) return; // stop before the log floods + } + } + } + } + } +} + +static void testStateFlushesToZeroWithoutStallingInDenormals() { + const double sr = 48000.0; + for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) { + VoiceFilter f; + f.prepare({mode, filterNormFromCutoffHz(1000.0f), 1.0f}, sr); + f.reset(); + + // Excite, then hard-cut to silence the way a released voice does. + for (int i = 0; i < 480; ++i) { + f.process(0, 0.5f * static_cast(std::sin(2.0 * kPi * 1000.0 * i / sr))); + } + + int subnormalSamples = 0; + int silentAt = -1; + for (int i = 0; i < 20000; ++i) { + f.process(0, 0.0f); + const VoiceFilter::State& s = f.state(0); + const float vals[5] = {s.x1, s.x2, s.y1, s.y2, s.fb}; + for (float v : vals) { + if (v != 0.0f && std::fabs(v) < FLT_MIN) { ++subnormalSamples; break; } + } + if (silentAt < 0 && f.isSilent()) silentAt = i; + } + // Without the flush the state grinds down through the subnormal range for thousands + // of samples; a stray sample or two at a zero crossing is not a stall. + CHECK(subnormalSamples <= 2); + CHECK(silentAt >= 0); + CHECK(silentAt < 20000); + // And it stays silent — a flush that perturbs the feedback loop would re-excite it. + for (int i = 0; i < 1000; ++i) CHECK(f.process(0, 0.0f) == 0.0f); + CHECK(f.isSilent()); + } +} + +// --------------------------------------------------------------------------- +// Impulse / step sanity and saturation +// --------------------------------------------------------------------------- + +static void testImpulseResponseMatchesDifferenceEquation() { + const double sr = 48000.0; + VoiceFilter f; + f.prepare({FilterMode::LowPass, filterNormFromCutoffHz(1000.0f), 0.5f}, sr); + f.reset(); + const BiquadCoeffs c = f.coeffs(); + + // First three impulse-response taps follow directly from the coefficients. + const float h0 = f.process(0, 1.0f); + const float h1 = f.process(0, 0.0f); + const float h2 = f.process(0, 0.0f); + CHECK_NEAR(h0, c.b0, 1e-6); + CHECK_NEAR(h1, c.b1 - c.a1 * c.b0, 1e-6); + CHECK_NEAR(h2, c.b2 - c.a1 * h1 - c.a2 * h0, 1e-6); +} + +static void testLowpassStepSettlesToUnity() { + const double sr = 48000.0; + VoiceFilter f; + f.prepare({FilterMode::LowPass, filterNormFromCutoffHz(1000.0f), 0.0f}, sr); + f.reset(); + float y = 0.0f; + for (int i = 0; i < 48000; ++i) y = f.process(0, 1.0f); + CHECK_NEAR(y, 1.0, 1e-3); // DC passes a lowpass at unity + + VoiceFilter hp; + hp.prepare({FilterMode::HighPass, filterNormFromCutoffHz(1000.0f), 0.0f}, sr); + hp.reset(); + for (int i = 0; i < 48000; ++i) y = hp.process(0, 1.0f); + CHECK_NEAR(y, 0.0, 1e-3); // and is fully rejected by a highpass +} + +static void testResetClearsHistoryButPrepareKeepsIt() { + VoiceFilter f; + f.prepare({FilterMode::LowPass, 0.5f, 0.5f}, 48000.0); + f.process(0, 1.0f); + CHECK(!f.isSilent()); + + // A live parameter move must not zero the history — that is what would click. + f.prepare({FilterMode::LowPass, 0.6f, 0.5f}, 48000.0); + CHECK(!f.isSilent()); + + f.reset(); + CHECK(f.isSilent()); +} + +static void testChannelStateIsIndependent() { + VoiceFilter f; + f.prepare({FilterMode::LowPass, 0.5f, 0.5f}, 48000.0); + f.reset(); + f.process(0, 1.0f); + CHECK(f.state(0).x1 == 1.0f); + CHECK(f.state(1).x1 == 0.0f); + + float frame[2] = {1.0f, -1.0f}; + f.processFrame(frame, 2); + CHECK(f.state(1).x1 == -1.0f); + CHECK(frame[0] != frame[1]); +} + +static void testFeedbackSaturationIsContinuousAndBounded() { + CHECK_NEAR(feedbackSaturate(0.0f), 0.0, 1e-9); + // Odd symmetry. + CHECK_NEAR(feedbackSaturate(1.5f), -feedbackSaturate(-1.5f), 1e-6); + // Continuous across the threshold at +/-2. + CHECK_NEAR(feedbackSaturate(2.0f - 1e-4f), feedbackSaturate(2.0f + 1e-4f), 1e-4); + // Compresses hard: a 100x input does not give a 100x output. + CHECK(std::fabs(feedbackSaturate(100.0f)) < 12.0f); + CHECK(feedbackSaturate(100.0f) > feedbackSaturate(50.0f)); +} + +int main() { + testCutoffMapsThreeDecadesLogarithmically(); + testCutoffNormRoundTrips(); + testQSpansPointOneToTenWithRootTwoAtCenter(); + testQNormRoundTrips(); + testCoefficientsMatchPinnedRbjValues(); + testCoefficientsTrackSampleRateAndClampBelowNyquist(); + testPassbandGainIsUnity(); + testHighQPeaksAtCutoffInBothModes(); + testMeasuredResponsePeaksAtCutoffInBothModes(); + testFullRangeCutoffSweepAtAudioRateStaysBounded(); + testStateFlushesToZeroWithoutStallingInDenormals(); + testImpulseResponseMatchesDifferenceEquation(); + testLowpassStepSettlesToUnity(); + testResetClearsHistoryButPrepareKeepsIt(); + testChannelStateIsIndependent(); + testFeedbackSaturationIsContinuousAndBounded(); + + if (g_fail == 0) std::printf("filter_tests: all passed\n"); + else std::printf("filter_tests: %d FAILED\n", g_fail); + return g_fail == 0 ? 0 : 1; +}