opus: ARMv4 assembly for the backward MDCT inner loops

Cuts realtime decode on the Sansa e200v1 from 52.1 MHz to 50.8 MHz.
clt_mdct_backward was the largest remaining item at 13.5% of decode.

Only the three inner loops move to assembly.  The setup stays in C, so
mdct.c remains readable and the assembly needs no knowledge of
mdct_lookup.

What the compiled loops lose is registers.  Each needs more live values
than gcc can hold, so it spills the loop-invariant pointers, strides and
limits and reloads them every pass: five stack accesses per iteration in
the post-rotation alone.  Holding the twiddle as a 16-bit value and
accumulating the product pair with smull/smlal is what makes the
bookkeeping fit, needing seven live registers where the shifted
MULT16_32_Q15 form needs nine.

ldm/stm helps only where the addressing allows.  The post-rotation walks
the buffer from both ends and so reads and writes contiguous pairs.  The
pre-rotation reads the spectrum through a runtime stride and writes
through the bitrev table, so only its 8-byte output pair merges, and the
TDAC mirror merges nothing.

Over 160 ms of stereo music, traced under qemu:

  clt_mdct_backward  1,037,962 ->   900,982   -13.2%
  whole decode       7,695,876 -> 7,558,896    -1.8%
  loads                650,157 ->   611,667    -5.9%
  stores               350,605 ->   323,605    -7.7%
  multiplies           337,493 ->   337,493   unchanged

Accuracy improves substantially, because all three loops keep 32 bits of
each Q15 product where MULT16_32_Q15_armv4 drops the low bit, and the
backward MDCT applies three such rounds per sample.  The rounding SNR of
the backward transform rises about 9.5 dB, and its worst case error falls
from 708 to 186.  Decoded output differs from the previous build in 90 of
15,360 samples, each by one LSB.

Build with OPUS_ARM_NO_MDCT_ASM to select the C loops instead.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Change-Id: I3c4404b4dbe581d8bcf1f357266a658a068fdb50
This commit is contained in:
Michael Giacomelli 2026-09-11 22:04:48 -04:00 • committed by Solomon Peachy
parent b18f5d6d65
commit d514ee9282
4 changed files with 264 additions and 0 deletions

View file

@ -10,6 +10,7 @@ celt/entenc.c
celt/kiss_fft.c
#if defined(CPU_ARM) && (ARM_ARCH == 4)
celt/arm/kiss_fft_armv4_asm.S
celt/arm/mdct_armv4_asm.S
#endif
celt/laplace.c
celt/mathops.c

View file

@ -0,0 +1,60 @@
/* Copyright (c) 2025 Xiph.Org Foundation and contributors
Copyright (c) 2026 Michael Giacomelli
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#ifndef MDCT_ARMv4_H
#define MDCT_ARMv4_H
/* Hand-written inner loops for the backward MDCT. Only the three loops are
replaced; the setup around them stays in C, so mdct.c remains readable and
the assembly needs no knowledge of mdct_lookup. See
celt/arm/mdct_armv4_asm.S for why the compiled loops lose.
Building with OPUS_ARM_NO_MDCT_ASM selects the C loops instead, which is
how the two are compared. */
#if defined(OPUS_ARM_INLINE_ASM) && defined(FIXED_POINT) \
&& !defined(OPUS_ARM_NO_MDCT_ASM)
#define OVERRIDE_MDCT_PREROT
#define OVERRIDE_MDCT_POSTROT
#define OVERRIDE_MDCT_MIRROR
/* step is a byte stride, so the caller scales by sizeof(kiss_fft_scalar). */
void mdct_prerot_armv4(const kiss_fft_scalar *xp1,
const kiss_fft_scalar *xp2,
const kiss_twiddle_scalar *t,
const opus_int16 *bitrev,
kiss_fft_scalar *yp, int N4, int step);
void mdct_postrot_armv4(kiss_fft_scalar *yp0, kiss_fft_scalar *yp1,
const kiss_twiddle_scalar *t, int N4, int count);
void mdct_mirror_armv4(kiss_fft_scalar *xp1, kiss_fft_scalar *yp1,
const opus_val16 *wp1, const opus_val16 *wp2,
int count);
#endif
#endif /* MDCT_ARMv4_H */

View file

