Break up the de-emphasis dependency chain in the float decoder

[fippo]

The de-emphasis filter is a first-order IIR that runs per channel, each
with its own filter memory, so the scalar loop is latency-bound: every
output waits for the multiply-add that produced the previous memory
value.

Expand the recurrence over blocks of four samples instead. The serial
chain is one multiply-add per sample, m = coef0 * (x[j] + m). After
the expansion only the carry coef0^4 * m is left on it, one per four
samples. The rest becomes independent work the compiler can pipeline and
vectorize. Remaining samples are handled by the existing scalar loop.

The transform relies on a*(x+m) being interchangeable with a*x + a*m
which does not hold in fixed point due to saturation so this is
float-only. The result is not bit-exact but the pole is contractive:

• over 120 s of speech 30 of 5.76M samples differ
• over 64 s of stereo music 612 of 6.19M
  in both cases by one LSB at most.

opus_demo -d on x86-64 with gcc -O3, split by the mode in the TOC
byte:
120 s of speech at 48 kHz mono:

mode	deemphasis	total
SILK-only, WB	n/a, the CELT decoder does not run in this mode
hybrid, FB	8.6% -> 3.0%	-6.2%
CELT-only, FB	8.1% -> 3.1%	-4.8%

64 s of music at 48 kHz stereo:

mode	deemphasis	total
CELT-only, FB	5.3% -> 3.6%	-0.5%

Bit-identical in fixed-point (obviously)

[jmvalin]

I just landed commit 8207ebb that does something similar with Neon assembly. Not sure what arch you're targeting (if any in particular), but I suspect the asm/intrinsics route may just be the way to go instead. I suspect it shouldn't be too hard to do for SSE.

[fippo]

Rebased, results still stand and I had Claude give the SSE a try because why not.

Hybrid stereo mode (which happens in the wild world of WebRTC I am aiming at) doesn't change at all vs this change alone. CELT-only gets an additional -0.5% (mono) to 0.9% (stereo).
For CELT-only stereo that is half of what this yeilds (-1.8%) so might be worth it for a 30 line function.

[fippo]

specializing C for 1 or 2 would let the compiler unroll this but yielded <1% at the expense of duplicate code so I didn't do it.

[fippo]

drive-by: "common" here is celt-only stereo, notably not hybrid. opus_decoder.c:370 sets celt_accum = (mode != MODE_CELT_ONLY).
Dropping !accum would have been ~6% but becomes vanishingly small with this patch (which the neon assembly might have gone after?)

[jmvalin]

I'd say my main issue here is just how it makes the deemphasis code even harder to follow than the mess it already was. Any way you can make this more self-contained? Also, just curious, how does your code compare to the Neon version that @cyanreg just landed?

[fippo]

Switching to math writing, the old form is basically

  for n = 0 .. N-1:
      t[n] = x[n] + m
      m    = c * t[n]
      y[n] = t[n]

which in this change becomes

  for j = 0; j + 4 <= N; j += 4:
      x0, x1, x2, x3 = x[j..j+3]

      y[j+0] = x0                                +      m
      y[j+1] = x1 + c*x0                         + c  * m
      y[j+2] = x2 + c*x1 + c^2*x0                + c^2* m
      y[j+3] = x3 + c*x2 + c^2*x1 + c^3*x0       + c^3* m

      m = c*x3 + c^2*x2 + c^3*x1 + c^4*x0        + c^4* m    # = c * y[j+3]
  # tail omittted

The neon code does the same but writes it slightly differently with m as a vector.

  (c, c^2, c^3, c^4) precomputed
  mv[0..3] = (m, c*m, c^2*m, c^3*m)

  for j = 0; j + 4 <= N; j += 4:
      x0, x1, x2, x3 = x[j..j+3]

      y[j+0] = x0                                + mv[0]
      y[j+1] = x1 + c*x0                         + mv[1]
      y[j+2] = x2 + c*x1 + c^2*x0                + mv[2]
      y[j+3] = x3 + c*x2 + c^2*x1 + c^3*x0       + mv[3]

      mv[0..3] = y[j+3] * (c, c^2, c^3, c^4)

  m = mv[0]
  # tail omittted

Updated in the vector carry commit. I like it slightly less but that is just because my brain is already wired for the other way of writing it.

Now LLMs are great because they don't complain when asked to do a line-by-line annotation of that:

  (c, c^2, c^3, c^4) precomputed                    # C: L316-319   NEON: L144-156

  mv[0..3] = (m, c*m, c^2*m, c^3*m)                 # C: L290-296   NEON: L162-164

  for j = 0; j + 4 <= N; j += 4:                    # C: L352       NEON: L169,L183-184 (x8)
      x0, x1, x2, x3 = x[j..j+3] (+ VERY_SMALL)     # C: L321-324   NEON: L169 ld1 (no VERY_SMALL)

      y[j+0] = x0                          + mv[0]  # C: L325,L329  NEON: L170-176 (three
      y[j+1] = x1 + c*x0                   + mv[1]  # C: L326,L330  NEON:  fmla + one fadd;
      y[j+2] = x2 + c*x1 + c^2*x0          + mv[2]  # C: L327,L331  NEON:  all four lanes
      y[j+3] = x3 + c*x2 + c^2*x1 + c^3*x0 + mv[3]  # C: L328,L332  NEON:  at once)

      mv[0..3] = y[j+3] * (c, c^2, c^3, c^4)        # C: L333-336   NEON: L179 (pair: L178)

      store SIG2RES(y[j..j+3])                      # C: L357-358   NEON: L180-182

  m = mv[0]                                         # C: L360       NEON: L196 (s0 = lane 0)

  for ; j < N; j += 1:                              # C: L363       NEON: L187-189
      t[j] = x[j] + m                               # C: L365       NEON: L190
      m    = c * t[j]                               # C: L366       NEON: L191
      y[j] = t[j]  (SIG2RES)                        # C: L367       NEON: L192-193

which shows one can align the two even more (in a separate commit; lines above refer to the second one)
(will update the commit message with the math when squashing)

Any way you can make this more self-contained?

The deemph_loop reorganises it a bit.