From d2364eb5ac0c18b913570be40c6632294a51c2ea Mon Sep 17 00:00:00 2001 From: daniel-c-harvey Date: Thu, 30 Jul 2026 10:22:49 -0400 Subject: [PATCH] Add an Oberheim-SEM morph law to the SVF filter: HP->notch->LP alongside HP->BP->LP, selected at prepare() time, free on the per-sample path --- CMakeLists.txt | 30 +- src/core/instrument/engine/filter/CLAUDE.md | 68 +- .../instrument/engine/filter/filter_morph.cpp | 12 +- .../instrument/engine/filter/filter_morph.h | 57 +- .../instrument/engine/filter/voice_filter.cpp | 2 +- .../instrument/engine/filter/voice_filter.h | 5 +- tests/test_filter.cpp | 704 +++++++----------- tests/test_filter_morph.cpp | 208 ++++++ tests/test_filter_params.cpp | 121 +++ tests/test_filter_state.cpp | 247 ++++++ 10 files changed, 971 insertions(+), 483 deletions(-) create mode 100644 tests/test_filter_morph.cpp create mode 100644 tests/test_filter_params.cpp create mode 100644 tests/test_filter_state.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e770077..65434a6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -995,7 +995,8 @@ 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 TPT/SVF with a continuous HP->BP->LP morph and an in-loop drive stage. +# filter — the per-voice TPT/SVF with a continuous morph (HP->BP->LP or HP->notch->LP, selected +# at prepare() time) and an in-loop drive stage. # The Cortex-M4 source's virtual FilterBase/Filter/Biquad hierarchy dispatched per channel per # sample, which the per-voice per-sample path forbids, so none of it came across. Control # mapping, SVF coefficients, morph weights, and the filter type each get their own file; @@ -1108,10 +1109,29 @@ 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 SVF coefficients against an independent -# derivation, asserts the cutoff/Q control mappings at their anchors, holds the morph endpoints -# to the analytic 2-pole targets, and measures rate/level invariance and drive stability by -# driving real sines. NEITHER SDK. +# filter: four targets along the module's own seams, so each asserts one domain. +# filter_params_tests — the rate-free control mappings (cutoff/Q/drive) and their inverses. +# filter_morph_tests — the pure morph-weight algebra under both morph laws; no DSP is run. +# filter_state_tests — numerical stability, the denormal flush, and the state lifecycle. +# filter_tests — the frequency response: pins the SVF coefficients against an independent +# derivation, holds the morph endpoints to the analytic 2-pole targets, and measures the +# HP-BP-LP corner flatness, the HP-notch-LP null, and rate/level invariance and drive +# stability by driving real sines. The seams above were chosen so this file alone owns the +# analytic reference and the steady-state gain measurement — a forked copy of a measurement +# reference is a worse defect than a long file, which is why it sits over the ~600-line bar. +# NEITHER SDK. +add_executable(filter_params_tests tests/test_filter_params.cpp) +target_link_libraries(filter_params_tests PRIVATE filter) +add_test(NAME filter_params_tests COMMAND filter_params_tests) + +add_executable(filter_morph_tests tests/test_filter_morph.cpp) +target_link_libraries(filter_morph_tests PRIVATE filter) +add_test(NAME filter_morph_tests COMMAND filter_morph_tests) + +add_executable(filter_state_tests tests/test_filter_state.cpp) +target_link_libraries(filter_state_tests PRIVATE filter) +add_test(NAME filter_state_tests COMMAND filter_state_tests) + add_executable(filter_tests tests/test_filter.cpp) target_link_libraries(filter_tests PRIVATE filter) add_test(NAME filter_tests COMMAND filter_tests) diff --git a/src/core/instrument/engine/filter/CLAUDE.md b/src/core/instrument/engine/filter/CLAUDE.md index 577ff38..3c016e3 100644 --- a/src/core/instrument/engine/filter/CLAUDE.md +++ b/src/core/instrument/engine/filter/CLAUDE.md @@ -3,7 +3,8 @@ ## Scope The pure per-voice filter a sounding voice runs: a Zavalishin TPT/SVF with a continuous -HP→BP→LP morph and a drive stage. No REAPER, no VST3, no allocation, no I/O. Everything +morph under one of two laws — HP→BP→LP or HP→notch→LP — and a drive stage. No REAPER, no +VST3, no allocation, no I/O. Everything here lives in `reasampler::instrument::engine::filter`, nested per the directory-mirrors-namespace convention — this keeps `FilterSettings` and friends out of `reasampler::instrument::engine` proper, where `zone_params.h` lives, since this module has @@ -14,8 +15,9 @@ responsibility each: drive depth, plus the exact inverses for cutoff and Q. - `filter_coeffs` — the DSP domain: `SvfCoeffs` and the TPT coefficient solve from (cutoff Hz, Q, sample rate). -- `filter_morph` — the morph domain: normalized position → per-tap weights, and the fold of - those weights into the three multipliers the kernel applies. +- `filter_morph` — the morph domain: `MorphLaw`, normalized position → per-tap weights under + the selected law, and the fold of those weights into the three multipliers the kernel + applies. - `filter_saturate` — `softLimit`, the drive stage's shaper. Header-only inline; it sits inside the per-sample recursion. - `voice_filter` — `FilterSettings` and `VoiceFilter`, the concrete per-voice type. @@ -63,21 +65,50 @@ now an explicit user-controlled stage instead of an emergent side effect. ### The morph is a blend of taps, never a coefficient switch An SVF produces high, band, and low from the same state, which is the reason this topology -was chosen. `FilterMode` as a discrete enum is retired. HP at 0.0, BP at 0.5, LP at 1.0, -continuous throughout, and the three endpoints are exact. +was chosen. `FilterMode` as a discrete enum is retired. HP at 0.0, LP at 1.0, continuous +throughout, and both endpoints are exact under either law — only the centre differs. -The crossfade is **equal-power between adjacent taps**, and both halves of that are forced -by the topology rather than picked by ear: +The crossfade is **equal-power** in both laws, and that is forced by the topology rather +than picked by ear. At the corner the taps are `HP = jQ`, `BP = Q`, `LP = -jQ` — adjacent +taps in exact quadrature and HP/LP in exact antiphase, relationships the bilinear transform +preserves exactly at the prewarped corner. A `cos`/`sin` pair therefore holds the crossfaded +power at unity across the whole sweep; a linear crossfade of a quadrature pair would sag to +`1/sqrt(2)` mid-leg, a 3 dB hole that reads as a defect rather than as character. -- At the corner the taps are `HP = jQ`, `BP = Q`, `LP = -jQ` — adjacent taps in exact - quadrature, which the bilinear transform preserves exactly at the prewarped corner. A - `cos`/`sin` pair therefore holds the corner magnitude at exactly `Q*sqrt(cos² + sin²) = Q` - at every morph position. A linear crossfade of a quadrature pair would sag to `Q/sqrt(2)` - mid-leg — a 3 dB hole that reads as a defect, not as character. -- **Adjacent only.** HP and LP are exactly antiphase at the corner, so any law giving both - simultaneous weight cancels there and cuts a notch. That notch is the Oberheim SEM's - centre tap. This control's centre is a band-pass, per the explicit HP/BP/LP enumeration — - do not "simplify" the two legs into one three-way weighting, which silently builds the SEM. +### The two morph laws, and why only one of them has a flat corner + +`MorphLaw` is a two-value selector on `FilterSettings`, **defaulting to `HighBandLow`** — +that is the reviewed-and-measured law, and it is enumerator 0 so a zero-initialized or absent +persisted field lands on it rather than on the SEM leg. + +- **`HighBandLow` (HP→BP→LP, the default).** Two equal-power legs crossfading **adjacent taps + only**, BP at the centre. Because adjacent taps are in quadrature, the corner magnitude is + algebraically `Q*sqrt(cos² + sin²) = Q` at every position — measured flat to 4e-6 across 65 + positions. **That flatness guarantee is specific to this law.** Do not weaken the assertion + that pins it in order to accommodate the other law. +- **`HighNotchLow` (HP→notch→LP, the Oberheim SEM).** One equal-power crossfade weighting HP + and LP **together** across the whole sweep, `bp == 0` throughout. The notch is not tuned in: + HP and LP sit at exactly +90° and −90° at the corner, so equal weights cancel there by + construction. Here the corner magnitude deliberately goes to **zero** at the centre — + measured worst case −88 dB across every rate/cutoff/Q, typically −110 to −145 dB. The fold + makes that structural rather than a runtime near-miss: `m2 = lp - hp` is **exactly** `0.0f` + at the centre, because `cos` and `sin` of π/4 differ by about an ulp of *double*, nine + orders below float's spacing there, so they narrow to one float. + +SEM's zero is at the **notch frequency**, not a broadband level sag — off the corner the pair +is still equal-power, so neither law's legs dip. Measuring that requires dividing by each +tap's own analytic response first: at `Q = 0.1` a 2-pole approaches its passband so slowly +that the pure low tap still reads 0.896 at 50 Hz, and a raw reading would report a 20% "sag" +that is the Q, not the morph. + +**The toggle is free on the hot path, and must stay that way.** `morphWeights` runs at +`prepare()` cadence; the law is consumed there and nowhere else. The kernel, `svfCoeffs`, and +`morphMix`'s fold are identical between the laws — all a law selects is three floats the +kernel was already multiplying by. Verified at the machine-code level, not by inspection: the +same TU compiled `/O2` against the pre-toggle and post-toggle headers emits byte-identical +assembly for `process()` and `processFrame()`. `VoiceFilter` gained no member and `process()` +gained no branch. A design that puts the law selector inside the per-sample path is wrong — +rework it rather than paying for it. ### Drive is a contraction inside the loop, which is what makes it unconditionally stable @@ -183,6 +214,11 @@ topology. - **The morph endpoints are asserted on the folded mix, exactly.** `morphWeights` snaps the leg endpoints instead of trusting `cos`/`sin` to land on 0 and 1, which they miss by ~1e-17 — enough to leave a -324 dB neighbour tap in what is specified as a pure response. +- **A NaN morph position falls back per law, not to one shared value.** Every comparison + against NaN is false, so it clamps to neither endpoint: `HighBandLow` lands on pure + band-pass, `HighNotchLow` on pure high-pass, since it has no band tap to land on. +- **Measuring a null needs a ring-time-adequate settle window.** At `Q = 10` the leftover + transient alone reads as −52 dB after 0.15 s and would be mistaken for the noise floor. - **No call site yet.** Wiring the filter into the voice path is a separate track; nothing in `sampler_core` references this module today. - **Decay to the denormal floor is a fixed wall-clock time, not a sample count.** A test diff --git a/src/core/instrument/engine/filter/filter_morph.cpp b/src/core/instrument/engine/filter/filter_morph.cpp index e06ad69..6ed05c8 100644 --- a/src/core/instrument/engine/filter/filter_morph.cpp +++ b/src/core/instrument/engine/filter/filter_morph.cpp @@ -23,10 +23,20 @@ Pair equalPower(double t) { } // namespace -MorphWeights morphWeights(float norm) { +MorphWeights morphWeights(float norm, MorphLaw law) { const double n = norm < 0.0 ? 0.0 : (norm > 1.0 ? 1.0 : static_cast(norm)); MorphWeights w; + if (law == MorphLaw::HighNotchLow) { + // ONE crossfade across the whole sweep rather than two legs, so HP and LP carry weight + // together everywhere between the endpoints and are equal at the centre. + const Pair p = equalPower(n); + w.hp = static_cast(p.a); + w.bp = 0.0f; + w.lp = static_cast(p.b); + return w; + } + if (n <= 0.5) { const Pair p = equalPower(2.0 * n); // HP -> BP w.hp = static_cast(p.a); diff --git a/src/core/instrument/engine/filter/filter_morph.h b/src/core/instrument/engine/filter/filter_morph.h index 15b484f..269ae2b 100644 --- a/src/core/instrument/engine/filter/filter_morph.h +++ b/src/core/instrument/engine/filter/filter_morph.h @@ -1,34 +1,55 @@ -// filter_morph.h — the continuous HP -> BP -> LP morph: normalized position to tap weights, -// and the fold of those weights into the three multipliers the kernel actually applies. An SVF -// produces all three taps from one state, so the morph is a blend, never a coefficient switch. +// filter_morph.h — the continuous morph: normalized position to tap weights under one of two +// laws, and the fold of those weights into the three multipliers the kernel actually applies. +// An SVF produces all three taps from one state, so the morph is a blend, never a coefficient +// switch. Weights are computed at prepare() cadence; the law never reaches the per-sample path. #pragma once namespace reasampler::instrument::engine::filter { -// Weight on each SVF tap. Exactly one of hp/lp is nonzero at a time — the morph crossfades -// between ADJACENT taps only, never HP against LP. +// Which shape the sweep traces between its two fixed endpoints. This selects CHARACTER, not +// topology — same SVF, same coefficients, same kernel under either law; only the centre differs. +// +// HighBandLow is enumerator 0 deliberately: a zero-initialized or absent persisted field then +// lands on the default rather than on the SEM leg. +enum class MorphLaw { + // HP -> BP -> LP. Crossfades ADJACENT taps only, so the corner magnitude is flat at Q the + // whole way across. The default. + HighBandLow, + // HP -> notch -> LP, the Oberheim SEM. One crossfade weighting HP and LP together, bp == 0 + // throughout; the notch falls out of the antiphase cancellation rather than being tuned in. + HighNotchLow, +}; + +// Weight on each SVF tap. Under HighBandLow exactly one of hp/lp is nonzero at a time — that law +// crossfades adjacent taps only, never HP against LP. Under HighNotchLow bp is always zero and +// hp/lp carry weight together, which is precisely what cuts the notch. struct MorphWeights { float hp = 0.0f; float bp = 0.0f; float lp = 1.0f; }; -// HP at 0.0, BP at 0.5, LP at 1.0. Out-of-range norm clamps to the endpoints; NaN clamps to -// neither and lands on pure band-pass instead (every comparison against it is false). +// HP at 0.0, LP at 1.0 under BOTH laws; the centre is a band-pass under HighBandLow and a notch +// under HighNotchLow. Out-of-range norm clamps to the endpoints; NaN clamps to neither (every +// comparison against it is false) and lands on the law's degenerate — pure band-pass under +// HighBandLow, pure high-pass under HighNotchLow, which has no band tap to land on. // -// Equal-power (cos/sin) rather than linear, and that choice is forced by the topology rather -// than picked by ear. At the corner frequency the three taps are HP = jQ, BP = Q, LP = -jQ, so -// adjacent taps are in exact QUADRATURE there (and the bilinear transform preserves that exactly -// at the prewarped corner). Under a cos/sin pair the corner magnitude is therefore -// Q*sqrt(cos^2 + sin^2) = Q at every morph position — algebraically flat across the whole sweep. -// A linear crossfade of the same quadrature pair would sag to Q/sqrt(2), a 3 dB hole mid-leg. +// Equal-power (cos/sin) in both laws rather than linear, and that choice is forced by the +// topology rather than picked by ear. At the corner frequency the three taps are HP = jQ, +// BP = Q, LP = -jQ, so ADJACENT taps are in exact QUADRATURE there (and the bilinear transform +// preserves that exactly at the prewarped corner). Under HighBandLow's cos/sin pair the corner +// magnitude is therefore Q*sqrt(cos^2 + sin^2) = Q at every morph position — algebraically flat +// across the whole sweep. A linear crossfade of the same quadrature pair would sag to Q/sqrt(2), +// a 3 dB hole mid-leg. // -// Crossfading adjacent taps only is the other half of it: HP and LP are exactly ANTIPHASE at the -// corner, so any law giving both simultaneous weight cancels there and cuts a notch. That notch -// is the Oberheim SEM's center tap; this control's center is a band-pass, per the explicit -// HP/BP/LP enumeration. -MorphWeights morphWeights(float norm); +// HP and LP are exactly ANTIPHASE at the corner (+90 and -90 degrees), so a law giving both +// simultaneous weight cancels there. HighBandLow avoids that by staying adjacent; HighNotchLow +// uses it — one equal-power crossfade of HP against LP over the whole sweep puts equal weights +// at the centre and the null is exact by construction, not tuned. That is why the corner-flat-at-Q +// guarantee is specific to HighBandLow: on the SEM leg the corner magnitude deliberately goes to +// zero at the centre. Equal power still holds off the notch frequency, so neither law's legs sag. +MorphWeights morphWeights(float norm, MorphLaw law); // The kernel applies out = m0*v0 + m1*v1 + m2*v2, where v0 is the input and v1/v2 are the SVF's // band and low outputs. Folding hp = v0 - k*v1 - v2 into the weights here keeps the per-sample diff --git a/src/core/instrument/engine/filter/voice_filter.cpp b/src/core/instrument/engine/filter/voice_filter.cpp index 2b0f369..b22278e 100644 --- a/src/core/instrument/engine/filter/voice_filter.cpp +++ b/src/core/instrument/engine/filter/voice_filter.cpp @@ -6,7 +6,7 @@ void VoiceFilter::prepare(const FilterSettings& settings, double sampleRate) { coeffs_ = svfCoeffs(filterCutoffHzFromNorm(settings.cutoffNorm), filterQFromNorm(settings.resonanceNorm), sampleRate); if (sampleRate > 0.0) { - mix_ = morphMix(morphWeights(settings.morphNorm), coeffs_.k); + mix_ = morphMix(morphWeights(settings.morphNorm, settings.morphLaw), coeffs_.k); } else { // Bypass: a1=1, a2=a3=0 makes both state updates the exact identity, and bypassMix() // reads only the input, never the state -- so clearing here is audibly free (the state diff --git a/src/core/instrument/engine/filter/voice_filter.h b/src/core/instrument/engine/filter/voice_filter.h index 079df92..f9b83cb 100644 --- a/src/core/instrument/engine/filter/voice_filter.h +++ b/src/core/instrument/engine/filter/voice_filter.h @@ -15,11 +15,14 @@ namespace reasampler::instrument::engine::filter { // Normalized control positions, as the editor moves them and the persisted state carries them. +// morphLaw is the one discrete control here — a two-value selector, not a normalized position — +// because its two values are characters to choose between, not points on a continuum. struct FilterSettings { float cutoffNorm = 1.0f; float resonanceNorm = 0.0f; - float morphNorm = 1.0f; // 0 = high-pass, 0.5 = band-pass, 1 = low-pass + float morphNorm = 1.0f; // 0 = high-pass, 1 = low-pass; the centre is set by morphLaw float driveNorm = 0.0f; + MorphLaw morphLaw = MorphLaw::HighBandLow; }; // Below this the recursion has decayed past -600 dB. Flushing keeps the state out of the diff --git a/tests/test_filter.cpp b/tests/test_filter.cpp index 390bf7b..8a4bbf4 100644 --- a/tests/test_filter.cpp +++ b/tests/test_filter.cpp @@ -1,8 +1,11 @@ -// Standalone tests for the per-voice TPT/SVF 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 derivation that -// shares no code with the implementation, and the responses against the analog 2-pole prototype -// evaluated at the bilinear-warped frequency. +// Standalone tests for the RUNNING per-voice TPT/SVF 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 derivation +// that shares no code with the implementation, and the responses against the analog 2-pole +// prototype evaluated at the bilinear-warped frequency. Sibling targets own the neighbouring +// domains: test_filter_params.cpp the control mappings, test_filter_morph.cpp the pure morph-weight +// algebra, test_filter_state.cpp the numerical/state behaviour. This file owns the analytic +// reference and the steady-state gain measurement, and everything here uses them. #include "../src/core/instrument/engine/filter/filter_coeffs.h" #include "../src/core/instrument/engine/filter/filter_morph.h" @@ -10,11 +13,9 @@ #include "../src/core/instrument/engine/filter/filter_saturate.h" #include "../src/core/instrument/engine/filter/voice_filter.h" -#include #include #include #include -#include using namespace reasampler::instrument::engine::filter; @@ -28,11 +29,18 @@ static int g_fail = 0; static constexpr double kPi = 3.14159265358979323846; -// Morph positions of the three pure taps. +// Morph positions. The endpoints are the same pure taps under both laws; only the centre differs +// — a band-pass under HighBandLow, a notch under HighNotchLow. static constexpr float kHighPass = 0.0f; static constexpr float kBandPass = 0.5f; +static constexpr float kCentre = 0.5f; static constexpr float kLowPass = 1.0f; +static const MorphLaw kBothLaws[] = {MorphLaw::HighBandLow, MorphLaw::HighNotchLow}; +static const char* lawName(MorphLaw law) { + return law == MorphLaw::HighBandLow ? "HP-BP-LP" : "HP-notch-LP"; +} + // The rates the invariance claims are made over. static const double kRates[] = {44100.0, 48000.0, 88200.0, 96000.0, 192000.0}; static constexpr int kRateCount = 5; @@ -75,101 +83,9 @@ static double measuredGain(const FilterSettings& fs, double sr, double freq, dou return std::sqrt(sumSq / measure) / (amp / std::sqrt(2.0)); } -static FilterSettings at(double fcHz, float res, float morph, float drive = 0.0f) { - return {filterNormFromCutoffHz(static_cast(fcHz)), res, morph, drive}; -} - -// --------------------------------------------------------------------------- -// Control mappings (carried over — the cutoff and Q laws are unchanged) -// --------------------------------------------------------------------------- - -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); - - 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(filterNormFromCutoffHz(1.0f) == 0.0f); - CHECK(filterNormFromCutoffHz(48000.0f) == 1.0f); -} - -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); - - // Pins the single quadratic-in-log-Q curve at two interior points, derived independently by - // solving log Q = a + b*n + c*n^2 through the three anchors above rather than read out of - // the implementation. A two-spliced-log-segments curve (log-linear on each half, the design - // this module doc explicitly rejects for its center-detent slope kink) would give 0.376 and - // 3.761 here instead — both comfortably outside this tolerance. - { - const double lo = std::log(static_cast(kFilterQMin)); - const double mid = std::log(static_cast(kFilterQCenter)); - const double hi = std::log(static_cast(kFilterQMax)); - const double c = 2.0 * lo + 2.0 * hi - 4.0 * mid; - const double b = hi - lo - c; - const double a = lo; - auto qLaw = [&](double n) { return std::exp(a + b * n + c * n * n); }; - CHECK_NEAR(filterQFromNorm(0.25f), qLaw(0.25), 1e-5); - CHECK_NEAR(filterQFromNorm(0.75f), qLaw(0.75), 1e-5); - } - - // 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); -} - -static void testDriveDepthIsZeroAtRestAndRisesMonotonically() { - // Exactly zero, not nearly: the limiter is the identity only at depth 0. - CHECK(filterDriveDepthFromNorm(0.0f) == 0.0f); - CHECK(filterDriveDepthFromNorm(-1.0f) == 0.0f); - CHECK_NEAR(filterDriveDepthFromNorm(1.0f), kFilterDriveDepthMax, 1e-6); - CHECK_NEAR(filterDriveDepthFromNorm(2.0f), kFilterDriveDepthMax, 1e-6); - - // Pins the SQUARE law at an interior point, not just the anchors: a linear law would give - // kFilterDriveDepthMax/2 (2.0) here, not kFilterDriveDepthMax/4 (1.0). - CHECK_NEAR(filterDriveDepthFromNorm(0.5f), kFilterDriveDepthMax * 0.25, 1e-6); - - float prev = -1.0f; - for (int i = 0; i <= 100; ++i) { - const float d = filterDriveDepthFromNorm(static_cast(i) / 100.0f); - CHECK(d > prev); - prev = d; - } +static FilterSettings at(double fcHz, float res, float morph, float drive = 0.0f, + MorphLaw law = MorphLaw::HighBandLow) { + return {filterNormFromCutoffHz(static_cast(fcHz)), res, morph, drive, law}; } // --------------------------------------------------------------------------- @@ -253,88 +169,47 @@ static void testNonPositiveRatePassesSignalThroughAtEveryMorph() { } // --------------------------------------------------------------------------- -// Morph +// Morph — measured, under both laws // --------------------------------------------------------------------------- -// The endpoints are pure taps EXACTLY, not to within a rounding of cos/sin. Asserted on the -// folded mix, where "pure" is an exact statement about three floats. -static void testMorphEndpointMixesAreExactlyPureTaps() { - const float k = 1.0f / filterQFromNorm(0.5f); - - const MorphMix hp = morphMix(morphWeights(kHighPass), k); - CHECK(hp.m0 == 1.0f && hp.m1 == -k && hp.m2 == -1.0f); // v0 - k*v1 - v2 - - const MorphMix bp = morphMix(morphWeights(kBandPass), k); - CHECK(bp.m0 == 0.0f && bp.m1 == 1.0f && bp.m2 == 0.0f); // v1 - - const MorphMix lp = morphMix(morphWeights(kLowPass), k); - CHECK(lp.m0 == 0.0f && lp.m1 == 0.0f && lp.m2 == 1.0f); // v2 - - // Out-of-range clamps to the endpoints rather than extrapolating. - CHECK(morphWeights(-1.0f).hp == 1.0f); - CHECK(morphWeights(2.0f).lp == 1.0f); - - // NaN clamps to neither endpoint (every comparison against it is false) and lands on pure - // band-pass instead -- no crash, a sane fallback rather than an extrapolation. - const MorphWeights nanW = morphWeights(std::numeric_limits::quiet_NaN()); - CHECK(nanW.hp == 0.0f && nanW.bp == 1.0f && nanW.lp == 0.0f); -} - -// Pins the cos/sin curve at an interior point, not just the endpoints and the quadrature -// identity (hp^2+bp^2+lp^2=1, which any equal-power reparameterization would also satisfy). -// theta=0.5*pi*t^2 (quadratic in the leg fraction, still equal-power, still exact at both -// ends) would give hp=0.9239/bp=0.3827 here instead of the cos/sin pair's 0.7071/0.7071. -static void testMorphInteriorPointMatchesCosSinNotAnAlternateEqualPowerCurve() { - const MorphWeights w = morphWeights(0.25f); // HP->BP leg, t = 2*0.25 = 0.5 - const double theta = 0.5 * kPi * 0.5; - CHECK_NEAR(w.hp, std::cos(theta), 1e-6); - CHECK_NEAR(w.bp, std::sin(theta), 1e-6); - CHECK(w.lp == 0.0f); -} - -// HP and LP never carry weight at the same time. That is what keeps the centre a band-pass -// instead of the Oberheim SEM's notch: the two are antiphase at the corner and would cancel. -static void testMorphNeverBlendsHighAgainstLowPass() { - for (int i = 0; i <= 200; ++i) { - const MorphWeights w = morphWeights(static_cast(i) / 200.0f); - CHECK(w.hp == 0.0f || w.lp == 0.0f); - CHECK(w.hp >= 0.0f && w.bp >= 0.0f && w.lp >= 0.0f); - // Equal power: the active pair sums in quadrature to unity. - CHECK_NEAR(w.hp * w.hp + w.bp * w.bp + w.lp * w.lp, 1.0, 1e-6); - } -} - +// The endpoints are exact 2-pole HP and LP under BOTH laws; only the centre is law-specific, so +// the centre is asserted here only for the law that has a pure tap there. static void testMorphEndpointsMatchTheAnalyticTwoPoleTargets() { const double sr = 48000.0, fc = 1000.0; - for (float res : {0.0f, 0.5f, 1.0f}) { - const double q = filterQFromNorm(res); - for (float morph : {kHighPass, kBandPass, kLowPass}) { - for (double f : {125.0, 500.0, 1000.0, 2000.0, 8000.0}) { - const double got = measuredGain(at(fc, res, morph), sr, f); - const double want = analyticMag(morph, f, fc, q, sr); - if (!(std::fabs(got / want - 1.0) <= kAgreement)) { - std::printf("FAIL line %d: morph %.1f res %.1f at %.0f Hz: %.6f vs analytic " - "%.6f (%.3f%%)\n", - __LINE__, morph, res, f, got, want, - (got / want - 1.0) * 100.0); - ++g_fail; + for (MorphLaw law : kBothLaws) { + for (float res : {0.0f, 0.5f, 1.0f}) { + const double q = filterQFromNorm(res); + for (float morph : {kHighPass, kBandPass, kLowPass}) { + if (morph == kBandPass && law != MorphLaw::HighBandLow) continue; + for (double f : {125.0, 500.0, 1000.0, 2000.0, 8000.0}) { + const double got = measuredGain(at(fc, res, morph, 0.0f, law), sr, f); + const double want = analyticMag(morph, f, fc, q, sr); + if (!(std::fabs(got / want - 1.0) <= kAgreement)) { + std::printf("FAIL line %d: %s morph %.1f res %.1f at %.0f Hz: %.6f vs " + "analytic %.6f (%.3f%%)\n", + __LINE__, lawName(law), morph, res, f, got, want, + (got / want - 1.0) * 100.0); + ++g_fail; + } } } } } } -// The reason the blend is equal-power rather than linear. At the corner the three taps are -// HP = jQ, BP = Q, LP = -jQ — adjacent taps in exact quadrature — so a cos/sin pair holds the -// corner magnitude at exactly Q the whole way across. A linear crossfade would sag to Q/sqrt(2) -// mid-leg, a 3 dB hole that would read as a defect rather than as character. -static void testCornerMagnitudeIsFlatAcrossTheWholeMorphSweep() { +// LAW-SPECIFIC, and deliberately not generalized: this guarantee belongs to HighBandLow alone. +// At the corner the three taps are HP = jQ, BP = Q, LP = -jQ — ADJACENT taps in exact quadrature +// — so a cos/sin pair holds the corner magnitude at exactly Q the whole way across. A linear +// crossfade would sag to Q/sqrt(2) mid-leg, a 3 dB hole that would read as a defect rather than +// as character. HighNotchLow deliberately violates this (its corner magnitude goes to zero at the +// centre); weakening this assertion to accommodate that law would throw the guarantee away. +static void testCornerMagnitudeIsFlatAtQAcrossTheHighBandLowSweep() { const double sr = 48000.0, fc = 1000.0; for (float res : {0.0f, 0.5f, 1.0f}) { const double q = filterQFromNorm(res); for (int i = 0; i <= 16; ++i) { const float m = static_cast(i) / 16.0f; - const double got = measuredGain(at(fc, res, m), sr, fc); + const double got = measuredGain(at(fc, res, m, 0.0f, MorphLaw::HighBandLow), sr, fc); if (!(std::fabs(got / q - 1.0) <= kAgreement)) { std::printf("FAIL line %d: morph %.4f res %.1f corner gain %.6f, expected Q " "%.6f (%.3f%%)\n", @@ -345,38 +220,163 @@ static void testCornerMagnitudeIsFlatAcrossTheWholeMorphSweep() { } } -// Continuity as a control, not just at the corner: no step between adjacent morph positions at -// any fixed frequency. A coefficient switch at the centre — the thing an enum would have forced — -// shows up here as a jump. -static void testMorphSweepHasNoDiscontinuity() { - const double sr = 48000.0, fc = 1000.0; - constexpr int kSteps = 40; - for (float res : {0.0f, 0.5f, 1.0f}) { - for (double f : {250.0, 1000.0, 4000.0}) { - double prev = -1.0; - for (int i = 0; i <= kSteps; ++i) { - const float m = static_cast(i) / kSteps; - const double got = measuredGain(at(fc, res, m), sr, f); - if (prev >= 0.0) { - // Scaled by the response's own magnitude at this setting — the passband is - // unity and the corner is Q, so below Q=1 the passband is what a step has to - // be small against, not Q. - const double scale = std::fmax(1.0, filterQFromNorm(res)); - // One step is 1/40 of the travel; the steepest leg moves well under a tenth - // of that scale over one step (measured worst case is 0.03). - const double jump = std::fabs(got - prev) / scale; - if (!(jump < 0.1)) { - std::printf("FAIL line %d: morph %.4f res %.1f at %.0f Hz jumps %.4f\n", - __LINE__, m, res, f, jump); - ++g_fail; - } +// The SEM's centre is a genuine null, not merely a dip: the corner magnitude falls to the float +// noise floor because HP and LP sit at exactly +90 and -90 degrees there, so equal weights cancel +// by construction. Measured worst case across this whole grid is 3.8e-05 (-88 dB); the typical +// figure is -110 to -145 dB. The settle window has to clear the resonator's ring-down before the +// residual means anything — at 0.15 s and Q=10 the leftover transient alone reads as -52 dB and +// would be mistaken for the floor. +static void testHighNotchLowCentreIsATrueNullAtTheCorner() { + for (int r = 0; r < kRateCount; ++r) { + for (double fc : {250.0, 1000.0, 4000.0}) { + for (float res : {0.0f, 0.5f, 1.0f}) { + const double got = measuredGain(at(fc, res, kCentre, 0.0f, MorphLaw::HighNotchLow), + kRates[r], fc, 0.25, 2.0, 0.5); + if (!(got < 2e-4)) { + std::printf("FAIL line %d: SEM notch at sr %.0f fc %.0f res %.1f is %.3e " + "(%.1f dB) — not a null\n", + __LINE__, kRates[r], fc, res, got, + 20.0 * std::log10(got + 1e-300)); + ++g_fail; } - prev = got; } } } } +// The null sits AT the cutoff, not merely somewhere nearby: the response falls monotonically into +// fc from both sides and is orders of magnitude below its own immediate neighbours. At fc=1 kHz, +// +/-5% off the notch already reads -20 dB while the notch itself reads -127 dB. +static void testHighNotchLowNullIsLocatedAtTheCutoff() { + const double sr = 48000.0, fc = 1000.0; + for (float res : {0.0f, 0.5f, 1.0f}) { + const FilterSettings fs = at(fc, res, kCentre, 0.0f, MorphLaw::HighNotchLow); + const double below[] = {0.5, 0.8, 0.95}; + double prev = 1e30; + for (double ratio : below) { + const double got = measuredGain(fs, sr, fc * ratio); + CHECK(got < prev); + prev = got; + } + const double atCorner = measuredGain(fs, sr, fc, 0.25, 2.0, 0.5); + CHECK(atCorner < prev); + + prev = atCorner; + for (double ratio : {1.05, 1.25, 2.0}) { + const double got = measuredGain(fs, sr, fc * ratio); + CHECK(got > prev); + prev = got; + } + // Against its own immediate neighbours, so this is a null rather than a broad scoop. + CHECK(atCorner < 1e-3 * measuredGain(fs, sr, fc * 0.95)); + } +} + +// The SEM's zero is AT the notch frequency, not a broadband level sag: away from the corner the +// two taps are still an equal-power pair, so the sweep holds constant power on its legs. Measured +// deep in each tap's own passband — 50 Hz for the low tap, 20 kHz for the high tap, both far from +// a 1 kHz corner — and divided by that tap's OWN analytic response there, so what is left is the +// weight the law applied. That normalization is load-bearing, not cosmetic: at Q = 0.1 a 2-pole +// approaches its passband so slowly that the pure low tap still reads 0.896 at 50 Hz, and a raw +// reading would report a 20% "sag" that is the Q, not the morph. A LINEAR crossfade would give +// 0.5 at the centre instead of 1.0, so this tolerance discriminates equal-power from linear +// decisively rather than merely confirming a plausible shape. +static void testHighNotchLowLegsHoldConstantPowerAwayFromTheNotch() { + const double sr = 48000.0, fc = 1000.0; + for (float res : {0.0f, 0.5f, 1.0f}) { + const double q = filterQFromNorm(res); + const double lowRef = analyticMag(kLowPass, 50.0, fc, q, sr); + const double highRef = analyticMag(kHighPass, 20000.0, fc, q, sr); + for (int i = 0; i <= 8; ++i) { + const float m = static_cast(i) / 8.0f; + const FilterSettings fs = at(fc, res, m, 0.0f, MorphLaw::HighNotchLow); + const double low = measuredGain(fs, sr, 50.0) / lowRef; + const double high = measuredGain(fs, sr, 20000.0) / highRef; + const double power = low * low + high * high; + if (!(std::fabs(power - 1.0) <= 0.02)) { + std::printf("FAIL line %d: SEM morph %.3f res %.1f leg power %.6f (low %.6f, " + "high %.6f) — expected 1.0\n", + __LINE__, m, res, power, low, high); + ++g_fail; + } + } + } +} + +// Continuity as a control, not just at the corner: no step between adjacent morph positions at +// any fixed frequency, under either law. A coefficient switch at the centre — the thing an enum +// over TOPOLOGIES would have forced — shows up here as a jump. Measured off the SEM's notch +// frequency, since the null itself is a legitimate near-step in the response. +static void testMorphSweepHasNoDiscontinuity() { + const double sr = 48000.0, fc = 1000.0; + constexpr int kSteps = 40; + for (MorphLaw law : kBothLaws) { + for (float res : {0.0f, 0.5f, 1.0f}) { + for (double f : {250.0, 1000.0, 4000.0}) { + if (f == fc && law == MorphLaw::HighNotchLow) continue; + double prev = -1.0; + for (int i = 0; i <= kSteps; ++i) { + const float m = static_cast(i) / kSteps; + const double got = measuredGain(at(fc, res, m, 0.0f, law), sr, f); + if (prev >= 0.0) { + // Scaled by the response's own magnitude at this setting — the passband is + // unity and the corner is Q, so below Q=1 the passband is what a step has + // to be small against, not Q. + const double scale = std::fmax(1.0, filterQFromNorm(res)); + // One step is 1/40 of the travel; the steepest leg moves well under a + // tenth of that scale over one step (measured worst case is 0.03). + const double jump = std::fabs(got - prev) / scale; + if (!(jump < 0.1)) { + std::printf("FAIL line %d: %s morph %.4f res %.1f at %.0f Hz jumps " + "%.4f\n", + __LINE__, lawName(law), m, res, f, jump); + ++g_fail; + } + } + prev = got; + } + } + } + } +} + +// The law selects a MIX, computed once per prepare(); it must not reach the coefficient solve at +// all. Asserted bit-exactly rather than by tolerance — the cutoff, the damping term, and the +// zero-delay-loop solution are the same floats under either law, so no cutoff/Q/rate behaviour +// can differ between them by construction. +static void testMorphLawDoesNotDisturbTheCoefficients() { + for (int r = 0; r < kRateCount; ++r) { + for (int ci = 0; ci <= 8; ++ci) { + for (float res : {0.0f, 0.5f, 1.0f}) { + for (int mi = 0; mi <= 4; ++mi) { + VoiceFilter band, sem; + const float m = mi / 4.0f; + band.prepare({ci / 8.0f, res, m, 0.5f, MorphLaw::HighBandLow}, kRates[r]); + sem.prepare({ci / 8.0f, res, m, 0.5f, MorphLaw::HighNotchLow}, kRates[r]); + const SvfCoeffs& a = band.coeffs(); + const SvfCoeffs& b = sem.coeffs(); + CHECK(a.g == b.g && a.k == b.k); + CHECK(a.a1 == b.a1 && a.a2 == b.a2 && a.a3 == b.a3); + } + } + } + } +} + +// The default is the reviewed-and-measured law, not the SEM leg. The editor and any persisted- +// state codec read this default, so a preset saved before the selector existed must still sound +// exactly as it did — asserted on the folded mix, which is the only thing the kernel sees. +static void testFilterSettingsDefaultsToTheHighBandLowLaw() { + CHECK(FilterSettings{}.morphLaw == MorphLaw::HighBandLow); + + VoiceFilter defaulted, explicitLaw; + defaulted.prepare({0.5f, 0.5f, kCentre, 0.0f}, 48000.0); + explicitLaw.prepare({0.5f, 0.5f, kCentre, 0.0f, MorphLaw::HighBandLow}, 48000.0); + CHECK(defaulted.mix().m0 == explicitLaw.mix().m0); + CHECK(defaulted.mix().m1 == explicitLaw.mix().m1); + CHECK(defaulted.mix().m2 == explicitLaw.mix().m2); +} + // --------------------------------------------------------------------------- // Drive // --------------------------------------------------------------------------- @@ -411,37 +411,41 @@ static float nextNoise(unsigned& rng) { static_cast(1 << 22); } +static void checkDriveZeroBitIdentity(float morph, MorphLaw law) { + VoiceFilter f; + f.prepare(at(1000.0, 1.0f, morph, 0.0f, law), 48000.0); + f.reset(); + LinearKernelRef ref0{f.coeffs(), f.mix()}; + LinearKernelRef ref1{f.coeffs(), f.mix()}; + + unsigned rng0 = 0x13579bdfu; + for (int i = 0; i < 4096; ++i) { + const float x = nextNoise(rng0); + CHECK(f.process(0, x) == ref0.step(x)); + } + + // process(1, ...): channel 1's state is independent of channel 0's above. + unsigned rng1 = 0x2468acefu; + for (int i = 0; i < 4096; ++i) { + const float x = nextNoise(rng1); + CHECK(f.process(1, x) == ref1.step(x)); + } + + // processFrame(): both channels advanced together through the frame entry point, + // continuing from the state each channel already has. + for (int i = 0; i < 4096; ++i) { + float frame[2] = {nextNoise(rng0), nextNoise(rng1)}; + const float want0 = ref0.step(frame[0]); + const float want1 = ref1.step(frame[1]); + f.processFrame(frame, 2); + CHECK(frame[0] == want0); + CHECK(frame[1] == want1); + } +} + static void testDriveZeroIsBitIdenticalToTheLinearKernel() { - for (float morph : {kHighPass, kBandPass, kLowPass}) { - VoiceFilter f; - f.prepare(at(1000.0, 1.0f, morph, 0.0f), 48000.0); - f.reset(); - LinearKernelRef ref0{f.coeffs(), f.mix()}; - LinearKernelRef ref1{f.coeffs(), f.mix()}; - - unsigned rng0 = 0x13579bdfu; - for (int i = 0; i < 4096; ++i) { - const float x = nextNoise(rng0); - CHECK(f.process(0, x) == ref0.step(x)); - } - - // process(1, ...): channel 1's state is independent of channel 0's above. - unsigned rng1 = 0x2468acefu; - for (int i = 0; i < 4096; ++i) { - const float x = nextNoise(rng1); - CHECK(f.process(1, x) == ref1.step(x)); - } - - // processFrame(): both channels advanced together through the frame entry point, - // continuing from the state each channel already has. - for (int i = 0; i < 4096; ++i) { - float frame[2] = {nextNoise(rng0), nextNoise(rng1)}; - const float want0 = ref0.step(frame[0]); - const float want1 = ref1.step(frame[1]); - f.processFrame(frame, 2); - CHECK(frame[0] == want0); - CHECK(frame[1] == want1); - } + for (MorphLaw law : kBothLaws) { + for (float morph : {kHighPass, kBandPass, kLowPass}) checkDriveZeroBitIdentity(morph, law); } } @@ -478,20 +482,22 @@ static void testFullDriveStaysBoundedAtEveryCutoffResonanceAndRate() { for (int r = 0; r < kRateCount; ++r) { const double sr = kRates[r]; - for (int ci = 0; ci <= 8; ++ci) { - for (int mi = 0; mi <= 4; ++mi) { - for (float res : {0.0f, 0.5f, 1.0f}) { - VoiceFilter f; - f.prepare({ci / 8.0f, res, mi / 4.0f, 1.0f}, sr); - f.reset(); - for (int i = 0; i < 4000; ++i) { - const float y = f.process(0, noise()); - if (!std::isfinite(y) || std::fabs(y) > 8.0f) { - std::printf("FAIL line %d: sr=%.0f cutoff=%.2f morph=%.2f res=%.1f " - "full drive produced %g\n", - __LINE__, sr, ci / 8.0, mi / 4.0, res, y); - ++g_fail; - return; + for (MorphLaw law : kBothLaws) { + for (int ci = 0; ci <= 8; ++ci) { + for (int mi = 0; mi <= 4; ++mi) { + for (float res : {0.0f, 0.5f, 1.0f}) { + VoiceFilter f; + f.prepare({ci / 8.0f, res, mi / 4.0f, 1.0f, law}, sr); + f.reset(); + for (int i = 0; i < 4000; ++i) { + const float y = f.process(0, noise()); + if (!std::isfinite(y) || std::fabs(y) > 8.0f) { + std::printf("FAIL line %d: %s sr=%.0f cutoff=%.2f morph=%.2f " + "res=%.1f full drive produced %g\n", + __LINE__, lawName(law), sr, ci / 8.0, mi / 4.0, res, y); + ++g_fail; + return; + } } } } @@ -505,16 +511,18 @@ static void testFullDriveStaysBoundedAtEveryCutoffResonanceAndRate() { static void testFullDriveDoesNotSelfOscillate() { for (int r = 0; r < kRateCount; ++r) { const double sr = kRates[r]; - for (float morph : {kHighPass, kBandPass, kLowPass}) { - VoiceFilter f; - f.prepare(at(1000.0, 1.0f, morph, 1.0f), sr); - f.reset(); - const int excite = static_cast(sr * 0.01); - for (int i = 0; i < excite; ++i) { - f.process(0, static_cast(std::sin(2.0 * kPi * 1000.0 * i / sr))); + for (MorphLaw law : kBothLaws) { + for (float morph : {kHighPass, kBandPass, kLowPass}) { + VoiceFilter f; + f.prepare(at(1000.0, 1.0f, morph, 1.0f, law), sr); + f.reset(); + const int excite = static_cast(sr * 0.01); + for (int i = 0; i < excite; ++i) { + f.process(0, static_cast(std::sin(2.0 * kPi * 1000.0 * i / sr))); + } + for (int i = 0; i < static_cast(sr * 0.5); ++i) f.process(0, 0.0f); + CHECK(f.isSilent()); } - for (int i = 0; i < static_cast(sr * 0.5); ++i) f.process(0, 0.0f); - CHECK(f.isSilent()); } } } @@ -621,198 +629,20 @@ static void testLowCutoffHighRateCornerHoldsTheAnalyticPeak() { } } -// --------------------------------------------------------------------------- -// Stability, denormals, and state -// --------------------------------------------------------------------------- - -static void testFullRangeCutoffSweepAtAudioRateStaysBounded() { - unsigned rng = 0x13579bdfu; - auto noise = [&rng]() { - rng = rng * 1664525u + 1013904223u; - return static_cast(static_cast(rng >> 9) - (1 << 22)) / - static_cast(1 << 22); - }; - - for (int r = 0; r < kRateCount; ++r) { - const double sr = kRates[r]; - for (float morph : {kHighPass, kBandPass, kLowPass}) { - for (float res : {0.0f, 1.0f}) { - for (float drive : {0.0f, 1.0f}) { - VoiceFilter f; - f.reset(); - // A fixed WALL-CLOCK sweep: the same cutoff travel per second at every rate, - // so the per-sample coefficient step gets no gentler as the rate rises. - const int n = static_cast(sr * 0.25); - for (int i = 0; i < n; ++i) { - const float t = static_cast(i) / static_cast(n - 1); - f.prepare({t, res, morph, drive}, 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 - } - } - } - } - } -} - -// The flush tests the ENVELOPE — both integrators — not one sample. ic1 and ic2 are in -// quadrature, so a resonator swings each through zero twice a cycle; flushing on a single one -// injects a step in phase with the resonance, which the resonance amplifies, and the filter -// limit-cycles at the floor forever instead of going quiet. Re-verified for TPT rather than -// assumed to carry over from the retired Direct Form I state. -static void testStateFlushesToZeroWithoutStallingInDenormals() { - for (int r = 0; r < kRateCount; ++r) { - const double sr = kRates[r]; - // The decay to the floor is a fixed WALL-CLOCK time, so the budget scales with the rate. - const int budget = static_cast(sr * 0.5); - for (float morph : {kHighPass, kBandPass, kLowPass}) { - for (float drive : {0.0f, 1.0f}) { - VoiceFilter f; - f.prepare(at(1000.0, 1.0f, morph, drive), sr); - f.reset(); - - // Excite, then hard-cut to silence the way a released voice does. - const int excite = static_cast(sr * 0.01); - for (int i = 0; i < excite; ++i) { - f.process(0, 0.5f * static_cast(std::sin(2.0 * kPi * 1000.0 * i / sr))); - } - - int subnormalSamples = 0, silentAt = -1; - for (int i = 0; i < budget; ++i) { - f.process(0, 0.0f); - const VoiceFilter::State& s = f.state(0); - if ((s.ic1 != 0.0f && std::fabs(s.ic1) < FLT_MIN) || - (s.ic2 != 0.0f && std::fabs(s.ic2) < FLT_MIN)) { - ++subnormalSamples; - } - 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 < budget); - // And it stays silent — a flush that perturbs the loop would re-excite it. - for (int i = 0; i < 1000; ++i) CHECK(f.process(0, 0.0f) == 0.0f); - CHECK(f.isSilent()); - } - } - } -} - -// A high-pass under sustained DC must settle to zero and STAY there. Sampling only the final -// value is not enough: a resonator swings through zero twice a cycle, so a single late sample -// can land near zero while the envelope still rings well above it. This regressed a click train -// on the retired topology, where flushing the FIR history discarded the pinned DC and the next -// sample recomputed a full-amplitude step. TPT has no FIR history to discard, so the hazard is -// structural rather than a tuning — but the assertion is cheap and pins the outcome. -static void testHighPassSustainedDCDoesNotReRing() { - for (int r = 0; r < kRateCount; ++r) { - const double sr = kRates[r]; - for (float drive : {0.0f, 1.0f}) { - VoiceFilter f; - f.prepare(at(1000.0, 1.0f, kHighPass, drive), sr); - f.reset(); - const int settle = static_cast(sr * 0.05); - float worstAfterSettle = 0.0f; - for (int i = 0; i < static_cast(sr * 0.5); ++i) { - const float y = f.process(0, 1.0f); - if (i >= settle) worstAfterSettle = std::fmax(worstAfterSettle, std::fabs(y)); - } - CHECK(worstAfterSettle < 1e-3f); - } - } -} - -static void testImpulseResponseMatchesTheKernel() { - VoiceFilter f; - f.prepare(at(1000.0, 0.5f, kLowPass), 48000.0); - f.reset(); - const SvfCoeffs c = f.coeffs(); - // From a cleared state the first sample reduces to the coefficients alone: v1 == a2, v2 == a3. - CHECK_NEAR(f.process(0, 1.0f), c.a3, 1e-7); - - VoiceFilter bp; - bp.prepare(at(1000.0, 0.5f, kBandPass), 48000.0); - bp.reset(); - CHECK_NEAR(bp.process(0, 1.0f), c.a2, 1e-7); - - VoiceFilter hp; - hp.prepare(at(1000.0, 0.5f, kHighPass), 48000.0); - hp.reset(); - CHECK_NEAR(hp.process(0, 1.0f), 1.0 - c.k * c.a2 - c.a3, 1e-7); -} - -static void testLowpassStepSettlesToUnityAndHighpassRejectsDC() { - const double sr = 48000.0; - VoiceFilter f; - f.prepare(at(1000.0, 0.0f, kLowPass), 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(at(1000.0, 0.0f, kHighPass), sr); - hp.reset(); - float worstAfterSettle = 0.0f; - for (int i = 0; i < 48000; ++i) { - y = hp.process(0, 1.0f); - if (i >= 200) worstAfterSettle = std::fmax(worstAfterSettle, std::fabs(y)); - } - CHECK(worstAfterSettle < 1e-3f); -} - -static void testResetClearsStateButPrepareKeepsIt() { - VoiceFilter f; - f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 48000.0); - f.process(0, 1.0f); - CHECK(!f.isSilent()); - - // A live parameter move must not zero the state — that is what would click. - f.prepare({0.6f, 0.5f, kLowPass, 0.0f}, 48000.0); - CHECK(!f.isSilent()); - f.prepare({0.6f, 0.5f, kBandPass, 1.0f}, 48000.0); - CHECK(!f.isSilent()); - - f.reset(); - CHECK(f.isSilent()); -} - -static void testChannelStateIsIndependent() { - VoiceFilter f; - f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 48000.0); - f.reset(); - f.process(0, 1.0f); - CHECK(f.state(0).ic2 != 0.0f); - CHECK(f.state(1).ic2 == 0.0f); - - float frame[2] = {1.0f, -1.0f}; - f.processFrame(frame, 2); - CHECK(f.state(1).ic2 < 0.0f); - CHECK(frame[0] != frame[1]); -} - int main() { - testCutoffMapsThreeDecadesLogarithmically(); - testCutoffNormRoundTrips(); - testQSpansPointOneToTenWithRootTwoAtCenter(); - testQNormRoundTrips(); - testDriveDepthIsZeroAtRestAndRisesMonotonically(); - testSvfCoefficientsMatchPinnedValues(); testTheSampleRateEntersOnlyThroughG(); testNonPositiveRatePrepareClearsStaleStateAndReportsSilent(); testNonPositiveRatePassesSignalThroughAtEveryMorph(); - testMorphEndpointMixesAreExactlyPureTaps(); - testMorphInteriorPointMatchesCosSinNotAnAlternateEqualPowerCurve(); - testMorphNeverBlendsHighAgainstLowPass(); testMorphEndpointsMatchTheAnalyticTwoPoleTargets(); - testCornerMagnitudeIsFlatAcrossTheWholeMorphSweep(); + testCornerMagnitudeIsFlatAtQAcrossTheHighBandLowSweep(); + testHighNotchLowCentreIsATrueNullAtTheCorner(); + testHighNotchLowNullIsLocatedAtTheCutoff(); + testHighNotchLowLegsHoldConstantPowerAwayFromTheNotch(); testMorphSweepHasNoDiscontinuity(); + testMorphLawDoesNotDisturbTheCoefficients(); + testFilterSettingsDefaultsToTheHighBandLowLaw(); testDriveZeroIsBitIdenticalToTheLinearKernel(); testDriveZeroResponseIsLevelInvariant(); @@ -824,14 +654,6 @@ int main() { testResponseIsRateInvariantAtEveryMorph(); testLowCutoffHighRateCornerHoldsTheAnalyticPeak(); - testFullRangeCutoffSweepAtAudioRateStaysBounded(); - testStateFlushesToZeroWithoutStallingInDenormals(); - testHighPassSustainedDCDoesNotReRing(); - testImpulseResponseMatchesTheKernel(); - testLowpassStepSettlesToUnityAndHighpassRejectsDC(); - testResetClearsStateButPrepareKeepsIt(); - testChannelStateIsIndependent(); - 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; diff --git a/tests/test_filter_morph.cpp b/tests/test_filter_morph.cpp new file mode 100644 index 0000000..aec6601 --- /dev/null +++ b/tests/test_filter_morph.cpp @@ -0,0 +1,208 @@ +// Standalone tests for the pure morph domain: normalized position -> tap weights under both +// morph laws, and the fold of those weights into the kernel's three multipliers. Algebra only — +// no filter is run here. Interior expectations are derived from the intended law (in radicals, +// so they share not even a trig call with the implementation) rather than read back out of it. +// The MEASURED consequences of each law — HP-BP-LP's flat corner, HP-notch-LP's null — live in +// test_filter.cpp, where a filter is actually driven. + +#include "../src/core/instrument/engine/filter/filter_morph.h" +#include "../src/core/instrument/engine/filter/filter_params.h" + +#include +#include +#include +#include +#include + +using namespace reasampler::instrument::engine::filter; + +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; + +static constexpr float kHighPass = 0.0f; +static constexpr float kCentre = 0.5f; +static constexpr float kLowPass = 1.0f; + +// cos and sin of pi/8, from the half-angle identity in radicals: cos(pi/8) = sqrt((1+cos(pi/4))/2) +// with cos(pi/4) = sqrt(2)/2. No trig call, so nothing here is shared with filter_morph's cos/sin. +static double cosPi8() { return std::sqrt((1.0 + std::sqrt(2.0) / 2.0) / 2.0); } // 0.9238795325 +static double sinPi8() { return std::sqrt((1.0 - std::sqrt(2.0) / 2.0) / 2.0); } // 0.3826834324 + +// --------------------------------------------------------------------------- +// Shared across both laws +// --------------------------------------------------------------------------- + +// HighBandLow is enumerator 0 by design: a zero-initialized field, or one absent from an older +// persisted blob and left default-constructed, must land on the default law rather than the SEM +// leg. A codec written against this enum depends on that. +static void testHighBandLowIsTheZeroEnumerator() { + CHECK(static_cast>(MorphLaw::HighBandLow) == 0); + CHECK(MorphLaw{} == MorphLaw::HighBandLow); +} + +// The endpoints are pure taps EXACTLY under BOTH laws, not to within a rounding of cos/sin — the +// laws differ only in the interior. Asserted on the folded mix, where "pure" is an exact +// statement about three floats. +static void testMorphEndpointMixesAreExactlyPureTapsUnderBothLaws() { + const float k = 1.0f / filterQFromNorm(0.5f); + + for (MorphLaw law : {MorphLaw::HighBandLow, MorphLaw::HighNotchLow}) { + const MorphMix hp = morphMix(morphWeights(kHighPass, law), k); + CHECK(hp.m0 == 1.0f && hp.m1 == -k && hp.m2 == -1.0f); // v0 - k*v1 - v2 + + const MorphMix lp = morphMix(morphWeights(kLowPass, law), k); + CHECK(lp.m0 == 0.0f && lp.m1 == 0.0f && lp.m2 == 1.0f); // v2 + + // Out-of-range clamps to the endpoints rather than extrapolating. + CHECK(morphWeights(-1.0f, law).hp == 1.0f); + CHECK(morphWeights(2.0f, law).lp == 1.0f); + } + + // Only the centre differs: a band-pass under one law, an HP+LP sum under the other. + const MorphMix bp = morphMix(morphWeights(kCentre, MorphLaw::HighBandLow), k); + CHECK(bp.m0 == 0.0f && bp.m1 == 1.0f && bp.m2 == 0.0f); // v1 + + const MorphMix notch = morphMix(morphWeights(kCentre, MorphLaw::HighNotchLow), k); + CHECK(notch.m0 != 0.0f && notch.m1 != 0.0f); +} + +// NaN clamps to neither endpoint (every comparison against it is false) and lands on each law's +// degenerate — no crash, a sane fallback rather than an extrapolation. HighNotchLow has no band +// tap to fall back to, so it lands on the leg-zero endpoint instead. +static void testNaNFallsBackToASaneTapPerLaw() { + const float nan = std::numeric_limits::quiet_NaN(); + + const MorphWeights band = morphWeights(nan, MorphLaw::HighBandLow); + CHECK(band.hp == 0.0f && band.bp == 1.0f && band.lp == 0.0f); + + const MorphWeights sem = morphWeights(nan, MorphLaw::HighNotchLow); + CHECK(sem.hp == 1.0f && sem.bp == 0.0f && sem.lp == 0.0f); +} + +// --------------------------------------------------------------------------- +// HighBandLow — adjacent taps only +// --------------------------------------------------------------------------- + +// Pins the cos/sin curve at an interior point, not just the endpoints and the quadrature +// identity (hp^2+bp^2+lp^2=1, which any equal-power reparameterization would also satisfy). +// theta=0.5*pi*t^2 (quadratic in the leg fraction, still equal-power, still exact at both +// ends) would give hp=0.9239/bp=0.3827 here instead of the cos/sin pair's 0.7071/0.7071. +static void testHighBandLowInteriorMatchesCosSinNotAnAlternateEqualPowerCurve() { + const MorphWeights w = morphWeights(0.25f, MorphLaw::HighBandLow); // HP->BP leg, t = 0.5 + const double theta = 0.5 * kPi * 0.5; + CHECK_NEAR(w.hp, std::cos(theta), 1e-6); + CHECK_NEAR(w.bp, std::sin(theta), 1e-6); + CHECK(w.lp == 0.0f); + + // Each leg is half the sweep, so a leg reaches at 0.125 what the SEM's single crossfade + // reaches at 0.25 — the crispest algebraic statement of how the two laws differ. + const MorphWeights eighth = morphWeights(0.125f, MorphLaw::HighBandLow); + CHECK_NEAR(eighth.hp, cosPi8(), 1e-6); + CHECK_NEAR(eighth.bp, sinPi8(), 1e-6); +} + +// HP and LP never carry weight at the same time under THIS law. That is what keeps its centre a +// band-pass: the two are antiphase at the corner and would otherwise cancel into a notch. This +// assertion is law-specific and is deliberately inverted for HighNotchLow below — do not relax +// it to cover both, which would give up the guarantee entirely. +static void testHighBandLowNeverBlendsHighAgainstLowPass() { + for (int i = 0; i <= 200; ++i) { + const MorphWeights w = morphWeights(static_cast(i) / 200.0f, MorphLaw::HighBandLow); + CHECK(w.hp == 0.0f || w.lp == 0.0f); + CHECK(w.hp >= 0.0f && w.bp >= 0.0f && w.lp >= 0.0f); + // Equal power: the active pair sums in quadrature to unity. + CHECK_NEAR(w.hp * w.hp + w.bp * w.bp + w.lp * w.lp, 1.0, 1e-6); + } +} + +// --------------------------------------------------------------------------- +// HighNotchLow — HP against LP, which is the whole mechanism +// --------------------------------------------------------------------------- + +// The exact inverse of the HighBandLow assertion above: blending HP against LP is not a defect +// to be avoided here, it is what produces the notch. The band tap is silent throughout. +static void testHighNotchLowBlendsHighAgainstLowPassWithNoBandTap() { + for (int i = 0; i <= 200; ++i) { + const float n = static_cast(i) / 200.0f; + const MorphWeights w = morphWeights(n, MorphLaw::HighNotchLow); + CHECK(w.bp == 0.0f); + CHECK(w.hp >= 0.0f && w.lp >= 0.0f); + // Both taps carry weight everywhere strictly between the endpoints. + if (i > 0 && i < 200) CHECK(w.hp > 0.0f && w.lp > 0.0f); + // Equal power, which is what keeps the legs from sagging away from the notch frequency. + CHECK_NEAR(w.hp * w.hp + w.lp * w.lp, 1.0, 1e-6); + } +} + +// Pins the single equal-power crossfade at interior points against radical-derived values, so a +// law that still hits both endpoints but bends differently between them fails. Discriminators at +// n=0.25: a LINEAR crossfade gives 0.75/0.25; a two-leg construction (HighBandLow's spacing +// applied to an HP/LP pair) gives 0.7071/0.7071. Both are far outside this tolerance. +static void testHighNotchLowInteriorWeightsMatchTheSingleEqualPowerCrossfade() { + // cos/sin of pi/8 and 3pi/8; the latter pair is the former swapped. + const double c8 = cosPi8(), s8 = sinPi8(); + CHECK_NEAR(c8, 0.9238795325112867, 1e-15); + CHECK_NEAR(s8, 0.3826834323650898, 1e-15); + + const MorphWeights quarter = morphWeights(0.25f, MorphLaw::HighNotchLow); + CHECK_NEAR(quarter.hp, c8, 1e-6); + CHECK_NEAR(quarter.lp, s8, 1e-6); + + const MorphWeights threeQuarters = morphWeights(0.75f, MorphLaw::HighNotchLow); + CHECK_NEAR(threeQuarters.hp, s8, 1e-6); + CHECK_NEAR(threeQuarters.lp, c8, 1e-6); + + const MorphWeights centre = morphWeights(kCentre, MorphLaw::HighNotchLow); + CHECK_NEAR(centre.hp, std::sqrt(2.0) / 2.0, 1e-6); + CHECK_NEAR(centre.lp, std::sqrt(2.0) / 2.0, 1e-6); + + // Symmetric about the centre, so the sweep reads the same in either direction. + for (int i = 0; i <= 100; ++i) { + const float n = static_cast(i) / 100.0f; + const MorphWeights a = morphWeights(n, MorphLaw::HighNotchLow); + const MorphWeights b = morphWeights(1.0f - n, MorphLaw::HighNotchLow); + CHECK_NEAR(a.hp, b.lp, 1e-6); + } +} + +// The centre's cancellation is STRUCTURAL, not a runtime near-miss of two large numbers. The fold +// is m2 = lp - hp, and at the centre the two weights are the same float — cos and sin of pi/4 +// differ by about an ulp of DOUBLE, ~1e-16, which is nine orders below float's ~6e-8 spacing +// there, so they round to one value. m2 is therefore exactly 0 and the output reduces to +// hp*(x - k*v1): the high and low taps cannot drift apart by a rounding. +static void testHighNotchLowCentreFoldsToAnExactlyCancellingMix() { + for (float res : {0.0f, 0.5f, 1.0f}) { + const float k = 1.0f / filterQFromNorm(res); + const MorphWeights w = morphWeights(kCentre, MorphLaw::HighNotchLow); + CHECK(w.hp == w.lp); + + const MorphMix m = morphMix(w, k); + CHECK(m.m2 == 0.0f); + CHECK(m.m0 == w.hp); + CHECK(m.m1 == -w.hp * k); + } +} + +int main() { + testHighBandLowIsTheZeroEnumerator(); + testMorphEndpointMixesAreExactlyPureTapsUnderBothLaws(); + testNaNFallsBackToASaneTapPerLaw(); + + testHighBandLowInteriorMatchesCosSinNotAnAlternateEqualPowerCurve(); + testHighBandLowNeverBlendsHighAgainstLowPass(); + + testHighNotchLowBlendsHighAgainstLowPassWithNoBandTap(); + testHighNotchLowInteriorWeightsMatchTheSingleEqualPowerCrossfade(); + testHighNotchLowCentreFoldsToAnExactlyCancellingMix(); + + if (g_fail == 0) std::printf("filter_morph_tests: all passed\n"); + else std::printf("filter_morph_tests: %d FAILED\n", g_fail); + return g_fail == 0 ? 0 : 1; +} diff --git a/tests/test_filter_params.cpp b/tests/test_filter_params.cpp new file mode 100644 index 0000000..4173098 --- /dev/null +++ b/tests/test_filter_params.cpp @@ -0,0 +1,121 @@ +// Standalone tests for the filter's control domain — normalized knob position to cutoff Hz, Q, +// and drive depth, plus the exact inverses. No DSP is run here and no sample rate appears, which +// is the point: filter_params is deliberately rate-free. Interior points are pinned against a +// derivation of the intended law written out in-test, so a curve that still hits the anchors but +// bends differently between them fails. + +#include "../src/core/instrument/engine/filter/filter_params.h" + +#include +#include + +using namespace reasampler::instrument::engine::filter; + +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 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); + + 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(filterNormFromCutoffHz(1.0f) == 0.0f); + CHECK(filterNormFromCutoffHz(48000.0f) == 1.0f); +} + +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); + + // Pins the single quadratic-in-log-Q curve at two interior points, derived independently by + // solving log Q = a + b*n + c*n^2 through the three anchors above rather than read out of + // the implementation. A two-spliced-log-segments curve (log-linear on each half, the design + // this module doc explicitly rejects for its center-detent slope kink) would give 0.376 and + // 3.761 here instead — both comfortably outside this tolerance. + { + const double lo = std::log(static_cast(kFilterQMin)); + const double mid = std::log(static_cast(kFilterQCenter)); + const double hi = std::log(static_cast(kFilterQMax)); + const double c = 2.0 * lo + 2.0 * hi - 4.0 * mid; + const double b = hi - lo - c; + const double a = lo; + auto qLaw = [&](double n) { return std::exp(a + b * n + c * n * n); }; + CHECK_NEAR(filterQFromNorm(0.25f), qLaw(0.25), 1e-5); + CHECK_NEAR(filterQFromNorm(0.75f), qLaw(0.75), 1e-5); + } + + // 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); +} + +static void testDriveDepthIsZeroAtRestAndRisesMonotonically() { + // Exactly zero, not nearly: the limiter is the identity only at depth 0. + CHECK(filterDriveDepthFromNorm(0.0f) == 0.0f); + CHECK(filterDriveDepthFromNorm(-1.0f) == 0.0f); + CHECK_NEAR(filterDriveDepthFromNorm(1.0f), kFilterDriveDepthMax, 1e-6); + CHECK_NEAR(filterDriveDepthFromNorm(2.0f), kFilterDriveDepthMax, 1e-6); + + // Pins the SQUARE law at an interior point, not just the anchors: a linear law would give + // kFilterDriveDepthMax/2 (2.0) here, not kFilterDriveDepthMax/4 (1.0). + CHECK_NEAR(filterDriveDepthFromNorm(0.5f), kFilterDriveDepthMax * 0.25, 1e-6); + + float prev = -1.0f; + for (int i = 0; i <= 100; ++i) { + const float d = filterDriveDepthFromNorm(static_cast(i) / 100.0f); + CHECK(d > prev); + prev = d; + } +} + +int main() { + testCutoffMapsThreeDecadesLogarithmically(); + testCutoffNormRoundTrips(); + testQSpansPointOneToTenWithRootTwoAtCenter(); + testQNormRoundTrips(); + testDriveDepthIsZeroAtRestAndRisesMonotonically(); + + if (g_fail == 0) std::printf("filter_params_tests: all passed\n"); + else std::printf("filter_params_tests: %d FAILED\n", g_fail); + return g_fail == 0 ? 0 : 1; +} diff --git a/tests/test_filter_state.cpp b/tests/test_filter_state.cpp new file mode 100644 index 0000000..91d7d7b --- /dev/null +++ b/tests/test_filter_state.cpp @@ -0,0 +1,247 @@ +// Standalone tests for the running filter's NUMERICAL behaviour and state lifecycle — bounded +// output under a live parameter sweep, the denormal flush, DC handling, the impulse response +// against the coefficients, and reset/prepare/per-channel state rules. Split from test_filter.cpp +// along the one seam that costs nothing: none of these need the frequency-response measurement +// harness, so the analytic reference lives in exactly one file and cannot fork. + +#include "../src/core/instrument/engine/filter/filter_coeffs.h" +#include "../src/core/instrument/engine/filter/filter_morph.h" +#include "../src/core/instrument/engine/filter/filter_params.h" +#include "../src/core/instrument/engine/filter/voice_filter.h" + +#include +#include +#include +#include + +using namespace reasampler::instrument::engine::filter; + +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; + +static constexpr float kHighPass = 0.0f; +static constexpr float kBandPass = 0.5f; +static constexpr float kCentre = 0.5f; +static constexpr float kLowPass = 1.0f; + +static const MorphLaw kBothLaws[] = {MorphLaw::HighBandLow, MorphLaw::HighNotchLow}; + +// The rates the invariance claims are made over. +static const double kRates[] = {44100.0, 48000.0, 88200.0, 96000.0, 192000.0}; +static constexpr int kRateCount = 5; + +static FilterSettings at(double fcHz, float res, float morph, float drive = 0.0f, + MorphLaw law = MorphLaw::HighBandLow) { + return {filterNormFromCutoffHz(static_cast(fcHz)), res, morph, drive, law}; +} + +static void testFullRangeCutoffSweepAtAudioRateStaysBounded() { + unsigned rng = 0x13579bdfu; + auto noise = [&rng]() { + rng = rng * 1664525u + 1013904223u; + return static_cast(static_cast(rng >> 9) - (1 << 22)) / + static_cast(1 << 22); + }; + + for (int r = 0; r < kRateCount; ++r) { + const double sr = kRates[r]; + for (MorphLaw law : kBothLaws) { + for (float morph : {kHighPass, kBandPass, kLowPass}) { + for (float res : {0.0f, 1.0f}) { + for (float drive : {0.0f, 1.0f}) { + VoiceFilter f; + f.reset(); + // A fixed WALL-CLOCK sweep: the same cutoff travel per second at every + // rate, so the per-sample coefficient step gets no gentler as the rate + // rises. + const int n = static_cast(sr * 0.25); + for (int i = 0; i < n; ++i) { + const float t = static_cast(i) / static_cast(n - 1); + f.prepare({t, res, morph, drive, law}, 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 + } + } + } + } + } + } +} + +// The flush tests the ENVELOPE — both integrators — not one sample. ic1 and ic2 are in +// quadrature, so a resonator swings each through zero twice a cycle; flushing on a single one +// injects a step in phase with the resonance, which the resonance amplifies, and the filter +// limit-cycles at the floor forever instead of going quiet. Re-verified for TPT rather than +// assumed to carry over from the retired Direct Form I state. +static void checkFlushGoesSilent(double sr, float morph, float drive, MorphLaw law) { + // The decay to the floor is a fixed WALL-CLOCK time, so the budget scales with the rate. + const int budget = static_cast(sr * 0.5); + VoiceFilter f; + f.prepare(at(1000.0, 1.0f, morph, drive, law), sr); + f.reset(); + + // Excite, then hard-cut to silence the way a released voice does. + const int excite = static_cast(sr * 0.01); + for (int i = 0; i < excite; ++i) { + f.process(0, 0.5f * static_cast(std::sin(2.0 * kPi * 1000.0 * i / sr))); + } + + int subnormalSamples = 0, silentAt = -1; + for (int i = 0; i < budget; ++i) { + f.process(0, 0.0f); + const VoiceFilter::State& s = f.state(0); + if ((s.ic1 != 0.0f && std::fabs(s.ic1) < FLT_MIN) || + (s.ic2 != 0.0f && std::fabs(s.ic2) < FLT_MIN)) { + ++subnormalSamples; + } + 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 < budget); + // And it stays silent — a flush that perturbs the loop would re-excite it. + for (int i = 0; i < 1000; ++i) CHECK(f.process(0, 0.0f) == 0.0f); + CHECK(f.isSilent()); +} + +static void testStateFlushesToZeroWithoutStallingInDenormals() { + for (int r = 0; r < kRateCount; ++r) { + for (MorphLaw law : kBothLaws) { + for (float morph : {kHighPass, kBandPass, kLowPass}) { + for (float drive : {0.0f, 1.0f}) checkFlushGoesSilent(kRates[r], morph, drive, law); + } + } + } +} + +// A high-pass under sustained DC must settle to zero and STAY there. Sampling only the final +// value is not enough: a resonator swings through zero twice a cycle, so a single late sample +// can land near zero while the envelope still rings well above it. This regressed a click train +// on the retired topology, where flushing the FIR history discarded the pinned DC and the next +// sample recomputed a full-amplitude step. TPT has no FIR history to discard, so the hazard is +// structural rather than a tuning — but the assertion is cheap and pins the outcome. Morph 0 is +// the same pure high-pass under either law, so this needs no law loop. +static void testHighPassSustainedDCDoesNotReRing() { + for (int r = 0; r < kRateCount; ++r) { + const double sr = kRates[r]; + for (float drive : {0.0f, 1.0f}) { + VoiceFilter f; + f.prepare(at(1000.0, 1.0f, kHighPass, drive), sr); + f.reset(); + const int settle = static_cast(sr * 0.05); + float worstAfterSettle = 0.0f; + for (int i = 0; i < static_cast(sr * 0.5); ++i) { + const float y = f.process(0, 1.0f); + if (i >= settle) worstAfterSettle = std::fmax(worstAfterSettle, std::fabs(y)); + } + CHECK(worstAfterSettle < 1e-3f); + } + } +} + +static void testImpulseResponseMatchesTheKernel() { + VoiceFilter f; + f.prepare(at(1000.0, 0.5f, kLowPass), 48000.0); + f.reset(); + const SvfCoeffs c = f.coeffs(); + // From a cleared state the first sample reduces to the coefficients alone: v1 == a2, v2 == a3. + CHECK_NEAR(f.process(0, 1.0f), c.a3, 1e-7); + + VoiceFilter bp; + bp.prepare(at(1000.0, 0.5f, kBandPass), 48000.0); + bp.reset(); + CHECK_NEAR(bp.process(0, 1.0f), c.a2, 1e-7); + + VoiceFilter hp; + hp.prepare(at(1000.0, 0.5f, kHighPass), 48000.0); + hp.reset(); + CHECK_NEAR(hp.process(0, 1.0f), 1.0 - c.k * c.a2 - c.a3, 1e-7); + + // Under HighNotchLow the centre's first sample is the SUM of the high and low taps, scaled by + // the shared weight — the same algebra the null rests on, seen one sample in. + VoiceFilter sem; + sem.prepare(at(1000.0, 0.5f, kCentre, 0.0f, MorphLaw::HighNotchLow), 48000.0); + sem.reset(); + const double w = morphWeights(kCentre, MorphLaw::HighNotchLow).hp; + const double highTap = 1.0 - c.k * c.a2 - c.a3; + const double lowTap = c.a3; + CHECK_NEAR(sem.process(0, 1.0f), w * (highTap + lowTap), 1e-6); +} + +static void testLowpassStepSettlesToUnityAndHighpassRejectsDC() { + const double sr = 48000.0; + VoiceFilter f; + f.prepare(at(1000.0, 0.0f, kLowPass), 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(at(1000.0, 0.0f, kHighPass), sr); + hp.reset(); + float worstAfterSettle = 0.0f; + for (int i = 0; i < 48000; ++i) { + y = hp.process(0, 1.0f); + if (i >= 200) worstAfterSettle = std::fmax(worstAfterSettle, std::fabs(y)); + } + CHECK(worstAfterSettle < 1e-3f); +} + +static void testResetClearsStateButPrepareKeepsIt() { + VoiceFilter f; + f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 48000.0); + f.process(0, 1.0f); + CHECK(!f.isSilent()); + + // A live parameter move must not zero the state — that is what would click. Switching the + // morph law is a parameter move like any other: it only recomputes the mix. + f.prepare({0.6f, 0.5f, kLowPass, 0.0f}, 48000.0); + CHECK(!f.isSilent()); + f.prepare({0.6f, 0.5f, kBandPass, 1.0f}, 48000.0); + CHECK(!f.isSilent()); + f.prepare({0.6f, 0.5f, kCentre, 1.0f, MorphLaw::HighNotchLow}, 48000.0); + CHECK(!f.isSilent()); + + f.reset(); + CHECK(f.isSilent()); +} + +static void testChannelStateIsIndependent() { + VoiceFilter f; + f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 48000.0); + f.reset(); + f.process(0, 1.0f); + CHECK(f.state(0).ic2 != 0.0f); + CHECK(f.state(1).ic2 == 0.0f); + + float frame[2] = {1.0f, -1.0f}; + f.processFrame(frame, 2); + CHECK(f.state(1).ic2 < 0.0f); + CHECK(frame[0] != frame[1]); +} + +int main() { + testFullRangeCutoffSweepAtAudioRateStaysBounded(); + testStateFlushesToZeroWithoutStallingInDenormals(); + testHighPassSustainedDCDoesNotReRing(); + testImpulseResponseMatchesTheKernel(); + testLowpassStepSettlesToUnityAndHighpassRejectsDC(); + testResetClearsStateButPrepareKeepsIt(); + testChannelStateIsIndependent(); + + if (g_fail == 0) std::printf("filter_state_tests: all passed\n"); + else std::printf("filter_state_tests: %d FAILED\n", g_fail); + return g_fail == 0 ? 0 : 1; +}