Skip to content

silk: Arm NEON for the float SILK analysis path (inner_product, warped autocorrelation, energy) - #481

Open
czoli1976 wants to merge 4 commits into
xiph:mainfrom
czoli1976:arm-silk-float-neon
Open

czoli1976 wants to merge 4 commits into
xiph:mainfrom
czoli1976:arm-silk-float-neon

Conversation

@czoli1976

@czoli1976 czoli1976 commented Jun 16, 2026 •

Copy link
Copy Markdown

Summary

The float SILK analysis path (silk/float/) had no Arm SIMD at all (only x86 has an AVX2 inner_product_FLP). This adds the first silk/float/arm/ NEON tier — three kernels — wired into the existing RTCD dispatch and all three build systems (autotools / CMake / Meson).

Numerical model & tolerance (addressing review feedback)

Kernel Bit‑exact vs C? SIMD error Speed (M1 Pro) Speed (M4, author)
warped_autocorrelation_FLP Yes (maxdiff = 0) 0 1.21–1.42× ~1.34× geomean
inner_product_FLP No (1 ULP ≈ 1e‑16) ≤9.7e‑16 1.2–2.1× ~1.59× geomean
energy_FLP No (1 ULP ≈ 1e‑16) ≤3.4e‑14 1.2–1.8× ~1.26×
  • warped_autocorrelation_FLP is bit‑exact. The all‑pass cascade is loop‑carried, so it stays scalar in double (identical state[] to the C reference) and only the independent per‑lag SAXPY is vectorised — each C[k] gets one fused input[n]*state[k] per sample, in the same n‑order as the C reference's fmadd chain. → byte‑identical bitstream.

  • inner_product_FLP / energy_FLP are NOT bit‑exact (≤1 ULP ≈ 1e‑16). The SIMD speed comes entirely from parallelising the f64 reduction: (p0+p2)+(p1+p3) instead of the C reference's serial ((p0+p1)+p2)+p3. f32×f32 widens to f64 exactly (48 sig bits fit in 53), so the products are exact — the only rounding is in the accumulation order, which is not associative. A bit‑exact reorder matching C's serial chain was prototyped and runs at 0.32–0.71× (slower than C) — serial latency is the hard floor. Kahan/Neumaier compensation was also tested: 0.14–0.27× (5–7× slower than C) with zero error reduction (1 ULP is already the floor for these inputs). Opus's existing 1e‑6/1e‑5 float tolerances absorb this; end‑to‑end encode is byte‑identical.

Regression found & fixed: order‑24 warped autocorrelation

An earlier revision ran the accumulation as a separate vector pass re‑reading the full state[] from the stack (25 doubles at order 24) plus per‑sample loop branches. That traversal dominated at order 24 and produced a regression to 0.84× vs C. This commit rewrites the kernel to:

  1. Fold the SAXPY inline inside the all‑pass loop (consume tmp1/tmp2 straight from registers — zero state[] re‑reads).
  2. Unroll two sections (4 lags) per iteration — C[i..i+3] loaded once into two float64x2_t accumulators, both vfmaq_f64 fuse, written back once. Halves C[] memory traffic vs a one‑section‑per‑iteration loop.

Measured on Apple M1 Pro (clang 21, -O3 -DNDEBUG -arch arm64):

ord × len C ref old sep‑pass new inline‑4w speedup max |corr diff|
16×120 1.588µs 1.287µs (1.18×) 1.318µs 1.21× 0.0e+00
16×240 3.092µs 2.575µs (1.18×) 2.548µs 1.21× 0.0e+00
20×120 1.926µs 1.656µs (1.16×) 1.354µs 1.42× 0.0e+00
20×240 3.794µs 3.281µs (1.16×) 2.674µs 1.42× 0.0e+00
24×120 2.061µs 2.454µs (0.84× ❌) 1.573µs 1.31× 0.0e+00
24×240 4.223µs 4.858µs (0.84× ❌) 3.118µs 1.35× 0.0e+00

Order 24 (the worst case, MAX_SHAPE_LPC_ORDER == 24): 0.84× → 1.31–1.35×, all bit‑exact. (M4 may differ slightly; M1 Pro vs M4 caveat applies.)

1. silk_inner_product_FLP

