Merge Θ-W1-T3: TPT/SVF filter with HP-BP-LP and HP-notch-LP morph laws and a drive stage

This commit is contained in:
2026-07-30 11:09:44 -04:00
20 changed files with 2067 additions and 398 deletions
+42
View File
@@ -995,6 +995,20 @@ 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 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;
# VoiceFilter::process is header-inline so the 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/filter_morph.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 +1109,34 @@ 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: 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, bounded-output/self-oscillation
# under full drive, and the state lifecycle — none of it needs the measurement harness below.
# 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.
# 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)
# ---------------------------------------------------------------------------
# 4) The REAPER extension — a loadable module (dlopen'd by REAPER, not linked).
# ---------------------------------------------------------------------------
+231
View File
@@ -0,0 +1,231 @@
# src/core/instrument/engine/filter — the per-voice resonant filter
## Scope
The pure per-voice filter a sounding voice runs: a Zavalishin TPT/SVF with a continuous
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
no call site yet to force a collision into the open at compile time. Five files, one
responsibility each:
- `filter_params` — the control domain: normalized [0,1] knob position → cutoff Hz, Q, and
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: `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.
`process()` is defined in the header.
## Invariants
### No vtable on the per-sample path
The Cortex-M4 source this began as 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, `process()` is inlined, and 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.
### The rate enters ONLY through `g = tan(pi*fc/sr)`
There is no reference sample rate, calibration rate, or fallback rate anywhere in this
module, and introducing one is the specific regression to guard against. An earlier design
carried a `kFilterFeedbackDelaySeconds = 1/48000` tuning constant for a feedback tap; that
tap, its ring buffer, and the constant are all deleted. A non-positive rate yields `g == 0`
and a bypass mix (signal passes through) — never an invented rate.
### Why the high-pass feedback tap was right on Q15 hardware and wrong here
The ported firmware fed a saturated share of an earlier output back into the high-pass
input. Its stated rationale — that the HP numerator collapses toward zero at low cutoff,
taking the resonance with it — is **inverted**, and the comment asserting it has been
removed rather than carried forward. Measurement: the HP `b0` approaches **1** as cutoff
falls (0.99987 at 20 Hz); it is the **low-pass** `b0` that collapses (1.7e06 at 20 Hz).
The tap was a Q15 fixed-point workaround. At 16-bit fixed point the low-cutoff biquad loses
a ~17-bit cancellation and the resonance really does die; the feedback injected it back by
another route. float32 survives that cancellation with 7 bits to spare, so on this target
the tap did not restore character — it *reduced* it (HP landed 0.4% off the analytic RBJ
target with the tap disabled, and 25% off with it enabled), and it introduced both level
dependence and rate dependence.
Daniel's ruling on the level-dependent resonance bloom it produced: *"was a feature on the
hardware (one knob colorful HP for master FX), wrong choice for this approach."* Drive is
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, LP at 1.0, continuous
throughout, and both endpoints are exact under either law — only the centre differs.
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.
### 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 on the shipped `{250, 1000, 4000}` Hz cutoff grid, typically 110 to
145 dB. Over the full control range (20 Hz 20 kHz, Q 0.1 10) the worst residual is
shallower — 69.8 dB at 192 kHz / 30 Hz / Q=10 — from float conditioning in the folded
`x k·v1` term as `fc/sr → 1e-4` at high Q; it is Q-dependent (Q=0.1 holds 110 dB everywhere)
and still an excellent notch, not a broadband defect. `test_filter.cpp`'s null test covers this
full range with a Q-scaled threshold rather than the flat 74 dB the shipped grid alone would
justify. The fold makes the centre's cancellation 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
`softLimit(u, depth) = u / sqrt(1 + (depth*u)²)` shapes the **band-pass integrator state**.
Three properties carry the design:
- `depth == 0` makes it algebraically the identity (`x / sqrt(1) == x`, exact in IEEE), so
drive 0 is **bit-exact** linear whether or not `softLimit` is actually called. The test
asserts bit-identity against the same kernel with the limiter deleted.
- `process()` gates the call on `driven_` (`driveDepth_ != 0`, cached at `prepare()`) rather
than calling `softLimit` unconditionally. `sqrt`/div sit on the per-sample recursive
dependency chain, so out-of-order execution can't hide their latency, and at drive 0 that
cost buys nothing. Measured: 11.2 ns/sample unconditional vs 4.1 ns gated — the gated form
lands at the limiter-removed floor. `driven_` only changes at `prepare()`, so the branch
predicts perfectly. The gate is a perf optimization on top of the bit-identity above, not a
substitute for it — deleting the gate would still be correct, just 2.7x slower at rest.
- `|softLimit(u, d)| <= |u|` for every depth, so the state update can only shrink the state.
The filter cannot gain energy from the drive stage: stability at any Q and any cutoff is
structural, and self-oscillation is impossible. This is why the shaper must keep unit slope
at the origin — a shaper with gain above 1 there turns the resonator into an oscillator.
- It shapes the **state**, not the zero-delay loop. A nonlinearity inside the loop would
break the closed-form `a1`/`a2`/`a3` solve and need per-sample Newton iteration.
Placement is the resonance path because that is where the firmware's character came from,
and because the band-pass state sits at zero in the passband and at DC — so drive colours
the resonance and leaves the passband transparent (measured 0.98 at max drive). It is not a
distortion box in series with the signal; a caller wanting that has every other plugin.
**Drive × resonance interact by design.** What reaches the shaper is the resonance state,
already multiplied by roughly `2*Q`, so the same drive setting bites harder the more
resonance is dialled in — and harder on a hotter input. That level dependence is the
*point* of an explicit drive control; what Daniel rejected was level dependence nobody
asked for. At drive 0 there is none, to 0.0004% over a 1000:1 level range.
`kFilterDriveDepthMax` (4.0) was set against measurement, not feel: at max drive, full-scale
input and max resonance the resonant peak lands ~10 dB under the passband — plainly
crushed, which is the asked-for "extreme". Raising it further inverts the filter's shape
(21 dB under passband at depth 64), turning the peak the user dialled in into a notch.
There is deliberately **no makeup gain** — any law for it would be invented rather than
derived, and drive is due an ear pass against the radial dial.
### 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 `svfCoeffs` where the rate
is already a parameter. 20 kHz is under 0.48·sr at 44.1k and above, so the clamp never eats
live knob travel there; the source's hardcoded 23 kHz endpoint did exactly that at 44.1k.
### 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 SVF consumes it as
`k = 1/Q`.
### Denormal flushing: why conjunctive, honestly
`process()` flushes **both** integrators to exact zero once both are below
`kFilterDenormalFloor` (1e-30). The honest reason is narrower than it sounds: `isSilent()`
means "both integrators are exactly zero," so both have to reach zero for that check to mean
anything, and the conjunctive test is the cheapest way to guarantee it.
The stronger claim — that a per-variable flush limit-cycles at the floor — does **not**
reproduce on this topology. Measured (Q=10, fc=1kHz, 48k): shipped conjunctive goes silent at
sample 10783 with 0 subnormals; a per-variable independent flush goes silent ~180 samples
earlier and an either-below-zero-both flush ~970 samples earlier, both also 0 subnormals, no
limit cycle, and the same excited RMS. That claim WAS real on the retired Direct Form I state,
where the flushed variables (`y1`/`y2`) were the actual filter OUTPUT, so zeroing one injected
a discontinuity the resonance then amplified. Here `ic1`/`ic2` are integrator STATE, not
output: zeroing one only removes energy, a contraction rather than an injection, so the hazard
is structurally absent. The only demonstrable hazard is no flush at all, which never reaches
exact zero and grinds through subnormals for thousands of samples on a released voice.
Keep the conjunctive test regardless — it costs nothing extra and is the right guarantee for
`isSilent()` — but don't cite the limit-cycle rationale for TPT; it belongs to the retired
topology.
## Gotchas
- **TPT is what fixed the low-cutoff conditioning defect** — this is a topology change, not
a relocation. Direct Form I encoded pole proximity in `a1 → -2`, `a2 → +1` and cancelled
them against each other every sample; at `fc/sr ≈ 1e-4` that ~17-bit cancellation moved the
measured 20 Hz / 192 kHz LP peak by **-27% on a true-peak scan, -57% measured at the
analytic peak frequency** (the degraded pole itself moves, so the two methods diverge), and
the error is non-monotone with rate rather than a fixed percentage (+5% high at 96 kHz).
TPT encodes the same proximity in `a1`'s small deviation from 1, which float32 resolves:
checked against an exact-double evaluation of the same difference equation (which matches
the analytic target to within measurement noise), TPT's float32-narrowed coefficients are
genuinely ~0.02% low at 48 kHz, widening to ~0.03% low at 192 kHz — real coefficient
narrowing, not measurement-window noise, and comfortably inside the test's 0.4% tolerance
either way. Do not reintroduce a direct-form kernel.
- **`prepare()` deliberately does not clear state** — a live parameter move must glide, not
click. Call `reset()` at note-on. **Exception: the non-positive-rate bypass path.** There,
`a1=1, a2=a3=0` makes both state updates the exact identity and `bypassMix()` never reads
the state at all, so a stale nonzero `ic1`/`ic2` would otherwise latch `isSilent()` false
forever with no audible effect either way — `prepare()` clears state on that path only,
which costs nothing audibly since bypass ignores it.
- **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
budget expressed in samples is therefore itself a rate assumption — a fixed 20000 samples
is ample at 48k and expires mid-decay at 96k and above.
@@ -0,0 +1,40 @@
#include "core/instrument/engine/filter/filter_coeffs.h"
#include <cmath>
namespace reasampler::instrument::engine::filter {
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
SvfCoeffs svfCoeffs(float cutoffHz, float q, double sampleRate) {
const double qq = clampd(q, kFilterQMin, kFilterQMax);
const double k = 1.0 / qq;
double g = 0.0;
if (sampleRate > 0.0) {
const double fc = clampd(cutoffHz, kFilterCutoffMinHz, kFilterNyquistFraction * sampleRate);
g = std::tan(kPi * fc / sampleRate);
}
// Solved in double and narrowed once. The intermediate g*(g+k) is the term that carries the
// pole proximity, so forming it in float would throw away the conditioning TPT just bought.
const double a1 = 1.0 / (1.0 + g * (g + k));
const double a2 = g * a1;
const double a3 = g * a2;
SvfCoeffs c;
c.g = static_cast<float>(g);
c.k = static_cast<float>(k);
c.a1 = static_cast<float>(a1);
c.a2 = static_cast<float>(a2);
c.a3 = static_cast<float>(a3);
return c;
}
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,42 @@
// filter_coeffs.h — Zavalishin topology-preserving-transform state-variable coefficients.
// The rate enters ONLY through g = tan(pi*fc/sr); there is no reference or calibration rate
// anywhere in this module, and reintroducing one would restore the rate-dependent resonance
// the TPT rewrite exists to remove.
#pragma once
#include "core/instrument/engine/filter/filter_params.h"
namespace reasampler::instrument::engine::filter {
// The two-integrator SVF's per-sample constants. a1/a2/a3 are the algebraic solution of the
// zero-delay feedback loop, so the kernel needs no iteration.
struct SvfCoeffs {
float g = 0.0f; // tan(pi*fc/sr) — the ONLY place the sample rate appears
float k = 1.0f; // 1/Q, the damping term
float a1 = 1.0f;
float a2 = 0.0f;
float a3 = 0.0f;
};
// Highest fraction of the sample rate the pre-warp stays well-conditioned at: tan() diverges
// as fc approaches sr/2.
inline constexpr double kFilterNyquistFraction = 0.48;
// cutoffHz is clamped into [kFilterCutoffMinHz, kFilterNyquistFraction*sampleRate] and q into
// [kFilterQMin, kFilterQMax]. A non-positive sampleRate yields g == 0 — we refuse to invent a
// rate rather than assume 44.1k.
//
// Float storage is safe HERE in a way it was not for the retired Direct Form I path. DF1 encoded
// pole proximity in a1 -> -2, a2 -> +1 and cancelled them against each other every sample; at
// fc/sr ~ 1e-4 that ~17-bit cancellation moved the measured 20 Hz/192 kHz LP peak by -27%
// (true-peak scan) to -57% (point measurement at the analytic peak frequency, since the
// degraded pole itself moves) -- and the error is non-monotone with rate, not a fixed percentage
// (+5% high at 96 kHz). TPT encodes the same proximity in a1's small DEVIATION from 1, which
// float resolves: measured against an exact-double evaluation of the same difference equation
// (which matches the analytic target to within measurement noise), TPT's float32-narrowed
// coefficients land genuinely ~0.02% low at 48 kHz, widening to ~0.03% low at 192 kHz -- both
// comfortably inside the test's 0.4% tolerance.
SvfCoeffs svfCoeffs(float cutoffHz, float q, double sampleRate);
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,64 @@
#include "core/instrument/engine/filter/filter_morph.h"
#include <cmath>
namespace reasampler::instrument::engine::filter {
namespace {
constexpr double kPi = 3.14159265358979323846;
struct Pair {
double a, b;
};
// Equal-power crossfade, EXACT at both ends by construction rather than by rounding: cos and sin
// of the leg's quarter turn are only 1e-17 from 0/1 at the endpoints, and the endpoints have to
// be pure taps, not a pure tap plus a -324 dB neighbour.
Pair equalPower(double t) {
if (!(t > 0.0)) return {1.0, 0.0};
if (t >= 1.0) return {0.0, 1.0};
const double theta = 0.5 * kPi * t;
return {std::cos(theta), std::sin(theta)};
}
} // namespace
MorphWeights morphWeights(float norm, MorphLaw law) {
const double n = norm < 0.0 ? 0.0 : (norm > 1.0 ? 1.0 : static_cast<double>(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<float>(p.a);
w.bp = 0.0f;
w.lp = static_cast<float>(p.b);
return w;
}
if (n <= 0.5) {
const Pair p = equalPower(2.0 * n); // HP -> BP
w.hp = static_cast<float>(p.a);
w.bp = static_cast<float>(p.b);
w.lp = 0.0f;
} else {
const Pair p = equalPower(2.0 * n - 1.0); // BP -> LP
w.hp = 0.0f;
w.bp = static_cast<float>(p.a);
w.lp = static_cast<float>(p.b);
}
return w;
}
MorphMix morphMix(const MorphWeights& w, float k) {
MorphMix m;
m.m0 = w.hp;
m.m1 = w.bp - w.hp * k;
m.m2 = w.lp - w.hp;
return m;
}
MorphMix bypassMix() { return MorphMix{1.0f, 0.0f, 0.0f}; }
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,70 @@
// 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 {
// 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, 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) 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.
//
// 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
// path at three multiplies and spares it ever forming the high tap.
struct MorphMix {
float m0 = 0.0f;
float m1 = 0.0f;
float m2 = 1.0f;
};
MorphMix morphMix(const MorphWeights& w, float k);
// Passes the input through untouched, whatever the morph position asks for. Reserved for a
// sample rate we cannot form a filter from: silencing an instrument is a worse failure than
// ignoring the morph, and at g == 0 a low-pass tap is analytically silent.
MorphMix bypassMix();
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,64 @@
#include "core/instrument/engine/filter/filter_params.h"
#include <cmath>
namespace reasampler::instrument::engine::filter {
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<double>(kFilterQMin));
const double mid = std::log(static_cast<double>(kFilterQCenter));
const double hi = std::log(static_cast<double>(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<double>(kFilterCutoffMinHz));
const double hi = std::log(static_cast<double>(kFilterCutoffMaxHz));
return static_cast<float>(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<double>(kFilterCutoffMinHz));
const double hi = std::log(static_cast<double>(kFilterCutoffMaxHz));
return static_cast<float>(clamp01((std::log(static_cast<double>(hz)) - lo) / (hi - lo)));
}
float filterQFromNorm(float norm) {
const QCurve k = qCurve();
const double n = clamp01(norm);
return static_cast<float>(std::exp(k.a + n * (k.b + k.c * n)));
}
float filterDriveDepthFromNorm(float norm) {
const double n = clamp01(norm);
return static_cast<float>(kFilterDriveDepthMax * n * 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<double>(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<float>(clamp01((-k.b + std::sqrt(d)) / (2.0 * k.c)));
}
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,52 @@
// filter_params.h — control-domain mapping for the voice filter: normalized [0,1] knob
// positions to cutoff Hz, Q, and drive depth. 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::filter {
// 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 44.1 kHz and above; below that (e.g. 32 kHz, 22.05 kHz)
// the clamp still handles it correctly, it just eats the top of the knob travel at those rates.
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;
// Depth at the top of the drive control. The limiter's knee is at 1/depth, and the resonance
// swings the state to roughly 2*Q*level, so this is the range over which drive bites. Chosen
// against measurement rather than by feel: at max drive, full-scale input and max resonance the
// resonant peak lands ~10 dB under the passband — plainly crushed, which is the asked-for
// "extreme". Raising it further inverts the filter's shape (measured 21 dB under passband at
// depth 64), turning the peak the user dialled in into a notch.
inline constexpr float kFilterDriveDepthMax = 4.0f;
// 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);
// Drive depth for the in-loop limiter. Square law, not linear: the knee is 1/depth, so a linear
// depth would spend most of the audible travel in the first tenth of the knob. Exactly 0 at
// norm 0 — the limiter is then algebraically the identity, which is what makes drive=0 bit-exact
// linear rather than merely close.
float filterDriveDepthFromNorm(float norm);
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,34 @@
// filter_saturate.h — the drive stage's soft limiter. Header-inline: it sits inside the
// per-voice per-sample recursion.
#pragma once
#include <cmath>
namespace reasampler::instrument::engine::filter {
// Odd, smooth, strictly monotone, bounded by 1/depth, with unit slope at the origin.
//
// Three properties are load-bearing and none of them are tuning:
// - depth == 0 makes this ALGEBRAICALLY the identity (x / sqrt(1) == x, exact in IEEE), so
// drive = 0 is bit-exact linear whether or not the caller special-cases it. (voice_filter.h
// gates the call on drive != 0 anyway, but as a perf optimization, not because correctness
// needs it.)
// - |softLimit(x, d)| <= |x| for every d, so dropping it into the resonance state update can
// only ever shrink the state. The filter therefore cannot gain energy from the drive stage:
// stability at any Q and any cutoff is structural, not a tuned margin, and it can never
// self-oscillate.
// - Unit slope at the origin, so the shaper adds no gain of its own at any depth. What reaches
// it is the resonance state, already multiplied by roughly 2*Q, which is why drive and
// resonance interact: the same drive setting bites harder the more resonance is dialled in.
//
// The retired feedbackSaturate() is deliberately not carried forward: it had 0.75 slope at the
// origin, a fixed +/-2.0 threshold calibrated for firmware excursion levels, and turned over
// (non-monotone) past x = 6. That absolute threshold is the origin of the level-dependent
// resonance this rewrite removes — do not reintroduce it.
inline float softLimit(float x, float depth) {
const float s = depth * x;
return x / std::sqrt(1.0f + s * s);
}
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,33 @@
#include "core/instrument/engine/filter/voice_filter.h"
namespace reasampler::instrument::engine::filter {
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, 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
// was already going to be ignored) and prevents a stale nonzero ic1/ic2 from latching
// isSilent() false forever, which prepare() otherwise deliberately never does.
mix_ = bypassMix();
for (State& s : state_) s = State{};
}
driveDepth_ = filterDriveDepthFromNorm(settings.driveNorm);
driven_ = driveDepth_ != 0.0f;
}
void VoiceFilter::reset() {
for (State& s : state_) s = State{};
}
bool VoiceFilter::isSilent() const {
for (const State& s : state_) {
if (s.ic1 != 0.0f || s.ic2 != 0.0f) return false;
}
return true;
}
} // namespace reasampler::instrument::engine::filter
@@ -0,0 +1,125 @@
// voice_filter.h — per-voice TPT state-variable filter with a continuous HP->BP->LP morph and
// an in-loop drive stage. Concrete type, no vtable: this sits on the per-voice per-sample path,
// so process() is header-inline. No allocation, no virtual dispatch, no I/O in process().
#pragma once
#include <cassert>
#include <type_traits>
#include "core/instrument/engine/filter/filter_coeffs.h"
#include "core/instrument/engine/filter/filter_morph.h"
#include "core/instrument/engine/filter/filter_params.h"
#include "core/instrument/engine/filter/filter_saturate.h"
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, 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
// 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 integrator pair per channel.
static constexpr int kMaxChannels = 2;
struct State {
float ic1 = 0.0f; // band-pass integrator
float ic2 = 0.0f; // low-pass integrator
};
// Recomputes coefficients from the control positions. State 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];
const float v3 = x - s.ic2;
const float v1 = coeffs_.a1 * s.ic1 + coeffs_.a2 * v3;
const float v2 = s.ic2 + coeffs_.a2 * s.ic1 + coeffs_.a3 * v3;
// The drive stage, and the only nonlinearity. It shapes the BAND-PASS integrator state
// rather than the input because that state IS the resonance: in the passband and at DC
// it sits at zero, so drive colours the resonance and leaves the passband transparent.
// Placing it on the state rather than inside the zero-delay loop keeps a1/a2/a3 an exact
// algebraic solve — a nonlinearity inside the loop would need per-sample Newton
// iteration. softLimit is a contraction, so this cannot destabilize the filter.
//
// Gated on driven_ rather than called unconditionally: sqrt and div sit on this
// recursive dependency chain, so out-of-order execution can't hide them, and at drive 0
// (the default) that cost buys nothing — softLimit(x, 0) == x algebraically. Measured:
// 11.2 ns/sample unconditional vs 4.1 ns gated, matching the limiter-removed floor.
// driven_ only changes at prepare(), so the branch predicts perfectly. Bit-identity at
// drive 0 holds either way, by algebra — the gate is a perf optimization, not what makes
// it exact.
const float u = 2.0f * v1 - s.ic1;
s.ic1 = driven_ ? softLimit(u, driveDepth_) : u;
s.ic2 = 2.0f * v2 - s.ic2;
// Snap the state once the whole resonator has decayed past -600 dB. isSilent() means
// "both integrators are exactly zero," so both must reach zero for that check to be
// meaningful — the conjunctive test is the cheapest guarantee of that, not a defense
// against a demonstrated limit cycle on this topology (measured: a per-variable flush
// and an either-below-zero-both flush both go silent here too, no limit cycle, no
// subnormals). That risk was real on the retired Direct Form I state, where a per-sample
// flush zeroed y1/y2 — the actual OUTPUT — injecting a step the resonance then amplified.
// ic1/ic2 are integrator STATE, not output; zeroing one only removes energy, a
// contraction rather than an injection. The only demonstrable hazard here is no flush at
// all, which never reaches exact zero and stalls in subnormals for thousands of samples.
if (s.ic1 > -kFilterDenormalFloor && s.ic1 < kFilterDenormalFloor &&
s.ic2 > -kFilterDenormalFloor && s.ic2 < kFilterDenormalFloor) {
s.ic1 = 0.0f;
s.ic2 = 0.0f;
}
return mix_.m0 * x + mix_.m1 * v1 + mix_.m2 * v2;
}
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 integrator 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 SvfCoeffs& coeffs() const { return coeffs_; }
const MorphMix& mix() const { return mix_; }
private:
SvfCoeffs coeffs_{};
MorphMix mix_{};
float driveDepth_ = 0.0f;
bool driven_ = false; // driveDepth_ != 0, cached so process() branches on a bool, not a float compare
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<VoiceFilter>, "no vtable on the per-sample path");
static_assert(std::is_trivially_copyable_v<VoiceFilter>, "state is plain values, never owned");
} // namespace reasampler::instrument::engine::filter
-73
View File
@@ -1,73 +0,0 @@
#pragma once
#include "util.hpp"
#include "filter.hpp"
#include "filter_params.hpp"
template <int k_channels, typename TUIParams>
class Biqaud : public Filter<k_channels, FeedbackLine, NormalCoefficients, TUIParams, FilterParameters>
{
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 <int k_channels, typename TUIParams>
class BiquadHP : public Biqaud<k_channels, TUIParams>
{
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 <int k_channels, typename TUIParams>
class BiquadLP : public Biqaud<k_channels, TUIParams>
{
public:
BiquadLP(const uint32_t& sample_rate, FilterParameters *params);
NormalCoefficients prepare_coefficients() override;
};
template <int k_channels, typename TUIParams>
class Biquad1PoleLP : public BiquadLP<k_channels, TUIParams>
{
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"
-195
View File
@@ -1,195 +0,0 @@
#pragma once
#include "basicmaths.h"
#include "biquad.hpp"
template <int k_channels, typename TUIParams>
Biqaud<k_channels, TUIParams>::Biqaud(const uint32_t& p_sample_rate, FilterParameters *p_params)
: Filter<k_channels, FeedbackLine, NormalCoefficients, TUIParams, FilterParameters>(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 <int k_channels, typename TUIParams>
void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
void Biqaud<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
BiquadHP<k_channels, TUIParams>::BiquadHP(const uint32_t& p_sample_rate, FilterParameters *p_params)
: Biqaud<k_channels, TUIParams>(p_sample_rate, p_params) {}
template <int k_channels, typename TUIParams>
void BiquadHP<k_channels, TUIParams>::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<k_channels, TUIParams>::process_channel_frame(state, coeff, input, y);
}
template <int k_channels, typename TUIParams>
NormalCoefficients BiquadHP<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
BiquadLP<k_channels, TUIParams>::BiquadLP(const uint32_t& p_sample_rate, FilterParameters *p_params)
: Biqaud<k_channels, TUIParams>(p_sample_rate, p_params) {}
template <int k_channels, typename TUIParams>
NormalCoefficients BiquadLP<k_channels, TUIParams>::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 <int k_channels, typename TUIParams>
Biquad1PoleLP<k_channels, TUIParams>::Biquad1PoleLP(const uint32_t& p_sample_rate, FilterParameters *p_params)
: BiquadLP<k_channels, TUIParams>(p_sample_rate, p_params) {}
template <int k_channels, typename TUIParams>
void Biquad1PoleLP<k_channels, TUIParams>::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<k_channels, TUIParams>::process_channel_frame(state, coeff, resonant_input, y);
}
template <int k_channels, typename TUIParams>
NormalCoefficients Biquad1PoleLP<k_channels, TUIParams>::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;
}
-62
View File
@@ -1,62 +0,0 @@
#pragma once
#include "util.hpp"
// #include "fdecorator.hpp"
template <int k_channels, typename TFeedbackLine, typename TCoefficients, typename TUIParams, typename TFilterParams>
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 <int, typename, typename, typename, typename, typename, typename>
friend class FilterDecorator;
};
template <int k_channels, typename TFeedbackLine, typename TCoefficients, typename TUIParams, typename TFilterParams>
class Filter : public FilterBase<k_channels, TFeedbackLine, TCoefficients, TUIParams, TFilterParams>
{
public:
Filter(TFilterParams *p)
: FilterBase<k_channels, TFeedbackLine, TCoefficients, TUIParams, TFilterParams>(p) {}
};
-22
View File
@@ -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;
-46
View File
@@ -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;
}
+593
View File
@@ -0,0 +1,593 @@
// 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"
#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 <cmath>
#include <cstdio>
#include <initializer_list>
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;
// 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;
// The measurement pass's bar, and the bar the rewrite exists to hold: peak and passband agree
// with the analytic target to better than this at every rate, level, and morph position.
static constexpr double kAgreement = 0.004;
// ---------------------------------------------------------------------------
// Independent references
// ---------------------------------------------------------------------------
// The analog 2-pole prototype |H(jW)| evaluated at the bilinear-warped frequency. The TPT maps
// the digital frequency onto the prototype EXACTLY at the prewarped corner, so this is the exact
// digital magnitude — derived from the continuous-time prototype and the transform rather than
// from anything filter_coeffs computes.
static double analyticMag(float morph, double freq, double fc, double q, double sr) {
const double w = std::tan(kPi * freq / sr) / std::tan(kPi * fc / sr);
const double dRe = 1.0 - w * w, dIm = w / q;
const double den = std::sqrt(dRe * dRe + dIm * dIm);
if (morph == kHighPass) return w * w / den;
if (morph == kBandPass) return w / den;
return 1.0 / den;
}
// Steady-state gain of the running filter at one frequency. Windows are wall-clock, not sample
// counts, so every rate integrates the same amount of signal.
static double measuredGain(const FilterSettings& fs, double sr, double freq, double amp = 0.25,
double settleSec = 0.15, double measureSec = 0.10) {
VoiceFilter f;
f.prepare(fs, sr);
f.reset();
const int settle = static_cast<int>(sr * settleSec);
const int measure = static_cast<int>(sr * measureSec);
double sumSq = 0.0;
for (int i = 0; i < settle + measure; ++i) {
const float y = f.process(0, static_cast<float>(amp * std::sin(2.0 * kPi * freq * i / sr)));
if (i >= settle) sumSq += static_cast<double>(y) * y;
}
return std::sqrt(sumSq / measure) / (amp / std::sqrt(2.0));
}
static FilterSettings at(double fcHz, float res, float morph, float drive = 0.0f,
MorphLaw law = MorphLaw::HighBandLow) {
return {filterNormFromCutoffHz(static_cast<float>(fcHz)), res, morph, drive, law};
}
// ---------------------------------------------------------------------------
// SVF coefficients — pinned literals plus an independent derivation
// ---------------------------------------------------------------------------
static void testSvfCoefficientsMatchPinnedValues() {
const double sr = 48000.0, fc = 1000.0, q = std::sqrt(2.0);
const SvfCoeffs c = svfCoeffs(static_cast<float>(fc), static_cast<float>(q), sr);
// Pinned literals: change the math and these fail.
CHECK_NEAR(c.g, 0.0655434653, 2e-9);
CHECK_NEAR(c.k, 0.7071067691, 2e-9);
CHECK_NEAR(c.a1, 0.9517988563, 2e-9);
CHECK_NEAR(c.a2, 0.0623841919, 2e-9);
CHECK_NEAR(c.a3, 0.0040888758, 2e-9);
// Independent derivation — proves the pins are the TPT solve and not just "what we emit".
const double g = std::tan(kPi * fc / sr);
const double k = 1.0 / q;
const double denom = 1.0 + g * g + g * k; // written out rather than factored as g*(g+k)
CHECK_NEAR(c.g, g, 1e-7);
CHECK_NEAR(c.k, k, 1e-7);
CHECK_NEAR(c.a1, 1.0 / denom, 1e-7);
CHECK_NEAR(c.a2, g / denom, 1e-7);
CHECK_NEAR(c.a3, g * g / denom, 1e-7);
}
static void testTheSampleRateEntersOnlyThroughG() {
// k and the cutoff mapping are rate-free; only g moves with the rate. A reference rate
// creeping back into the module would break this.
const SvfCoeffs a = svfCoeffs(1000.0f, 2.0f, 48000.0);
const SvfCoeffs b = svfCoeffs(1000.0f, 2.0f, 96000.0);
CHECK(a.k == b.k);
CHECK(a.g != b.g);
CHECK_NEAR(b.g, std::tan(kPi * 1000.0 / 96000.0), 1e-7);
// Requesting above 0.48*sr clamps rather than diverging through tan().
const SvfCoeffs clamped = svfCoeffs(20000.0f, 1.0f, 32000.0);
CHECK_NEAR(clamped.g, std::tan(kPi * 0.48), 1e-5);
CHECK(std::isfinite(clamped.a1) && std::isfinite(clamped.a3));
// A non-positive rate yields g == 0 instead of inventing 44.1k.
CHECK(svfCoeffs(1000.0f, 1.0f, 0.0).g == 0.0f);
CHECK(svfCoeffs(1000.0f, 1.0f, -48000.0).g == 0.0f);
}
// A voice re-prepared at a non-positive rate while still ringing must not latch isSilent()
// false forever -- a future voice allocator using isSilent() as its free condition would leak
// the voice. Bypass ignores state entirely (a1=1, a2=a3=0, bypassMix reads only the input), so
// clearing it here is audibly free.
static void testNonPositiveRatePrepareClearsStaleStateAndReportsSilent() {
VoiceFilter f;
f.prepare(at(1000.0, 1.0f, kLowPass), 48000.0);
f.reset();
for (int i = 0; i < 100; ++i) {
f.process(0, static_cast<float>(std::sin(2.0 * kPi * 1000.0 * i / 48000.0)));
}
CHECK(!f.isSilent()); // genuinely ringing before the rate goes bad
f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 0.0);
CHECK(f.isSilent());
for (int i = 0; i < 480000; ++i) {
const float x = static_cast<float>(std::sin(0.1 * i));
CHECK(f.process(0, x) == x);
}
CHECK(f.isSilent());
}
// An invalid rate must pass the signal, not silence the instrument, whatever the morph asks for.
static void testNonPositiveRatePassesSignalThroughAtEveryMorph() {
for (float morph : {kHighPass, kBandPass, kLowPass}) {
VoiceFilter f;
f.prepare({0.5f, 0.5f, morph, 0.0f}, 0.0);
f.reset();
for (int i = 0; i < 64; ++i) {
const float x = static_cast<float>(std::sin(0.1 * i));
CHECK(f.process(0, x) == x);
}
}
}
// ---------------------------------------------------------------------------
// Morph — measured, under both laws
// ---------------------------------------------------------------------------
// 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 (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;
}
}
}
}
}
}
// 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<float>(i) / 16.0f;
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",
__LINE__, m, res, got, q, (got / q - 1.0) * 100.0);
++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. Grid spans the full control range (20 Hz - 20 kHz), not just three interior
// cutoffs: the residual is worse near the low-cutoff/high-rate corner (float conditioning in the
// folded x - k*v1 term as fc/sr -> 1e-4 at high Q) and is Q-dependent, so the threshold scales
// with Q rather than repeating a flat bound sized off the shallow grid. Measured worst case on
// this wider grid: 2.6e-06 (-111.7 dB) at Q=0.1, 7.0e-05 (-83.1 dB) at Q=sqrt(2), 3.2e-04
// (-69.8 dB) at Q=10, all at 192 kHz / 30 Hz — still an excellent notch, not a broadband defect.
// 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 : {20.0, 30.0, 50.0, 250.0, 1000.0, 4000.0, 16000.0, 20000.0}) {
if (fc > kRates[r] * 0.48) continue;
for (float res : {0.0f, 0.5f, 1.0f}) {
const double q = filterQFromNorm(res);
// Sized against measurement (margins 6.6x/1.55x/2.2x at Q=0.1/sqrt(2)/10 on this
// grid), not copied from the corner figure alone.
const double threshold = 1e-5 + 7e-5 * q;
const double got = measuredGain(at(fc, res, kCentre, 0.0f, MorphLaw::HighNotchLow),
kRates[r], fc, 0.25, 2.0, 0.5);
if (!(got < threshold)) {
std::printf("FAIL line %d: SEM notch at sr %.0f fc %.0f res %.1f is %.3e "
"(%.1f dB) — not a null (threshold %.3e)\n",
__LINE__, kRates[r], fc, res, got,
20.0 * std::log10(got + 1e-300), threshold);
++g_fail;
}
}
}
}
}
// 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<float>(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<float>(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
// ---------------------------------------------------------------------------
// The hard acceptance criterion, in its strongest form: at drive 0 the kernel is BIT-IDENTICAL
// to the same kernel with the limiter deleted. softLimit(x, 0) is x / sqrt(1) == x exactly, so
// this holds by algebra rather than by tolerance. Both channels and both entry points
// (process() and processFrame()) are covered, not just channel 0 through process().
struct LinearKernelRef {
SvfCoeffs c;
MorphMix mix;
float ic1 = 0.0f, ic2 = 0.0f;
float step(float x) {
const float v3 = x - ic2;
const float v1 = c.a1 * ic1 + c.a2 * v3;
const float v2 = ic2 + c.a2 * ic1 + c.a3 * v3;
ic1 = 2.0f * v1 - ic1; // no limiter at all
ic2 = 2.0f * v2 - ic2;
if (ic1 > -kFilterDenormalFloor && ic1 < kFilterDenormalFloor &&
ic2 > -kFilterDenormalFloor && ic2 < kFilterDenormalFloor) {
ic1 = 0.0f;
ic2 = 0.0f;
}
return mix.m0 * x + mix.m1 * v1 + mix.m2 * v2;
}
};
static float nextNoise(unsigned& rng) {
rng = rng * 1664525u + 1013904223u;
return static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) /
static_cast<float>(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 (MorphLaw law : kBothLaws) {
for (float morph : {kHighPass, kBandPass, kLowPass}) checkDriveZeroBitIdentity(morph, law);
}
}
// The complaint the rewrite answers: resonance must not track how hard the sample hits the
// filter unless the user asked for it. At drive 0 the response is identical over a 1000:1 level
// range; the tap this replaced moved by 14% over the same span. Runs under both laws; the centre
// is skipped under HighNotchLow because analyticMag has no notch formula to compare against there
// — level invariance at drive 0 is structural for any linear combination of the SVF's taps, so
// skipping one morph position on one law loses no real coverage.
static void testDriveZeroResponseIsLevelInvariant() {
const double sr = 48000.0, fc = 1000.0;
for (MorphLaw law : kBothLaws) {
for (float morph : {kHighPass, kBandPass, kLowPass}) {
if (morph == kBandPass && law != MorphLaw::HighBandLow) continue;
const double q = filterQFromNorm(1.0f);
const double want = analyticMag(morph, fc, fc, q, sr);
for (double amp : {0.001, 0.01, 0.1, 1.0}) {
const double got = measuredGain(at(fc, 1.0f, morph, 0.0f, law), sr, fc, amp);
if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
std::printf("FAIL line %d: %s morph %.1f amp %g gain %.6f vs analytic %.6f "
"(%.3f%%)\n",
__LINE__, lawName(law), morph, amp, got, want,
(got / want - 1.0) * 100.0);
++g_fail;
}
}
}
}
}
// Drive has to actually do something at the top of its travel, and do it monotonically — the
// brief's "extreme, not politely warm". Measured at the corner, where the resonance state is
// what the limiter sees.
static void testDriveCompressesTheResonantPeakMonotonically() {
const double sr = 48000.0, fc = 1000.0;
double prev = 1e30;
for (int i = 0; i <= 8; ++i) {
const double got = measuredGain(at(fc, 1.0f, kLowPass, i / 8.0f), sr, fc, 1.0);
CHECK(got < prev);
prev = got;
}
// Full drive against no drive: a large, unmistakable reduction of the resonant peak.
CHECK(prev < 0.5 * filterQFromNorm(1.0f));
// And the passband is left alone at every drive setting — drive colours the resonance, it
// is not a distortion box in series with the signal.
for (int i = 0; i <= 4; ++i) {
CHECK_NEAR(measuredGain(at(fc, 1.0f, kLowPass, i / 4.0f), sr, 100.0, 1.0), 1.0, 0.05);
}
}
// ---------------------------------------------------------------------------
// Sample-rate invariance
// ---------------------------------------------------------------------------
// The rate must enter only through g = tan(pi*fc/sr), so the response at a given cutoff and Q is
// the same filter at every rate. The retired feedback tap made this false: it closed the loop
// once per SAMPLE, so emphasis ran 5.02 at 48k against 8.52 at 192k. Runs under both laws; the
// centre is skipped under HighNotchLow because analyticMag has no notch formula to compare
// against there — SEM centre behavior across rates is covered by
// testHighNotchLowCentreIsATrueNullAtTheCorner instead.
static void testResponseIsRateInvariantAtEveryMorph() {
for (MorphLaw law : kBothLaws) {
for (float morph : {kHighPass, kBandPass, kLowPass}) {
if (morph == kBandPass && law != MorphLaw::HighBandLow) continue;
for (float res : {0.2f, 0.5f, 1.0f}) {
const double q = filterQFromNorm(res);
for (double fc : {250.0, 1000.0, 4000.0}) {
for (int r = 0; r < kRateCount; ++r) {
const double got = measuredGain(at(fc, res, morph, 0.0f, law), kRates[r], fc);
const double want = analyticMag(morph, fc, fc, q, kRates[r]);
if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
std::printf("FAIL line %d: %s morph %.1f res %.1f fc %.0f at %.0f Hz: "
"%.6f vs analytic %.6f (%.3f%%)\n",
__LINE__, lawName(law), morph, res, fc, kRates[r], got, want,
(got / want - 1.0) * 100.0);
++g_fail;
}
}
}
}
}
}
}
// The conditioning corner: fc/sr ~ 1e-4. Float32 Direct Form I encoded pole proximity in
// a1 -> -2, a2 -> +1 and cancelled them every sample, costing ~17 bits and putting the measured
// peak 15% LOW at 20 Hz / 192 kHz. TPT encodes the same proximity in a1's small deviation from
// 1, which float resolves; this pins that the defect is gone at every rate.
static void testLowCutoffHighRateCornerHoldsTheAnalyticPeak() {
const double q = filterQFromNorm(1.0f);
// A 2-pole low-pass peaks at W = sqrt(1 - 1/(2Q^2)), where |H| = Q / sqrt(1 - 1/(4Q^2)).
const double wPeak = std::sqrt(1.0 - 1.0 / (2.0 * q * q));
const double want = q / std::sqrt(1.0 - 1.0 / (4.0 * q * q));
CHECK_NEAR(want, 10.012516, 1e-5); // the figure the measurement pass quoted
for (int r = 0; r < kRateCount; ++r) {
const double sr = kRates[r];
const double fPeak = sr / kPi * std::atan(wPeak * std::tan(kPi * 20.0 / sr));
// Q=10 at 20 Hz rings for ~0.16 s, so the settle window has to be seconds, not samples.
const double got = measuredGain(at(20.0, 1.0f, kLowPass), sr, fPeak, 0.25, 3.0, 1.0);
if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
std::printf("FAIL line %d: 20 Hz peak at %.0f Hz is %.6f vs analytic %.6f (%.3f%%)\n",
__LINE__, sr, got, want, (got / want - 1.0) * 100.0);
++g_fail;
}
}
}
int main() {
testSvfCoefficientsMatchPinnedValues();
testTheSampleRateEntersOnlyThroughG();
testNonPositiveRatePrepareClearsStaleStateAndReportsSilent();
testNonPositiveRatePassesSignalThroughAtEveryMorph();
testMorphEndpointsMatchTheAnalyticTwoPoleTargets();
testCornerMagnitudeIsFlatAtQAcrossTheHighBandLowSweep();
testHighNotchLowCentreIsATrueNullAtTheCorner();
testHighNotchLowNullIsLocatedAtTheCutoff();
testHighNotchLowLegsHoldConstantPowerAwayFromTheNotch();
testMorphSweepHasNoDiscontinuity();
testMorphLawDoesNotDisturbTheCoefficients();
testFilterSettingsDefaultsToTheHighBandLowLaw();
testDriveZeroIsBitIdenticalToTheLinearKernel();
testDriveZeroResponseIsLevelInvariant();
testDriveCompressesTheResonantPeakMonotonically();
testResponseIsRateInvariantAtEveryMorph();
testLowCutoffHighRateCornerHoldsTheAnalyticPeak();
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;
}
+210
View File
@@ -0,0 +1,210 @@
// 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. Most interior expectations are derived from the intended law in radicals,
// sharing not even a trig call with the implementation; one check evaluates std::cos/std::sin
// directly at the same argument the implementation does, but a radical-derived check of the same
// leg sits right beside it, so no coverage rests solely on the shared call.
// 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 <cmath>
#include <cstdio>
#include <initializer_list>
#include <limits>
#include <type_traits>
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<std::underlying_type_t<MorphLaw>>(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<float>::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<float>(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<float>(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<float>(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;
}
+121
View File
@@ -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 <cmath>
#include <cstdio>
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<float>(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<double>(kFilterQMin));
const double mid = std::log(static_cast<double>(kFilterQCenter));
const double hi = std::log(static_cast<double>(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<float>(i) / 1000.0f);
CHECK(q > prev);
prev = q;
}
}
static void testQNormRoundTrips() {
for (int i = 0; i <= 20; ++i) {
const float n = static_cast<float>(i) / 20.0f;
CHECK_NEAR(filterNormFromQ(filterQFromNorm(n)), n, 1e-5);
}
CHECK_NEAR(filterNormFromQ(static_cast<float>(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<float>(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;
}
+346
View File
@@ -0,0 +1,346 @@
// Standalone tests for the running filter's NUMERICAL behaviour and state lifecycle — bounded
// output under a live parameter sweep, full-drive stability and self-oscillation, the softLimit
// shaper's own properties, 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/filter_saturate.h"
#include "../src/core/instrument/engine/filter/voice_filter.h"
#include <cfloat>
#include <cmath>
#include <cstdio>
#include <initializer_list>
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};
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;
static FilterSettings at(double fcHz, float res, float morph, float drive = 0.0f,
MorphLaw law = MorphLaw::HighBandLow) {
return {filterNormFromCutoffHz(static_cast<float>(fcHz)), res, morph, drive, law};
}
static void testFullRangeCutoffSweepAtAudioRateStaysBounded() {
unsigned rng = 0x13579bdfu;
auto noise = [&rng]() {
rng = rng * 1664525u + 1013904223u;
return static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) /
static_cast<float>(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<int>(sr * 0.25);
for (int i = 0; i < n; ++i) {
const float t = static_cast<float>(i) / static_cast<float>(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
}
}
}
}
}
}
}
// Drive is bounded by construction, not by tuning: softLimit is a contraction, so the state
// update can only ever shrink the state and the filter cannot gain energy from it. This sweeps
// the corners that would expose a tuned margin instead.
static void testFullDriveStaysBoundedAtEveryCutoffResonanceAndRate() {
unsigned rng = 0x2468aceu;
auto noise = [&rng]() {
rng = rng * 1664525u + 1013904223u;
return static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) /
static_cast<float>(1 << 22);
};
for (int r = 0; r < kRateCount; ++r) {
const double sr = kRates[r];
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;
}
}
}
}
}
}
}
}
// Full drive at full resonance with no input must still go quiet. A nonlinearity in the loop is
// exactly where a self-oscillator would hide, and softLimit's sub-unit slope is what forbids it.
static void testFullDriveDoesNotSelfOscillate() {
for (int r = 0; r < kRateCount; ++r) {
const double sr = kRates[r];
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<int>(sr * 0.01);
for (int i = 0; i < excite; ++i) {
f.process(0, static_cast<float>(std::sin(2.0 * kPi * 1000.0 * i / sr)));
}
for (int i = 0; i < static_cast<int>(sr * 0.5); ++i) f.process(0, 0.0f);
CHECK(f.isSilent());
}
}
}
}
// The softLimit shaper's own properties, independent of any running filter: bit-exact at depth 0,
// odd, a contraction over the whole excursion range, bounded by the knee, and increasing where the
// shaping actually happens.
static void testSoftLimitIsOddMonotoneBoundedAndExactAtZeroDepth() {
for (double x : {-3.0, -0.5, 0.0, 1e-9, 0.25, 7.0}) {
// Depth 0 is the identity by algebra, so correctness doesn't require special-casing it —
// voice_filter.h gates the call anyway, but as a perf optimization (see its comment).
CHECK(softLimit(static_cast<float>(x), 0.0f) == static_cast<float>(x));
}
CHECK_NEAR(softLimit(1.5f, 2.0f), -softLimit(-1.5f, 2.0f), 1e-9);
for (float depth : {0.5f, 4.0f, 64.0f}) {
// The two properties the stability argument rests on, over the whole excursion range a
// resonating state can reach. Monotonicity is NOT asserted here: far past the knee the
// curve is asymptotically flat, so the true increment between adjacent samples falls
// below float epsilon and rounding can walk it backwards by an ulp.
for (int i = -400; i <= 400; ++i) {
const float x = static_cast<float>(i) * 0.05f;
const float y = softLimit(x, depth);
CHECK(std::fabs(y) <= std::fabs(x)); // a contraction — the stability argument
CHECK(std::fabs(y) < 1.0f / depth + 1e-6f); // bounded by the knee
}
// Strictly increasing across the knee, which is where the shaping actually happens.
const float knee = 1.0f / depth;
float prev = -1e30f;
for (int i = -20; i <= 20; ++i) {
const float y = softLimit(static_cast<float>(i) * 0.1f * knee, depth);
CHECK(y > prev);
prev = y;
}
}
}
// The flush tests the ENVELOPE — both integrators — not one sample: isSilent() means "both are
// exactly zero," so both have to reach zero for that check to mean anything, and the conjunctive
// test is the cheapest guarantee of that (see voice_filter.h's flush comment). The stronger
// limit-cycle rationale belongs to the retired Direct Form I state, where the flushed variables
// were the actual filter OUTPUT rather than integrator state — it does not reproduce here.
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<int>(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<int>(sr * 0.01);
for (int i = 0; i < excite; ++i) {
f.process(0, 0.5f * static_cast<float>(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<int>(sr * 0.05);
float worstAfterSettle = 0.0f;
for (int i = 0; i < static_cast<int>(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();
testFullDriveStaysBoundedAtEveryCutoffResonanceAndRate();
testFullDriveDoesNotSelfOscillate();
testSoftLimitIsOddMonotoneBoundedAndExactAtZeroDepth();
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;
}