audioecho.FeedbackDelay: a cross-fed stereo tail reaches exact zero (#170)

src/shared/audiodsp_feedback_delay.c

@@ -342,21 +342,55 @@ static int16_t to_s16(float value) {
 // was enough to hand that line back (audiodsp#157), so near the floor a state
 // that has stopped moving is set onto its input (E11).
 //
+// The cross-feed is convex only in exact arithmetic. In float32
+// `own * direct + other * crossed` can land an ulp above both lanes: with
+// both lanes of a stereo line at k, one or two float32 steps under
+// f = 1 - 0.5 / k, the test saw a gap a hair over 0.5, rounded, and wrote k
+// back on every lap (audiodsp#170, E14). So a gap within EDGE of 0.5 is not
+// trusted: there the write is worked out, and rounding stands only if what
+// comes back is smaller than the larger lane. That is the property the whole
+// argument needs, checked directly, so it holds however the sum was rounded
+// -- each product on its own or fused into a multiply-add, in either order --
+// and it never touches a write that was already smaller than both lanes,
+// which is every write in mono, at cross-feed 0 or 1, and in stereo wherever
+// the lanes differ. How far the sum can overshoot is bounded: a few float32
+// steps of a value no larger than 51 (the test only bites where |sent| <=
+// 0.5 / (1 - 0.99)), under 1.1e-5, so EDGE is 2^-14 = 6.1e-5 and a gap past
+// it goes straight to rounding, one compare, as before. (A margin on the
+// test alone, the same constant with no check behind it, also ends every
+// tail, but it moves the edge for every sample that crosses it, and loud
+// passages cross it at every zero crossing: it moved 90 of 286 audible
+// renders by a 1 LSB flip carried round the loop.)
+//
 // Why not truncate always, as audiodsp_tank.c's tank_quantize() does: that
 // takes half an LSB on average from every pass at every level, and a short
 // loop makes hundreds of passes a second. On a 3 ms flanger at the clamp it
 // ended the ring half a second early, 5 dB down at 1.5 s. Stepping only where
 // rounding would stall keeps every repeat above the stall region
 // byte-identical and changes only the part that used to circulate for ever.
-static float recirculated(float sent, float feedback) {
+#define AUDIODSP_FEEDBACK_DELAY_STALL_EDGE (1.0f / 16384.0f)
+
+static float recirculated(float sent, float own, float other,
+    float feedback) {
     const float back = feedback * sent;
     const float in_size = sent < 0.0f ? -sent : sent;
     const float back_size = back < 0.0f ? -back : back;
-    if (in_size - back_size > 0.5f) {
+    const float gap = in_size - back_size;
+    if (gap > 0.5f + AUDIODSP_FEEDBACK_DELAY_STALL_EDGE) {
         return back;
     }
-    // Here |back| <= |sent| <= 0.5 / (1 - feedback), at most 50 under the
-    // 0.99 clamp, so the cast is defined.
+    if (gap > 0.5f) {
+        // Rounded, `back` would write this (to_s16 of it alone); keep it only
+        // if that is smaller than the larger lane (audiodsp#170).
+        const float own_size = own < 0.0f ? -own : own;
+        const float other_size = other < 0.0f ? -other : other;
+        const float larger = own_size > other_size ? own_size : other_size;
+        if ((float)(int32_t)(back_size + 0.5f) < larger) {
+            return back;
+        }
+    }
+    // Here |back| <= |sent| <= (0.5 + EDGE) / (1 - feedback), about 50
+    // under the 0.99 clamp, so the cast is defined.
     return (float)(int32_t)back;
 }
 
@@ -586,7 +620,8 @@ void audiodsp_feedback_delay_process_s16(
             const float fed = source * config->feed_own[channel] +
                 other_source * config->feed_other[channel];
             state->line[(size_t)channel * length + state->write_frame] =
-                to_s16(fed + recirculated(sent, config->feedback));
+                to_s16(fed + recirculated(sent, loop[channel], other,
+                    config->feedback));
             // The dry path is the channel's own signal, never the panned
             // one: `input_pan` steers what goes round the loop, not what the
             // listener hears straight through.

tests/parity/feedback_delay_state_probe.py

@@ -1,6 +1,6 @@
 """Deterministic FeedbackDelay PCM for the node's own state: the damping
-low-pass near the floor, the loop filters taken out and put back, and a wow
-depth moved while playing.
+low-pass near the floor, the loop filters taken out and put back, a wow
+depth moved while playing, and a cross-fed stereo line near the floor.
 
     feedback_delay_state_probe.py audioecho
 
@@ -110,3 +110,27 @@ def emit(tag, node, blocks):
 emit("depth-b", node, 3)
 node.set(wow_depth_ms=0.0)
 emit("depth-c", node, 3)
+
+# audiodsp#170. A stereo line cross-fed at s / 127, at two of the feedbacks
+# where the float32 cross-feed sum landed an ulp above a landed line of k and
+# the stall test handed k back for ever (k = 9 at 1/127, k = 50 at 39/127).
+# A DC of 2k + 2 for one block, then silence, at 8 kHz with a 1 ms line, so a
+# lap is eight frames and the floor comes in 60 blocks without a buffer the
+# default MicroPython heap cannot hold; the approach in one number, then the
+# last four blocks.
+for feedback, cross_feed, k in ((0.9444443583, 1.0 / 127.0, 9),
+                                (0.9899999499, 39.0 / 127.0, 50)):
+    node = echo.FeedbackDelay(sample_rate=8000, max_delay_ms=3.0,
+                              delay_ms=1.0, feedback=feedback, mix=2.0,
+                              cross_feed=cross_feed)
+    values = array("h", [2 * k + 2] * (256 * 2))
+    values.extend(array("h", bytes(256 * 63 * 2 * 2)))
+    node.play(audiocore.RawSample(values, sample_rate=8000,
+                                  channel_count=2))
+    tag = "crossfed-%d" % k
+    whole = 2166136261
+    for _block in range(60):
+        for byte in bytes(audiocore.get_buffer(node)[1]):
+            whole = ((whole ^ byte) * 16777619) & 0xffffffff
+    print("fbds", tag, "approach", whole)
+    emit(tag, node, 4)

tests/parity/verify_dsp.py

@@ -99,7 +99,8 @@
     ("feedback_delay_probe.py", "audioecho", {}, None),
     ("feedback_delay_options_probe.py", "audioecho", {}, None),
     # The node's own state: the damping floor (#157), then the filters out
-    # and back (#158, #159) and a moved wow depth (#160). Its own file, so
+    # and back (#158, #159), a moved wow depth (#160) and the cross-fed
+    # floor (#170). Its own file, so
     # the two above keep rendering what they rendered.
     ("feedback_delay_state_probe.py", "audioecho", {}, None),
     ("ladder_probe.py", "audioladder", {}, None),

tests/test_cpython_audioecho.py

@@ -26,6 +26,7 @@
 | E11 | E10 with the damping low-pass in, at the feedbacks where 0.5 / (1 - feedback) is whole | exact, 0 LSB |
 | E12 | A loop filter set to 0 and back plays nothing stale: silence stays silence, a steady line stays put, a high-pass comes back as a fresh one | exact, 0 LSB |
 | E13 | A wow depth moved on a live node glides in over 20 ms, then renders as a node built with it | no step past 1.5 x the tone's own; exact after the ramp |
+| E14 | E10 in stereo with a cross-feed strictly between 0 and 1, at the feedbacks where the cross-fed sum lands an ulp above its lanes | exact, 0 LSB |
 
 **E10 is audiodsp#153.** The feedback write rounded to nearest, so a repeat x
 came back as round(feedback * x) and every |x| <= 0.5 / (1 - feedback) was its
@@ -53,6 +54,7 @@
 """
 
 from array import array
+import struct
 import unittest
 
 import audiocore
@@ -719,5 +721,115 @@ def test_after_the_ramp_it_is_the_node_built_with_that_depth(self):
                                  built[160 * CHANNELS:])
 
 
+def _f32(value):
+    return struct.unpack("f", struct.pack("f", value))[0]
+
+
+def _f32_step(value, steps):
+    bits = struct.unpack("I", struct.pack("f", value))[0]
+    return struct.unpack("f", struct.pack("I", bits + steps))[0]
+
+
+def crossfeed_stall_cells():
+    """The (s, k, feedback) cells where audiodsp main's float32 arithmetic
+    handed a landed stereo line of k back: cross_feed s / 127, both lanes at
+    k, feedback within three float32 steps of 1 - 0.5 / k. Worked out here in
+    the node's own float32 order (each product rounded, then the sum; the
+    stall test at 0.5 as it was), so the cells do not depend on the build
+    under test."""
+    def write(sent, feedback):
+        back = _f32(feedback * sent)
+        if _f32(abs(sent) - abs(back)) <= 0.5:
+            back = float(int(back))
+        return int(_f32(back + 0.5))
+
+    cells = []
+    for s in range(1, 127):
+        crossed = _f32(s / 127.0)
+        direct = _f32(1.0 - crossed)
+        for k in range(1, 51):
+            sent = _f32(_f32(k * direct) + _f32(k * crossed))
+            if not sent > k:
+                continue
+            centre = _f32(1.0 - 0.5 / k)
+            for steps in range(-3, 4):
+                feedback = _f32_step(centre, steps)
+                if write(sent, feedback) >= k > write(float(k), feedback):
+                    cells.append((s, k, feedback))
+    return cells
+
+
+class CrossFedTailReachesZeroTest(unittest.TestCase):
+    """E14 - audiodsp#170. E10 in stereo with the cross-feed strictly
+    between 0 and 1.
+
+    With both lanes of the line at the same whole k, each lane sends
+    `k * (1 - cross_feed) + k * cross_feed` round the loop. That is k, but
+    the two float32 products could sum to an ulp above it; one or two float32
+    steps under f = 1 - 0.5 / k the stall test then saw a gap a hair over 0.5,
+    rounded, and wrote k back on every lap: 117 cells on the 7-bit cross_feed
+    grid, 9 to 50 LSB held for ever, while the same feedback ended at
+    cross_feed 0 and in mono. Within 2^-14 of that edge the node now works
+    out the write and keeps the rounding only if it is smaller than the
+    larger lane, which holds however the sum is rounded and leaves every
+    other write as it was.
+    """
+
+    def _held(self, rate, delay_ms, feedback, cross_feed, level,
+              channels=2, silent_blocks=40):
+        lap = int(rate * delay_ms / 1000.0 + 0.5)
+        lead = (4 * lap + 255) // 256 * 256
+        data = array("h", [level] * (lead * channels))
+        data.extend(array("h", bytes(2 * channels * 256 * silent_blocks)))
+        node = audioecho.FeedbackDelay(
+            sample_rate=rate, channel_count=channels,
+            max_delay_ms=delay_ms + 2.0, delay_ms=delay_ms,
+            feedback=feedback, mix=2.0, cross_feed=cross_feed)
+        node.play(audiocore.RawSample(data, sample_rate=rate,
+                                      channel_count=channels))
+        rendered = words(node, lead // 256 + silent_blocks)
+        self.assertGreater(max(abs(v) for v in rendered[:lead * channels]),
+                           0)
+        return max(abs(v) for v in rendered[-lap * channels:])
+
+    def test_the_search_finds_the_cells_it_was_written_for(self):
+        """The cells below are the 117 AnalogDelay's audit found."""
+        cells = crossfeed_stall_cells()
+        self.assertEqual(len(cells), 117)
+        self.assertEqual(len(set(k for _s, k, _f in cells)), 30)
+
+    def test_a_cross_fed_stereo_tail_reaches_exact_zero(self):
+        """E14: every cell, both signs, a 1 ms line at 8 kHz so a lap is
+        eight frames: the last lap is silent."""
+        for s, k, feedback in crossfeed_stall_cells():
+            for level in (2 * k + 2, -(2 * k + 2)):
+                with self.subTest(s=s, k=k, feedback=feedback, level=level):
+                    self.assertEqual(self._held(8000, 1.0, feedback,
+                                                s / 127.0, level), 0)
+
+    def test_the_reported_cells_at_48k(self):
+        """E14 as it was found: 48 kHz, 20 ms, and the typed feedback
+        0.9899999 with cross_feed 39 / 127 that held 50 LSB. Four laps of
+        DC build the line to about 4 x the level, and at 0.99 that takes
+        some 300 laps to come down; 1 600 blocks is 426 laps."""
+        for feedback, s, k in ((0.9444443583, 1, 9), (0.9545453787, 3, 11),
+                               (0.9899999499, 2, 50), (0.9899999, 39, 50),
+                               (0.9666666, 21, 15)):
+            with self.subTest(feedback=feedback, s=s):
+                self.assertEqual(self._held(48000, 20.0, feedback, s / 127.0,
+                                            2 * k + 2, silent_blocks=1600), 0)
+
+    def test_the_controls_end_as_they_did(self):
+        """Cross-feed 0 and 1 and mono ended before the fix and still do."""
+        for s, k, feedback in crossfeed_stall_cells()[::9]:
+            for cross_feed, channels in ((0.0, 2), (1.0, 2),
+                                         (s / 127.0, 1)):
+                with self.subTest(s=s, k=k, cross_feed=cross_feed,
+                                  channels=channels):
+                    self.assertEqual(self._held(8000, 1.0, feedback,
+                                                cross_feed, 2 * k + 2,
+                                                channels=channels), 0)
+
+
 if __name__ == "__main__":
     unittest.main()