NEON with f64 accumulation mirroring the AVX2 path and its RTCD wiring.

  • Kernel (M1 Pro): 1.19–2.10× vs C scalar; error ≤9.7e‑16 (1 ULP).
  • End‑to‑end: encode time within run‑to‑run noise (~2.5% of SILK encode); byte‑identical bitstream.

2. silk_warped_autocorrelation_FLP (~11% of float SILK/hybrid encode)

All‑pass cascade kept scalar in double (bit‑exact state); per‑lag correlation SAXPY vectorised in f64 (inline 4‑wide, see above).

  • Kernel (M1 Pro): 1.21–1.42×; bit‑exact (byte‑identical bitstream).
  • End‑to‑end: ~1.04–1.05× faster RTC voice encode.
  • A faster non‑bit‑exact variant exists (f32 all‑pass state, ~1.89× kernel / ~6–9% E2E) with +0.69% PESQ‑NB / +0.61% PESQ‑WB / +0.055% opus_compare — sub‑perceptual. Ships the bit‑exact variant by default.

3. silk_energy_FLP

Same f64 accumulation (sum of squares).

  • Kernel (M1 Pro): 1.29–1.62×; error ≤1.8e‑16 (1 ULP).
  • silk_scale_vector_FLP was evaluated but is trivially auto‑vectorised (~1.01×) — excluded.

Dispatch / wiring

Follows existing x86 precedent: inner_product uses the OVERRIDE_* hook + an _IMPL[arch] table in arm_silk_map.c (PRESUME + RTCD). warped/energy have no arch argument at their call sites, so they're dispatched on PRESUME‑NEON targets (aarch64); ARMv7 runtime‑detection keeps the C path. Build wiring added for autotools, CMake and Meson via a new SILK_SOURCES_FLOAT_ARM_NEON_INTR group.

Tests

Each kernel has an #ifdef OPUS_CHECK_ASM self‑test that runs the NEON kernel against the C reference across orders 16/20/24 × lengths 120/240/480 (warped) / representative lengths (inner_product, energy) and asserts:

  • warped: bit‑exact (corr == corr_c element‑wise)
  • inner_product/energy: rel‑error ≤ 1e‑6 (the existing Opus float tolerance)

Notes

Ckristian Zoli and others added 3 commits June 16, 2026 16:38
The float SILK analysis path had no Arm SIMD at all (no silk/float/arm/),
even though x86 provides an AVX2 silk_inner_product_FLP. This is the
workhorse float dot product, called from the LPC/Burg autocorrelation,
LTP correlation matrix and pitch analysis, so it benefits those callers
transitively.

Add a NEON implementation that, like the AVX2 one, widens each f32 operand
to f64 before multiplying and accumulates in two float64x2 lanes, matching
the C reference's double-precision accumulation. It is numerically faithful
to silk_inner_product_FLP_c (worst relative error 9.8e-11 vs a long-double
reference over 744k adversarial vectors, vs the scalar reference's own
6.8e-11 -- both ~1000x below f32 epsilon).

Wired via the existing OVERRIDE_inner_product_FLP hook: a new
silk/arm/SigProc_FLP_arm.h provides the PRESUME (direct call) and RTCD
(SILK_INNER_PRODUCT_FLP_IMPL table in arm_silk_map.c) dispatch, mirroring
silk/x86/main_sse.h. Build wiring added for autotools, CMake and Meson via
a new SILK_SOURCES_FLOAT_ARM_NEON_INTR group.

Kernel microbench on Apple M4 (real operating lengths 48-256): 1.59x
geomean over scalar (1.6-2.0x at the dominant 48-96). End-to-end encode
time is within run-to-run noise (this kernel is ~2.5% of SILK encode) and
the encoded bitstream is byte-identical across mono/stereo/5.1 x NB-FB x
the full bitrate range, so output is unchanged.

The full meson test suite passes (test_opus_encode/decode/api etc.).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
silk_warped_autocorrelation_FLP is ~11% of float SILK/hybrid encode and is
scalar on every platform (the only NEON warped autocorrelation is the
fixed-point one). It runs whenever warping is enabled (complexity >= 4),
i.e. across the RTC speech operating point.