@ -0,0 +1,187 @@
/* ARMv4 inner loops for the CELT backward MDCT.
*
* Copyright (c) 2007-2008 CSIRO
* Copyright (c) 2007-2008 Xiph.Org Foundation
* Copyright (c) 2008 Gregory Maxwell
* Copyright (c) 2026 Michael Giacomelli
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the conditions stated in
* celt/mdct.c are met.
*
* Why this exists
* ---------------
* The three loops here are the whole cost of clt_mdct_backward, and what the
* compiled versions lose is registers. Each loop needs more live values
* than gcc can hold, so it spills the loop-invariant pointers, strides and
* limits and reloads them every pass: five stack accesses per iteration in
* the post-rotation alone. Here the bookkeeping stays resident, and the
* only reload is the pre-rotation's input stride, which buys the seventh
* working register the complex multiply needs.
*
* Holding the twiddle as a 16-bit value and accumulating the product pair
* with smull/smlal is what makes that fit: seven live registers, where the
* shifted MULT16_32_Q15 form needs nine. The ARM7TDMI datasheet also
* promises a 16-bit Rs will finish in two multiplier passes rather than
* four, but device measurement shows no sign of that paying here. Treat
* the form as register economy, not as a multiplier win.
*
* ldm/stm helps only where the addressing allows it. The post-rotation
* reads and writes contiguous pairs, so it bursts. The pre-rotation reads
* the spectrum through a runtime stride and writes through the bit-reversal
* table, so only its 8-byte output pair merges, and the TDAC mirror walks
* two pointers in opposite directions and merges nothing.
*
* All three keep all 32 bits of each Q15 product, where MULT16_32_Q15_armv4
* drops the low bit for speed. The backward MDCT applies three such rounds
* per sample, so this is worth about 9.5 dB rather than being incidental.
*/
#if defined(__thumb__) || defined(__thumb2__)
#error "mdct_armv4_asm.S must be assembled in ARM mode"
#endif
.text
.align 2
/* One Q15 complex multiply, the same register discipline as the C_MUL asm in
arm/kiss_fft_armv4.h. a is (\ar, \ai), b is (\br, \bi), \tt is scratch and
the 64-bit low word. Results: \mr = a.r*b.r - a.i*b.i, \mi = a.r*b.i +
a.i*b.r. \bi and \br are destroyed. \mr must differ from \ar and \br,
\mi from \ai and \tt. */
.macro CMULQ15 ar, ai, br, bi, tt, mr, mi
smull \tt, \mi, \ai, \br
smlal \tt, \mi, \ar, \bi
mov \tt, \tt, lsr #15
orr \mi, \tt, \mi, lsl #17
rsb \bi, \bi, #0
smull \br, \mr, \ar, \br
smlal \br, \mr, \ai, \bi
mov \br, \br, lsr #15
orr \mr, \br, \mr, lsl #17
.endm
/* ------------------------------------------------------------------------
* void mdct_prerot_armv4(const kiss_fft_scalar *xp1,
* const kiss_fft_scalar *xp2,
* const kiss_twiddle_scalar *t,
* const opus_int16 *bitrev,
* kiss_fft_scalar *yp, int N4, int step)
*
* yr = xp2*t[i] + xp1*t[N4+i]; yi = xp1*t[i] - xp2*t[N4+i]
* yp[2*rev+1] = yr; yp[2*rev] = yi; xp1 += step; xp2 -= step
*
* step is in bytes. Both halves of the output are one 8-byte pair, so they
* store together even though the address is scattered by the bitrev table.
* ------------------------------------------------------------------------ */
.global mdct_prerot_armv4
.type mdct_prerot_armv4, %function
mdct_prerot_armv4:
push {r4-r11, lr}
sub sp, sp, #12
ldr r5, [sp, #48] @ yp
ldr r6, [sp, #52] @ N4, also the loop count
ldr r4, [sp, #56] @ step, in bytes
str r4, [sp] @ the one loop-invariant reload
mov r4, r6, lsl #1 @ byte offset from t[i] to t[N4+i]
cmp r6, #0
ble .Lpre_done
.Lpre_loop:
ldr r7, [sp]
ldr r8, [r0], r7 @ x1 = *xp1++
ldr r9, [r1], -r7 @ x2 = *xp2--
ldrsh r11, [r2, r4] @ t[N4+i]
ldrsh r10, [r2], #2 @ t[i]
CMULQ15 r8, r9, r10, r11, lr, r7, ip
ldrsh r10, [r3], #2 @ rev
add r10, r5, r10, lsl #3
stm r10, {r7, ip} @ yi at 2*rev, yr at 2*rev+1
subs r6, r6, #1
bne .Lpre_loop
.Lpre_done:
add sp, sp, #12
pop {r4-r11, pc}
.size mdct_prerot_armv4, .-mdct_prerot_armv4
/* ------------------------------------------------------------------------
* void mdct_postrot_armv4(kiss_fft_scalar *yp0, kiss_fft_scalar *yp1,
* const kiss_twiddle_scalar *t, int N4, int count)
*
* Walks the buffer from both ends at once, so each end contributes a
* contiguous complex pair: two ldm in, and the writes pair up as well, apart
* from yp0[0], which is stored early to free a register before the second
* complex multiply.
*
* The read of yp1 must happen before yp1[1] is written, which is why the
* store order here looks less symmetric than the arithmetic.
* ------------------------------------------------------------------------ */
.align 2
.global mdct_postrot_armv4
.type mdct_postrot_armv4, %function
mdct_postrot_armv4:
push {r4-r11, lr}
ldr r5, [sp, #36] @ count
mov r4, r3, lsl #1 @ byte offset from t[k] to t[N4+k]
sub r6, r3, #1
add r3, r2, r6, lsl #1 @ tb = &t[N4-1], walks down
cmp r5, #0
ble .Lpost_done
.Lpost_loop:
ldrsh r7, [r2, r4] @ t[N4+i]
ldrsh r6, [r2], #2 @ t[i]
ldm r0, {r8, r9} @ im = yp0[0], re = yp0[1]
/* yr = re*t0 + im*t1, yi = re*t1 - im*t0 */
CMULQ15 r9, r8, r7, r6, lr, r11, ip
ldm r1, {r8, r9} @ im = yp1[0], re = yp1[1]
str ip, [r0] @ yp0[0] = yr
ldrsh r7, [r3, r4] @ t[N2-i-1]
ldrsh r6, [r3], #-2 @ t[N4-i-1]
/* mr is yi' and mi is yr', so the pair lands with yr' in the lower
register and stores straight into yp1. */
CMULQ15 r9, r8, r7, r6, lr, ip, r10
stm r1, {r10, r11} @ yp1[0] = yr', yp1[1] = yi
str ip, [r0, #4] @ yp0[1] = yi'
add r0, r0, #8
sub r1, r1, #8
subs r5, r5, #1
bne .Lpost_loop
.Lpost_done:
pop {r4-r11, pc}
.size mdct_postrot_armv4, .-mdct_postrot_armv4
/* ------------------------------------------------------------------------
* void mdct_mirror_armv4(kiss_fft_scalar *xp1, kiss_fft_scalar *yp1,
* const opus_val16 *wp1, const opus_val16 *wp2,
* int count)
*
* The TDAC fold. Two data pointers walk in opposite directions and two
* window pointers do the same, so nothing here is contiguous and nothing
* bursts; the win is entirely the multiplier operand form.
* ------------------------------------------------------------------------ */
.align 2
.global mdct_mirror_armv4
.type mdct_mirror_armv4, %function
mdct_mirror_armv4:
push {r4-r11, lr}
ldr r4, [sp, #36] @ count
cmp r4, #0
ble .Lmir_done
.Lmir_loop:
ldrsh r5, [r3], #-2 @ w2 = *wp2--
ldrsh r6, [r2], #2 @ w1 = *wp1++
ldr r7, [r0] @ x1 = *xp1
ldr r8, [r1] @ x2 = *yp1
/* lo = w2*x2 - w1*x1, hi = w1*x2 + w2*x1 */
CMULQ15 r8, r7, r5, r6, lr, r9, ip
str r9, [r1], #4 @ *yp1++ = lo
str ip, [r0], #-4 @ *xp1-- = hi
subs r4, r4, #1
bne .Lmir_loop
.Lmir_done:
pop {r4-r11, pc}
.size mdct_mirror_armv4, .-mdct_mirror_armv4
.section .note.GNU-stack,"",%progbits

View file

@ -51,6 +51,9 @@
#include <math.h>
#include "os_support.h"
#include "mathops.h"
#if defined(OPUS_ARM_ASM)
#include "arm/mdct_armv4.h"
#endif
#include "stack_alloc.h"
#if defined(MIPSr1_ASM)
@ -265,6 +268,10 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
kiss_fft_scalar * OPUS_RESTRICT yp = out+(overlap>>1);
const kiss_twiddle_scalar * OPUS_RESTRICT t = &trig[0];
const opus_int16 * OPUS_RESTRICT bitrev = l->kfft[shift]->bitrev;
#ifdef OVERRIDE_MDCT_PREROT
mdct_prerot_armv4(xp1, xp2, t, bitrev, yp, N4,
2*stride*(int)sizeof(kiss_fft_scalar));
#else
for(i=0;i<N4;i++)
{
int rev;
@ -279,6 +286,7 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
xp1+=2*stride;
xp2-=2*stride;
}
#endif
}
opus_fft_impl(l->kfft[shift], (kiss_fft_cpx*)(out+(overlap>>1)));
@ -291,6 +299,9 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
const kiss_twiddle_scalar *t = &trig[0];
/* Loop to (N4+1)>>1 to handle odd N4. When N4 is odd, the
middle pair will be computed twice. */
#ifdef OVERRIDE_MDCT_POSTROT
mdct_postrot_armv4(yp0, yp1, t, N4, (N4+1)>>1);
#else
for(i=0;i<(N4+1)>>1;i++)
{
kiss_fft_scalar re, im, yr, yi;
@ -319,6 +330,7 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
yp0 += 2;
yp1 -= 2;
}
#endif
}
/* Mirror on both sides for TDAC */
@ -328,6 +340,9 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
const opus_val16 * OPUS_RESTRICT wp1 = window;
const opus_val16 * OPUS_RESTRICT wp2 = window+overlap-1;
#ifdef OVERRIDE_MDCT_MIRROR
mdct_mirror_armv4(xp1, yp1, wp1, wp2, overlap/2);
#else
for(i = 0; i < overlap/2; i++)
{
kiss_fft_scalar x1, x2;
@ -338,6 +353,7 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
wp1++;
wp2--;
}
#endif
}
}
#endif /* OVERRIDE_clt_mdct_backward */