daniel/reasampler
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
902030bfba0dd6ce08075abd983f7effc759cb19
reasampler/tests/test_filter.cpp
T

745 lines
33 KiB
C++
Raw Blame History

// Standalone tests for the per-voice TPT/SVF filter — no VST3, no REAPER, no framework. Same
// fast assert loop as the sibling pure tests. The coefficient pins are literals so a refactor
// that changes the DSP fails loudly; they are cross-checked in-test against a derivation that
// shares no code with the implementation, and the responses against the analog 2-pole prototype
// evaluated at the bilinear-warped frequency.
#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;
// Morph positions of the three pure taps.
static constexpr float kHighPass = 0.0f;
static constexpr float kBandPass = 0.5f;
static constexpr float kLowPass = 1.0f;
// 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) {
return {filterNormFromCutoffHz(static_cast<float>(fcHz)), res, morph, drive};
}
// ---------------------------------------------------------------------------
// Control mappings (carried over — the cutoff and Q laws are unchanged)
// ---------------------------------------------------------------------------
static void testCutoffMapsThreeDecadesLogarithmically() {
CHECK_NEAR(filterCutoffHzFromNorm(0.0f), 20.0, 1e-3);
CHECK_NEAR(filterCutoffHzFromNorm(1.0f), 20000.0, 1e-2);
// Exactly three decades, so the decade midpoints land on round numbers.
CHECK_NEAR(filterCutoffHzFromNorm(1.0f / 3.0f), 200.0, 1e-3);
CHECK_NEAR(filterCutoffHzFromNorm(2.0f / 3.0f), 2000.0, 1e-2);
// Half-decade steps confirm the sweep is log, not linear.
CHECK_NEAR(filterCutoffHzFromNorm(1.0f / 6.0f), 20.0 * std::sqrt(10.0), 1e-3);
CHECK_NEAR(filterCutoffHzFromNorm(0.5f), 20.0 * std::sqrt(1000.0), 1e-2);
CHECK_NEAR(filterCutoffHzFromNorm(-1.0f), 20.0, 1e-3);
CHECK_NEAR(filterCutoffHzFromNorm(2.0f), 20000.0, 1e-2);
}
static void testCutoffNormRoundTrips() {
for (int i = 0; i <= 20; ++i) {
const float n = static_cast<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);
// 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);
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;
}
}
// ---------------------------------------------------------------------------
// 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);
}
// 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
// ---------------------------------------------------------------------------
// The endpoints are pure taps EXACTLY, not to within a rounding of cos/sin. Asserted on the
// folded mix, where "pure" is an exact statement about three floats.
static void testMorphEndpointMixesAreExactlyPureTaps() {
const float k = 1.0f / filterQFromNorm(0.5f);
const MorphMix hp = morphMix(morphWeights(kHighPass), k);
CHECK(hp.m0 == 1.0f && hp.m1 == -k && hp.m2 == -1.0f); // v0 - k*v1 - v2
const MorphMix bp = morphMix(morphWeights(kBandPass), k);
CHECK(bp.m0 == 0.0f && bp.m1 == 1.0f && bp.m2 == 0.0f); // v1
const MorphMix lp = morphMix(morphWeights(kLowPass), k);
CHECK(lp.m0 == 0.0f && lp.m1 == 0.0f && lp.m2 == 1.0f); // v2
// Out-of-range clamps to the endpoints rather than extrapolating.
CHECK(morphWeights(-1.0f).hp == 1.0f);
CHECK(morphWeights(2.0f).lp == 1.0f);
}
// HP and LP never carry weight at the same time. That is what keeps the centre a band-pass
// instead of the Oberheim SEM's notch: the two are antiphase at the corner and would cancel.
static void testMorphNeverBlendsHighAgainstLowPass() {
for (int i = 0; i <= 200; ++i) {
const MorphWeights w = morphWeights(static_cast<float>(i) / 200.0f);
CHECK(w.hp == 0.0f || w.lp == 0.0f);
CHECK(w.hp >= 0.0f && w.bp >= 0.0f && w.lp >= 0.0f);
// Equal power: the active pair sums in quadrature to unity.
CHECK_NEAR(w.hp * w.hp + w.bp * w.bp + w.lp * w.lp, 1.0, 1e-6);
}
}
static void testMorphEndpointsMatchTheAnalyticTwoPoleTargets() {
const double sr = 48000.0, fc = 1000.0;
for (float res : {0.0f, 0.5f, 1.0f}) {
const double q = filterQFromNorm(res);
for (float morph : {kHighPass, kBandPass, kLowPass}) {
for (double f : {125.0, 500.0, 1000.0, 2000.0, 8000.0}) {
const double got = measuredGain(at(fc, res, morph), sr, f);
const double want = analyticMag(morph, f, fc, q, sr);
if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
std::printf("FAIL line %d: morph %.1f res %.1f at %.0f Hz: %.6f vs analytic "
"%.6f (%.3f%%)\n",
__LINE__, morph, res, f, got, want,
(got / want - 1.0) * 100.0);
++g_fail;
}
}
}
}
}
// The reason the blend is equal-power rather than linear. At the corner the three taps are
// HP = jQ, BP = Q, LP = -jQ — adjacent taps in exact quadrature — so a cos/sin pair holds the
// corner magnitude at exactly Q the whole way across. A linear crossfade would sag to Q/sqrt(2)
// mid-leg, a 3 dB hole that would read as a defect rather than as character.
static void testCornerMagnitudeIsFlatAcrossTheWholeMorphSweep() {
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), 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;
}
}
}
}
// Continuity as a control, not just at the corner: no step between adjacent morph positions at
// any fixed frequency. A coefficient switch at the centre — the thing an enum would have forced —
// shows up here as a jump.
static void testMorphSweepHasNoDiscontinuity() {
const double sr = 48000.0, fc = 1000.0;
constexpr int kSteps = 40;
for (float res : {0.0f, 0.5f, 1.0f}) {
for (double f : {250.0, 1000.0, 4000.0}) {
double prev = -1.0;
for (int i = 0; i <= kSteps; ++i) {
const float m = static_cast<float>(i) / kSteps;
const double got = measuredGain(at(fc, res, m), sr, f);
if (prev >= 0.0) {
// Scaled by the response's own magnitude at this setting — the passband is
// unity and the corner is Q, so below Q=1 the passband is what a step has to
// be small against, not Q.
const double scale = std::fmax(1.0, filterQFromNorm(res));
// One step is 1/40 of the travel; the steepest leg moves well under a tenth
// of that scale over one step (measured worst case is 0.03).
const double jump = std::fabs(got - prev) / scale;
if (!(jump < 0.1)) {
std::printf("FAIL line %d: morph %.4f res %.1f at %.0f Hz jumps %.4f\n",
__LINE__, m, res, f, jump);
++g_fail;
}
}
prev = got;
}
}
}
}
// ---------------------------------------------------------------------------
// 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.
static void testDriveZeroIsBitIdenticalToTheLinearKernel() {
for (float morph : {kHighPass, kBandPass, kLowPass}) {
VoiceFilter f;
f.prepare(at(1000.0, 1.0f, morph, 0.0f), 48000.0);
f.reset();
const SvfCoeffs c = f.coeffs();
const MorphMix mix = f.mix();
float ic1 = 0.0f, ic2 = 0.0f;
unsigned rng = 0x13579bdfu;
for (int i = 0; i < 4096; ++i) {
rng = rng * 1664525u + 1013904223u;
const float x = static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) /
static_cast<float>(1 << 22);
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;
}
CHECK(f.process(0, x) == mix.m0 * x + mix.m1 * v1 + mix.m2 * v2);
}
}
}
// 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.
static void testDriveZeroResponseIsLevelInvariant() {
const double sr = 48000.0, fc = 1000.0;
for (float morph : {kHighPass, kBandPass, kLowPass}) {
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), sr, fc, amp);
if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
std::printf("FAIL line %d: morph %.1f amp %g gain %.6f vs analytic %.6f "
"(%.3f%%)\n",
__LINE__, morph, amp, got, want, (got / want - 1.0) * 100.0);
++g_fail;
}
}
}
}
// 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 (int ci = 0; ci <= 8; ++ci) {
for (int mi = 0; mi <= 4; ++mi) {
for (float res : {0.0f, 0.5f, 1.0f}) {
VoiceFilter f;
f.prepare({ci / 8.0f, res, mi / 4.0f, 1.0f}, sr);
f.reset();
for (int i = 0; i < 4000; ++i) {
const float y = f.process(0, noise());
if (!std::isfinite(y) || std::fabs(y) > 8.0f) {
std::printf("FAIL line %d: sr=%.0f cutoff=%.2f morph=%.2f res=%.1f "
"full drive produced %g\n",
__LINE__, sr, ci / 8.0, mi / 4.0, res, y);
++g_fail;
return;
}
}
}
}
}
}
}
// 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 (float morph : {kHighPass, kBandPass, kLowPass}) {
VoiceFilter f;
f.prepare(at(1000.0, 1.0f, morph, 1.0f), sr);
f.reset();
const int excite = static_cast<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());
}
}
}
// 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);
}
}
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 drive 0 needs no special case on the hot path.
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;
}
}
}
// ---------------------------------------------------------------------------
// 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.
static void testResponseIsRateInvariantAtEveryMorph() {
for (float morph : {kHighPass, kBandPass, kLowPass}) {
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), 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: morph %.1f res %.1f fc %.0f at %.0f Hz: %.6f "
"vs analytic %.6f (%.3f%%)\n",
__LINE__, 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;
}
}
}
// ---------------------------------------------------------------------------
// Stability, denormals, and state
// ---------------------------------------------------------------------------
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 (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}, sr);
const float y = f.process(0, noise());
CHECK(std::isfinite(y));
CHECK(std::fabs(y) < 100.0f);
if (!std::isfinite(y)) return; // stop before the log floods
}
}
}
}
}
}
// The flush tests the ENVELOPE — both integrators — not one sample. ic1 and ic2 are in
// quadrature, so a resonator swings each through zero twice a cycle; flushing on a single one
// injects a step in phase with the resonance, which the resonance amplifies, and the filter
// limit-cycles at the floor forever instead of going quiet. Re-verified for TPT rather than
// assumed to carry over from the retired Direct Form I state.
static void testStateFlushesToZeroWithoutStallingInDenormals() {
for (int r = 0; r < kRateCount; ++r) {
const double sr = kRates[r];
// The decay to the floor is a fixed WALL-CLOCK time, so the budget scales with the rate.
const int budget = static_cast<int>(sr * 0.5);
for (float morph : {kHighPass, kBandPass, kLowPass}) {
for (float drive : {0.0f, 1.0f}) {
VoiceFilter f;
f.prepare(at(1000.0, 1.0f, morph, drive), sr);
f.reset();
// Excite, then hard-cut to silence the way a released voice does.
const int excite = static_cast<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());
}
}
}
}
// A high-pass under sustained DC must settle to zero and STAY there. Sampling only the final
// value is not enough: a resonator swings through zero twice a cycle, so a single late sample
// can land near zero while the envelope still rings well above it. This regressed a click train
// on the retired topology, where flushing the FIR history discarded the pinned DC and the next
// sample recomputed a full-amplitude step. TPT has no FIR history to discard, so the hazard is
// structural rather than a tuning — but the assertion is cheap and pins the outcome.
static void testHighPassSustainedDCDoesNotReRing() {
for (int r = 0; r < kRateCount; ++r) {
const double sr = kRates[r];
for (float drive : {0.0f, 1.0f}) {
VoiceFilter f;
f.prepare(at(1000.0, 1.0f, kHighPass, drive), sr);
f.reset();
const int settle = static_cast<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);
}
static void testLowpassStepSettlesToUnityAndHighpassRejectsDC() {
const double sr = 48000.0;
VoiceFilter f;
f.prepare(at(1000.0, 0.0f, kLowPass), sr);
f.reset();
float y = 0.0f;
for (int i = 0; i < 48000; ++i) y = f.process(0, 1.0f);
CHECK_NEAR(y, 1.0, 1e-3); // DC passes a lowpass at unity
VoiceFilter hp;
hp.prepare(at(1000.0, 0.0f, kHighPass), sr);
hp.reset();
float worstAfterSettle = 0.0f;
for (int i = 0; i < 48000; ++i) {
y = hp.process(0, 1.0f);
if (i >= 200) worstAfterSettle = std::fmax(worstAfterSettle, std::fabs(y));
}
CHECK(worstAfterSettle < 1e-3f);
}
static void testResetClearsStateButPrepareKeepsIt() {
VoiceFilter f;
f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 48000.0);
f.process(0, 1.0f);
CHECK(!f.isSilent());
// A live parameter move must not zero the state — that is what would click.
f.prepare({0.6f, 0.5f, kLowPass, 0.0f}, 48000.0);
CHECK(!f.isSilent());
f.prepare({0.6f, 0.5f, kBandPass, 1.0f}, 48000.0);
CHECK(!f.isSilent());
f.reset();
CHECK(f.isSilent());
}
static void testChannelStateIsIndependent() {
VoiceFilter f;
f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 48000.0);
f.reset();
f.process(0, 1.0f);
CHECK(f.state(0).ic2 != 0.0f);
CHECK(f.state(1).ic2 == 0.0f);
float frame[2] = {1.0f, -1.0f};
f.processFrame(frame, 2);
CHECK(f.state(1).ic2 < 0.0f);
CHECK(frame[0] != frame[1]);
}
int main() {
testCutoffMapsThreeDecadesLogarithmically();
testCutoffNormRoundTrips();
testQSpansPointOneToTenWithRootTwoAtCenter();
testQNormRoundTrips();
testDriveDepthIsZeroAtRestAndRisesMonotonically();
testSvfCoefficientsMatchPinnedValues();
testTheSampleRateEntersOnlyThroughG();
testNonPositiveRatePassesSignalThroughAtEveryMorph();
testMorphEndpointMixesAreExactlyPureTaps();
testMorphNeverBlendsHighAgainstLowPass();
testMorphEndpointsMatchTheAnalyticTwoPoleTargets();
testCornerMagnitudeIsFlatAcrossTheWholeMorphSweep();
testMorphSweepHasNoDiscontinuity();
testDriveZeroIsBitIdenticalToTheLinearKernel();
testDriveZeroResponseIsLevelInvariant();
testFullDriveStaysBoundedAtEveryCutoffResonanceAndRate();
testFullDriveDoesNotSelfOscillate();
testDriveCompressesTheResonantPeakMonotonically();
testSoftLimitIsOddMonotoneBoundedAndExactAtZeroDepth();
testResponseIsRateInvariantAtEveryMorph();
testLowCutoffHighRateCornerHoldsTheAnalyticPeak();
testFullRangeCutoffSweepAtAudioRateStaysBounded();
testStateFlushesToZeroWithoutStallingInDenormals();
testHighPassSustainedDCDoesNotReRing();
testImpulseResponseMatchesTheKernel();
testLowpassStepSettlesToUnityAndHighpassRejectsDC();
testResetClearsStateButPrepareKeepsIt();
testChannelStateIsIndependent();
if (g_fail == 0) std::printf("filter_tests: all passed\n");
else std::printf("filter_tests: %d FAILED\n", g_fail);
return g_fail == 0 ? 0 : 1;
}
Reference in New Issue View Git Blame Copy Permalink