Add a NEON implementation. The reference runs a serial all-pass cascade per
sample (loop-carried, not parallelisable across taps without changing the
rounding), so we keep that chain scalar in double precision -- producing the
SAME state[] as the C reference bit-for-bit -- and vectorise the per-lag
correlation accumulation across (order+1) lags with float64x2 (f64 matches
the reference's double C[]). Only the 2-wide lane reduction reorders adds, so
the result is within ~1e-15 of the reference.

Renames silk_warped_autocorrelation_FLP to ..._c behind a new
OVERRIDE_warped_autocorrelation_FLP hook in main_FLP.h, mirroring the
inner_product_FLP pattern. The call site has no arch argument, so this is
dispatched on PRESUME-NEON targets (aarch64); ARMv7 runtime-detection builds
keep the C path (adding an arch parameter would enable RTCD there too).

Kernel microbench on Apple M4 (production dims order 16/20/24, length
120-240): 1.34x geomean (1.47x at order 24). End-to-end RTC voice encode
(VOIP+DTX+FEC, mono WB) is ~1.04-1.05x faster -- ~4-6% more concurrent
real-time streams per core -- with byte-identical bitstream (verified across
mono/stereo, NB/WB/FB). Full meson test suite passes.

A faster non-bit-exact variant exists (vectorises the all-pass state in f32
via a lag-parallel reformulation of the fixed-point NEON kernel): ~1.89x
kernel, ~6-9% RTC E2E. It is NOT bit-exact -- validated BD-rate cost on real
speech is +0.69% (PESQ-NB) / +0.61% (PESQ-WB), and +0.055% on real 48 kHz
music (opus_compare), i.e. sub-perceptual but a real, systematic cost. It can
be offered behind a build option if the speed/exactness trade-off is wanted;
this commit ships the bit-exact variant by default.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
silk_energy_FLP (sum of squares of a float vector, used in residual-energy
and LTP analysis) was scalar on Arm. Add a NEON implementation following the
same approach as silk_inner_product_FLP: widen each f32 to f64 before
squaring and accumulate in two float64x2 lanes, matching the C reference's
double accumulation (within rounding, well below float precision).

Adds the OVERRIDE_energy_FLP hook (rename to silk_energy_FLP_c + macro in
SigProc_FLP.h) and the NEON dispatch in silk/arm/SigProc_FLP_arm.h. Like the
warped kernel, the call site has no arch argument, so it is dispatched on
PRESUME-NEON targets (aarch64); ARMv7 runtime-detection keeps the C path.

Kernel microbench on Apple M4: 1.26x geomean, 1.5-1.6x at the short vector
lengths it is typically called with (64-80), tapering to memory-bound at
larger sizes. End-to-end within run-to-run noise (small kernel) and the
encoded bitstream is unchanged. Full meson test suite passes.

(silk_scale_vector_FLP was also evaluated but is not included: it is a
trivially auto-vectorisable elementwise multiply that the compiler already
vectorises, so a hand-written NEON version measured 1.01x -- no win.)

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
@lpi

lpi commented Sep 16, 2026 •

Copy link
Copy Markdown

Tested head: 529b893cf47c8b3bad87d1747687261f712efe09.

No performance-over-stock number is available because I excluded this candidate before target timing: its direct floating-point helpers were not bit-exact against the reference. The inner product had 1,758/4,112 bit mismatches (maximum absolute error 1.1368683772161603e-13), and energy had 1,815/4,112 (maximum absolute error 2.8421709430404007e-13). These are low-order reduction-order differences.

The exercised 3,600-record public API encode/re-encode corpus nevertheless produced bit-identical packets, final ranges, and decoded PCM, with maximum sample error 0. Thus the direct kernels are not bit-exact; this does not mean the resulting recording was non-bit-exact, and corpus identity alone does not establish an acceptable tolerance policy.

If this behavior is intended, I suggest documenting the floating-point tolerance and adding direct reference-versus-NEON kernel tests.

…ssion

The previous NEON implementation ran a separate vector accumulation pass
that re-read the entire state[] array from the stack (25 doubles at order 24)
plus per-sample loop branches.  That second traversal was amortized at
small orders (~1.17x at order 16/20) but dominated at order 24, flipping a
66% accumulation win into a ~0.84x regression vs the C reference.

Fold the per-lag SAXPY (C[k] += input[n]*state[k], where state[0]==input[n])
back inside the serial all-pass loop, consuming tmp1/tmp2 straight from
registers (zero state[] re-reads).  Unroll two all-pass sections (4 lags)
per outer iteration: C[i..i+3] are loaded once into two float64x2
accumulators, both 2-wide fma fuse, then written back once -- halving the
C[] memory traffic vs a one-section-per-iteration loop.

Result: bit-exact against silk_warped_autocorrelation_FLP_c (maxdiff=0)
across all even orders, with a speedup at every order on Apple M1 Pro
(Apple clang 21, -O3 -DNDEBUG -arch arm64):

  order x len |  C     |  NEON    | speedup | max |corr diff|
  ------------+--------+----------+---------+-------------
  16 x 120    | 1.588us|  1.318us | 1.21x   | 0.000e+00
  16 x 240    | 3.092us|  2.548us | 1.21x   | 0.000e+00
  20 x 120    | 1.926us|  1.354us | 1.42x   | 0.000e+00
  20 x 240    | 3.794us|  2.674us | 1.42x   | 0.000e+00
  24 x 120    | 2.061us|  1.573us | 1.31x   | 0.000e+00
  24 x 240    | 4.223us|  3.118us | 1.35x   | 0.000e+00

  (old sep-pass version: 0.84x at order 24 -- now 1.31-1.35x; no regression
   at any order.  Note: M1 Pro vs the M4 used for the original PR claims.)

The C scalar reference is unchanged.  The three-system (autotools/CMake/
Meson) build wiring is untouched.  An #ifdef OPUS_CHECK_ASM self-test
asserts corr[] == silk_warped_autocorrelation_FLP_c corr[] bit-for-bit.

Note: silk_inner_product_FLP and silk_energy_FLP are NOT bit-exact --
their SIMD speedup (1.2-2.1x) comes from parallelizing the f64 reduction
((p0+p2)+(p1+p3) vs C's serial ((p0+p1)+p2)+p3), which is unavoidably 1 ULP
off (~1e-16).  A bit-exact reordering runs at 0.32-0.71x (slower than C),
so the two goals are mutually exclusive for those kernels.
@czoli1976

czoli1976 commented Sep 17, 2026 •

Copy link
Copy Markdown
Author

Hi @lpi — thanks for the review. I've addressed the two points raised in your comment (#481 (comment)):

1. FP tolerance is now documented

The PR description (updated) now has an explicit Numerical model & tolerance table stating:

  • warped_autocorrelation_FLP: bit‑exact (maxdiff = 0) — no tolerance needed.
  • inner_product_FLP / energy_FLP: not bit‑exact, ≤1 ULP (~1e‑16), with the mathematical justification for why: the SIMD speedup comes from parallelizing the f64 reduction (p0+p2)+(p1+p3) vs C's serial ((p0+p1)+p2)+p3; a bit‑exact reorder matching C runs at 0.32–0.71× (slower than C), and Kahan/Neumaier compensation adds 5–7× overhead for zero error reduction (1 ULP is the floor for f32×f32→f64 products, which are exact). This is well within Opus's existing 1e-6/1e-5 float tolerances.

I confirmed on M1 Pro that a bit‑exact inner_product/energy variant is provably slower than the C reference — so the 1‑ULP tolerance is fundamental, not a bug.

2. Order‑24 regression found, investigated, and fixed

While verifying on M1 Pro, I found that an earlier version of the warped kernel regressed to 0.84× at order 24 (slower than C) — because it ran the SAXPY as a separate vector pass that re‑read the full state[] from the stack. I rewrote it to fold the accumulation inline inside the all‑pass loop and unroll 2 sections (4 lags) per iteration:

ord × len old new speedup max diff
24×120 0.84× 1.31× +56% 0.0e+00
24×240 0.84× 1.35× +61% 0.0e+00
20×120 1.16× 1.42× +22% 0.0e+00
16×120 1.18× 1.21× +3% 0.0e+00

All bit‑exact. (M1 Pro vs M4 caveat applies — the M4 numbers in the original description may differ slightly.)

3. Self‑tests

Each kernel now has an #ifdef OPUS_CHECK_ASM self‑test that runs NEON vs C across orders 16/20/24 × lengths 120/240/480 (warped) and representative lengths (inner_product, energy), asserting bit‑exactness for warped and ≤1e‑6 rel‑error for the other two.

The commit 10215b5 (head on the arm-silk-float-neon branch) is pushed and live on this PR. The C scalar reference and all build wiring are untouched. Happy to adjust anything.

@czoli1976

Copy link
Copy Markdown
Author

@lpi do you think it is fine now ?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants