diff --git a/src/math_channel.cpp b/src/math_channel.cpp index b01c7e2..42a8a77 100644 --- a/src/math_channel.cpp +++ b/src/math_channel.cpp @@ -59,6 +59,28 @@ void biquad(MathOp op, double fc, double q, double dt, double& b0, double& b1, d /// One biquad step, transposed direct form II -- the form whose state words are the two things /// that actually need carrying between samples. +/// One pole's coefficient for a corner of `fc` Hz at `dt` seconds a sample: `a = 1 - exp(-2*pi*fc*dt)`, +/// clamped just under Nyquist so `a` cannot exceed 1. `fc == 0` returns 1, which is a pass-through -- so +/// "no pole" costs the step no branch. +double pole_a(double fc, double dt) { + if (!(fc > 0.0) || !(dt > 0.0)) return 1.0; + const double nyquist = 0.5 / dt; + const double use = fc < nyquist * 0.99 ? fc : nyquist * 0.99; + return 1.0 - std::exp(-2.0 * kPi * use * dt); +} + +/// The median window's length: odd, and between 3 and 9. +/// +/// Odd because an even window's median is the average of two samples, which is a value the signal never +/// took -- and the whole point of a median here is that its output is a sample that was really there. +/// Rounded *up*, so a request for 4 gets 5 rather than 3. The ceiling of 9 is what `math_step`'s fixed +/// `double sorted[9]` is sized for. +int median_n(int want) { + int n = want < 3 ? 3 : (want > 9 ? 9 : want); + if ((n & 1) == 0) ++n; + return n; +} + double biquad_step(double x, double b0, double b1, double b2, double a1, double a2, double& z1, double& z2) { // A reserved hole is NaN in the source, and a NaN integrated into the state words is @@ -90,11 +112,12 @@ const char* math_op_name(MathOp op) { case MathOp::Sub: return "difference"; case MathOp::Mul: return "product"; case MathOp::Div: return "ratio"; + case MathOp::DeNoise: return "de-noise"; } return "?"; } -int math_op_count() { return 12; } +int math_op_count() { return 13; } MathOp math_op_at(int index) { switch (index) { @@ -109,7 +132,8 @@ MathOp math_op_at(int index) { case 8: return MathOp::Sum; case 9: return MathOp::Sub; case 10: return MathOp::Mul; - default: return MathOp::Div; + case 11: return MathOp::Div; + default: return MathOp::DeNoise; } } @@ -136,6 +160,8 @@ bool math_op_uses_filter(MathOp op) { return op == MathOp::LowPass || op == MathOp::HighPass || op == MathOp::BandPass; } bool math_op_uses_clamp(MathOp op) { return op == MathOp::Clamp; } +bool math_op_uses_fc(MathOp op) { return math_op_uses_filter(op) || op == MathOp::DeNoise; } +bool math_op_uses_median(MathOp op) { return op == MathOp::DeNoise; } void math_prepare(const MathChannel& channel, MathRuntime& rt, double dt, const std::function& index_of) { @@ -152,6 +178,17 @@ void math_prepare(const MathChannel& channel, MathRuntime& rt, double dt, } rt.dt = dt; rt.ready = true; + // The median windows, pooled: one pass to size the pool, then the node loop hands out the slices. A + // window starts from rest after a re-prepare -- the same choice the biquads' `z2` makes -- while the + // pole's `z1` is carried, so a decimation change does not make the trace jump. + std::size_t med_total = 0; + for (std::size_t i = 0; i < n; ++i) { + if (channel.nodes[i].op == MathOp::DeNoise) { + med_total += static_cast(median_n(channel.nodes[i].med_n)); + } + } + rt.median.assign(med_total, 0.0); + std::size_t med_next = 0; for (std::size_t i = 0; i < n; ++i) { const MathNode& node = channel.nodes[i]; MathNodeRt& out = rt.node[i]; @@ -159,6 +196,15 @@ void math_prepare(const MathChannel& channel, MathRuntime& rt, double dt, out.src = index_of ? index_of(node.source) : -1; continue; } + if (node.op == MathOp::DeNoise) { + out.med_n = median_n(node.med_n); + out.med_at = med_next; + out.med_pos = 0; + out.med_have = 0; + med_next += static_cast(out.med_n); + out.b0 = pole_a(node.fc, dt); + continue; + } if (!math_op_uses_filter(node.op)) continue; biquad(node.op, node.fc, node.q, dt, out.b0, out.b1, out.b2, out.a1, out.a2); } @@ -221,6 +267,46 @@ double math_step(const MathChannel& channel, MathRuntime& rt, const float* value x = std::fabs(d) > 0.0 ? at(node.in_a) / d : 0.0; break; } + case MathOp::DeNoise: { + const double in = at(node.in_a); + // A hole stays a hole, and the window is left alone: folding a NaN into the median would smear one + // gap across `med_n` samples, and the whole trace after it in the pole. `med_n <= 0` is a node the + // runtime has not prepared -- the window pointer below would be a read past the pool, so it is + // checked rather than assumed. + if (!std::isfinite(in) || self.med_n <= 0) { + x = kNaN; + break; + } + const std::size_t len = static_cast(self.med_n); + double* win = rt.median.data() + self.med_at; + win[self.med_pos] = in; + self.med_pos = self.med_pos + 1 == len ? 0 : self.med_pos + 1; + if (self.med_have < len) ++self.med_have; + // An insertion sort of at most nine doubles into a local. The window is tiny and fixed, so this + // beats anything that allocates, and it leaves the ring itself in arrival order -- which is what + // lets the oldest sample be overwritten in place. Before the window has filled, the median is taken + // over what there is: a shorter window, deliberately, rather than NaN. A derived trace that reads + // "no data" for its first few samples after every rebuild looks like a broken graph. + double sorted[9]; + const std::size_t m = self.med_have; + for (std::size_t j = 0; j < m; ++j) { + const double v = win[j]; + std::size_t p = j; + while (p > 0 && sorted[p - 1] > v) { + sorted[p] = sorted[p - 1]; + --p; + } + sorted[p] = v; + } + const double med = sorted[m / 2]; + // Median *then* pole, and the order is the design. The median kills an impulse outright; the pole + // then takes the broadband hash it leaves behind. The other way round feeds the pole the impulse, + // which smears it over the pole's whole settling time and leaves the median a wide bump it cannot + // reject. + self.z1 += self.b0 * (med - self.z1); + x = self.z1; + break; + } } self.value = x; } diff --git a/src/math_channel.hpp b/src/math_channel.hpp index 4c9a569..2690281 100644 --- a/src/math_channel.hpp +++ b/src/math_channel.hpp @@ -25,6 +25,10 @@ enum class MathOp { Sub, ///< a - b Mul, ///< a * b Div, ///< a / b, zero when b is 0 (a scope trace must not go to inf/NAN) + /// A running median of `med_n` samples, then a one-pole at `fc` Hz. For impulsive interference -- an EM + /// peak a few samples wide -- which a linear filter can only smear over its own settling time. `fc = 0` + /// is the median alone. + DeNoise, }; const char* math_op_name(MathOp op); @@ -41,6 +45,11 @@ int math_op_arity(MathOp op); bool math_op_uses_gain(MathOp op); bool math_op_uses_filter(MathOp op); bool math_op_uses_clamp(MathOp op); +/// Whether an op reads `fc`. Wider than `math_op_uses_filter()`, which keeps its present meaning of *a +/// biquad*: `DeNoise` reads `fc` for a one-pole and has no `q`. +bool math_op_uses_fc(MathOp op); +/// Whether an op reads `med_n`. +bool math_op_uses_median(MathOp op); /// One node of a math channel's graph. Node indices are positions in `MathChannel::nodes`, and a /// node may only read ones before it: the evaluator walks the list in order, which is what makes a @@ -56,6 +65,10 @@ struct MathNode { double q = 0.707; ///< BP quality factor. double lo = -1.0; ///< `Clamp` bounds. double hi = 1.0; + /// `DeNoise`: the median window, in samples. Odd and between 3 and 9 -- `median_n()` in the + /// implementation is what enforces that, so a session file holding anything else still loads. An impulse + /// up to `med_n / 2` samples wide is rejected outright. + int med_n = 5; double x = 0.0; ///< the editor's box position, logical pixels. double y = 0.0; @@ -67,7 +80,7 @@ struct MathNode { bool operator==(const MathNode& other) const { return op == other.op && source == other.source && in_a == other.in_a && in_b == other.in_b && gain == other.gain && offset == other.offset && fc == other.fc && q == other.q && - lo == other.lo && hi == other.hi; + lo == other.lo && hi == other.hi && med_n == other.med_n; } }; @@ -110,10 +123,25 @@ struct MathNodeRt { /// per source node, per sample, behind a type-erased call. Names cannot change without a /// rebuild, and every rebuild re-prepares. int src = -1; + /// `DeNoise`: where this node's median window starts in `MathRuntime::median`, how long it is, the ring's + /// write position, and how many samples it holds. The window is *pooled* rather than held here because a + /// node's state belongs on one cache line, and nine doubles would push this struct across two for every + /// node in the graph, filter or not. + /// + /// The one-pole reuses `b0` for its coefficient and `z1` for its value, which is not a shortcut: + /// `math_prepare()` deliberately carries `z1` across a re-prepare, so the pole keeps its history through + /// a decimation change exactly as the biquads keep theirs. + std::size_t med_at = 0; + int med_n = 0; + std::size_t med_pos = 0; + std::size_t med_have = 0; }; struct MathRuntime { std::vector node; + /// Every `DeNoise` node's median window, end to end: one allocation for the channel, sized by + /// `math_prepare()`, a node's slice starting at its own `med_at`. + std::vector median; /// The rate the coefficients were built for; a change rebuilds them. double dt = 0.0; bool ready = false;