len0rd/rockbox
1
0
Fork 1
mirror of https://github.com/Rockbox/rockbox.git synced 2026-10-10 08:03:04 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

opus: Good-Thomas FFT for the backward MDCT

Browse source 
Every 48 kHz CELT length is 15 times a power of two and the factors are
co-prime, so the inter-stage twiddles -- 73% of the FFT multiplies at
N=480 -- vanish.  New pfa_fft15 and mdct_postrot_pfa kernels on both
cores; accuracy also improves 0.3 to 0.4 dB against opusdec.

Modelled: -6.17% ARMv4, -1.82% ARMv5E.
Measured: e200v1 42.33 -> 40.10 MHz, Clip+ 29.30 -> 28.24 MHz.
Rescheduling these kernels, folded in here, measured a further -0.75% on
e200v1 and -0.39% on Clip+.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Change-Id: Ifa4a30d1045e905522828958949954615faa440d
This commit is contained in:
Michael Giacomelli 2026-09-18 11:59:28 -04:00 • committed by Solomon Peachy
parent 7b4d1a7f75
commit da9df96c30
14 changed files with 1488 additions and 14 deletions
Download patch file Download diff file Expand all files Collapse all files

3
lib/rbcodec/codecs/libopus/SOURCES
View file

@ -11,12 +11,14 @@ celt/kiss_fft.c
#if defined(CPU_ARM) && (ARM_ARCH == 4)
celt/arm/kiss_fft_armv4_asm.S
celt/arm/mdct_armv4_asm.S
celt/arm/pfa_armv4_asm.S
celt/arm/comb_filter_armv4_asm.S
celt/arm/denorm_armv4_asm.S
celt/arm/deemph_armv4_asm.S
#elif defined(CPU_ARM) && (ARM_ARCH >= 5) && (ARCH_PROFILE == ARM_PROFILE_CLASSIC)
celt/arm/kiss_fft_armv5e_asm.S
celt/arm/mdct_armv5e_asm.S
celt/arm/pfa_armv5e_asm.S
celt/arm/pitch_armv5e_asm.S
celt/arm/comb_filter_armv5e_asm.S
celt/arm/denorm_armv5e_asm.S
@ -28,6 +30,7 @@ celt/arm/deemph_armv5e_asm.S
celt/laplace.c
celt/mathops.c
celt/mdct.c
celt/pfa.c
celt/modes.c
celt/pitch.c
celt/quant_bands.c

9
lib/rbcodec/codecs/libopus/celt/arm/mdct_armv4.h
View file

@ -44,6 +44,10 @@
#define mdct_prerot_opt mdct_prerot_armv4
#define mdct_postrot_opt mdct_postrot_armv4
#ifdef OPUS_PFA
#define OVERRIDE_MDCT_POSTROT_PFA
#define mdct_postrot_pfa_opt mdct_postrot_pfa_armv4
#endif
#define mdct_mirror_opt mdct_mirror_armv4
/* step is a byte stride, so the caller scales by sizeof(kiss_fft_scalar). */
@ -56,6 +60,11 @@ void mdct_prerot_armv4(const kiss_fft_scalar *xp1,
void mdct_postrot_armv4(kiss_fft_scalar *yp0, kiss_fft_scalar *yp1,
const kiss_twiddle_scalar *t, int N4, int count);
#ifdef OPUS_PFA
void mdct_postrot_pfa_armv4(const kiss_fft_scalar *S, kiss_fft_scalar *yp0,
kiss_fft_scalar *yp1, const kiss_twiddle_scalar *t,
const opus_int16 *pmap, int N4);
#endif
void mdct_mirror_armv4(kiss_fft_scalar *xp1, kiss_fft_scalar *yp1,
const opus_val16 *wp1, const opus_val16 *wp2,

130
lib/rbcodec/codecs/libopus/celt/arm/mdct_armv4_asm.S
View file

@ -39,6 +39,16 @@
#if defined(__thumb__) || defined(__thumb2__)
#error "mdct_armv4_asm.S must be assembled in ARM mode"
#endif
/* These kernels are 2,892 bytes that carry 43% of the cycles an ARMv4 decode
* spends, which makes them the densest thing in the codec to put in IRAM.
* config.h sets OPUS_ARM_ICODE on the targets whose codec IRAM window has
* room; everywhere else this stays in .text. The firmware config.h it reads
* that from is assembly-safe, and libopus/config.h guards the rest of its
* includes with __ASSEMBLER__ so this include costs nothing here. */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
.text
@ -48,7 +58,9 @@
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. */
\mi from \ai and \tt. \tt and \mr may be the same register: \tt is dead
after the fourth instruction and \mr is not written until the sixth, and
sharing them is what lets these loops keep their invariants resident. */
.macro CMULQ15 ar, ai, br, bi, tt, mr, mi
smull \tt, \mi, \ai, \br
smlal \tt, \mi, \ar, \bi
@ -79,28 +91,27 @@
.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]
ldr r5, [sp, #36] @ yp
ldr r6, [sp, #40] @ N4, also the loop count
ldr lr, [sp, #44] @ step, in bytes
mov r4, r6, lsl #1
sub r4, r4, #2 @ t[i] is read post-incremented, so the
@ offset on to t[N4+i] is two less
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]
ldr r8, [r0], lr @ x1 = *xp1++
ldr r9, [r1], -lr @ x2 = *xp2--
ldrsh r10, [r2], #2 @ t[i]
CMULQ15 r8, r9, r10, r11, lr, r7, ip
ldrsh r11, [r2, r4] @ t[N4+i], and cover t[i]'s I cycle
CMULQ15 r8, r9, r10, r11, r7, r7, ip
ldrsh r10, [r3], #2 @ rev
subs r6, r6, #1 @ between the load and the address it
@ forms, where the counter is free
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
@ -118,6 +129,14 @@ mdct_prerot_armv4:
* store order here looks less symmetric than the arithmetic.
* ------------------------------------------------------------------------ */
.align 2
/* Unreachable once the prime factor transform is built, so it gives its IRAM
back rather than sitting resident and unread. Still assembled: a size the
Good-Thomas path does not cover falls through to it, and correct-but-slow
from DRAM is the right failure mode. */
.text
.align 2
.global mdct_postrot_armv4
.type mdct_postrot_armv4, %function
mdct_postrot_armv4:
@ -162,6 +181,10 @@ mdct_postrot_armv4:
* bursts; the win is entirely the multiplier operand form.
* ------------------------------------------------------------------------ */
.align 2
.text
.align 2
.global mdct_mirror_armv4
.type mdct_mirror_armv4, %function
mdct_mirror_armv4:
@ -184,4 +207,83 @@ mdct_mirror_armv4:
pop {r4-r11, pc}
.size mdct_mirror_armv4, .-mdct_mirror_armv4
#ifdef OPUS_PFA
.text
/* ------------------------------------------------------------------------
* void mdct_postrot_pfa_armv4(const kiss_fft_scalar *S,
* kiss_fft_scalar *yp0, kiss_fft_scalar *yp1,
* const kiss_twiddle_scalar *t,
* const opus_int16 *pmap, int N4)
*
* The post-rotation for the Good-Thomas transform. Same arithmetic as
* mdct_postrot_armv4, but the spectrum arrives permuted -- X[k] sits at
* M*(k mod 15) + (k mod M) -- so each end is gathered through pmap, which
* holds that index already scaled to a byte offset. One table serves both
* ends: the ascending half reads it forwards, the descending half backwards.
*
* Gathering is why this reads a different buffer than it writes, and that
* turns out to simplify the kernel rather than complicate it. The
* mixed-radix version has to store yp0[0] early and defer yp1[1] because it
* must read yp1 before overwriting it; here yp0 and yp1 are write-only, so
* both halves of each complex multiply store as soon as they exist and
* nothing stays live across the second one.
*
* No iteration counter: the two pmap pointers walk towards each other and
* the loop ends when they cross, which is the same (N4+1)>>1 iterations for
* the even N4 every CELT size has.
*
* N4*2, the offset from t[i] to t[N4+i], does fit in a register after all:
* the complex multiply's 64-bit scratch and its real result can share one,
* so the gather address register is free by the time the multiply needs a
* temporary. It used to be reloaded from the frame once per end, which cost
* six cycles an iteration for want of that observation.
* ------------------------------------------------------------------------ */
.align 2
.global mdct_postrot_pfa_armv4
.type mdct_postrot_pfa_armv4, %function
mdct_postrot_pfa_armv4:
push {r4-r11, lr}
ldr r5, [sp, #36] @ pmap
ldr r6, [sp, #40] @ N4
mov r7, r6, lsl #1 @ byte offset from t[k] to t[N4+k]
sub r6, r6, #1
add r4, r3, r6, lsl #1 @ tb = &t[N4-1], walks down
add r6, r5, r6, lsl #1 @ pd = &pmap[N4-1], walks down
cmp r5, r6
bhi .Lpost_pfa_done
.Lpost_pfa_loop:
/* ascending end: X[i]. The gather index leads, so the address add and
the ldm each have an unrelated load in front of them, and the pair
lands a cycle before the multiply that wants it. */
ldrsh r12, [r5], #2 @ byte offset of X[i]
ldrsh r11, [r3, r7] @ t1 = t[N4+i]
add r12, r0, r12
ldm r12, {r8, r9} @ im = S[2p], re = S[2p+1]
ldrsh r10, [r3], #2 @ t0 = t[i]
CMULQ15 r9, r8, r11, r10, r12, r12, lr
str lr, [r1] @ yp0[0] = yr
str r12, [r2, #4] @ yp1[1] = yi
/* descending end: X[N4-1-i] */
ldrsh r12, [r6], #-2 @ byte offset of X[N4-1-i]
ldrsh r11, [r4, r7] @ t1 = t[N2-i-1]
add r12, r0, r12
ldm r12, {r8, r9}
ldrsh r10, [r4], #-2 @ t0 = t[N4-i-1]
CMULQ15 r9, r8, r11, r10, r12, r12, lr
str lr, [r2] @ yp1[0] = yr'
str r12, [r1, #4] @ yp0[1] = yi'
add r1, r1, #8
sub r2, r2, #8
cmp r5, r6
bls .Lpost_pfa_loop @ bls, not blo, so an odd N4 still gets
@ its middle pair, as the C loop does
.Lpost_pfa_done:
pop {r4-r11, pc}
.size mdct_postrot_pfa_armv4, .-mdct_postrot_pfa_armv4
#endif /* OPUS_PFA */
.section .note.GNU-stack,"",%progbits

7
lib/rbcodec/codecs/libopus/celt/arm/mdct_armv5e.h
View file

@ -48,6 +48,13 @@
#define mdct_prerot_opt mdct_prerot_armv5e
#define mdct_postrot_opt mdct_postrot_armv5e
#ifdef OPUS_PFA
#define OVERRIDE_MDCT_POSTROT_PFA
#define mdct_postrot_pfa_opt mdct_postrot_pfa_armv5e
void mdct_postrot_pfa_armv5e(const kiss_fft_scalar *S, kiss_fft_scalar *yp0,
kiss_fft_scalar *yp1, const kiss_twiddle_scalar *t,
const opus_int16 *pmap, int N4);
#endif
#define mdct_mirror_opt mdct_mirror_armv5e
/* step is a byte stride, so the caller scales by sizeof(kiss_fft_scalar). */

165
lib/rbcodec/codecs/libopus/celt/arm/mdct_armv5e_asm.S
View file

@ -300,5 +300,170 @@ mdct_mirror_armv5e:
pop {r4-r11, pc}
.size mdct_mirror_armv5e, .-mdct_mirror_armv5e
#ifdef OPUS_PFA
.text
.align 2
/* One post-rotated output.
*
* \sel picks which half of the packed twiddle pair to use
* \pm the pmap pointer, post-incremented by \pi
* \b0/\o0 where yr goes, \b1/\o1 where yi goes
*
* r10 holds t[k] packed with its neighbour and r11 holds t[N4+k] with its
* neighbour, so two loads serve two outputs. The gather is an ldrsh of the
* byte offset and one ldm, which fetches both halves of the complex value in
* a single transaction.
*/
/* The gather index arrives in r12, and this body loads the next one before
it needs its own loaded pair -- an index is a halfword, which interlocks
for two cycles if the address add follows it, and the ldm interlocks for
one if the multiply does. Rotating the load one body ahead covers both,
and it is the only place to put it: every other register is live.
That in turn is why yr stores before yi is computed. Holding the index
in r12 costs the body its second scratch, so the accumulator has to be
free again before the imaginary half starts.
\npm and \npi name the pointer the NEXT body reads, or none at the end of
the chain, where the loop head reloads instead. */
.macro PR5_PFA sel, npm, npi, b0, o0, b1, o1
add r12, r0, r12
ldm r12, {r9, lr} @ r9 = im = S[2p], lr = re = S[2p+1];
@ ldm fills in register order, not the
@ order written, so r9 takes the lower
@ address whichever way the list reads
.ifnc \npm, none
ldrsh r12, [\npm], #\npi @ the next body's index
.endif
smulw\sel r8, lr, r10 @ re * t0
smlaw\sel r8, r9, r11, r8 @ + im * t1 -> yr at Q16
mov r8, r8, lsl #1
str r8, [\b0, #\o0]
smulw\sel r8, lr, r11 @ re * t1
smulw\sel lr, r9, r10 @ im * t0
sub r8, r8, lr
mov r8, r8, lsl #1
str r8, [\b1, #\o1]
.endm
/* The same body for the scalar fallback, which loads its own index: it runs
one output per end, so there is no next body to rotate against. */
.macro PR5_PFA_S sel, pm, pi, b0, o0, b1, o1
ldrsh r12, [\pm], #\pi
add r12, r0, r12
ldm r12, {r9, lr}
smulw\sel r8, lr, r10
smlaw\sel r8, r9, r11, r8
smulw\sel r12, lr, r11
smulw\sel lr, r9, r10
mov r8, r8, lsl #1
sub r12, r12, lr
str r8, [\b0, #\o0]
mov r12, r12, lsl #1
str r12, [\b1, #\o1]
.endm
/* ------------------------------------------------------------------------
* void mdct_postrot_pfa_armv5e(const kiss_fft_scalar *S,
* kiss_fft_scalar *yp0, kiss_fft_scalar *yp1,
* const kiss_twiddle_scalar *t,
* const opus_int16 *pmap, int N4)
*
* The post-rotation for the Good-Thomas transform. Same arithmetic as
* mdct_postrot_armv5e; what changes is that the spectrum arrives permuted,
* X[k] at M*(k mod 15) + (k mod M), so each end is gathered through pmap,
* which already holds that index as a byte offset. The ascending half reads
* the table forwards and the descending half backwards.
*
* Gathering makes yp0 and yp1 write-only, and that is what lets the twiddles
* stay packed. A kiss_twiddle_scalar is 16 bits, so t[k] and t[k+1] are one
* aligned word; the mixed-radix kernel exploits that by holding all four
* twiddle words -- both ends at once -- across a body that does two outputs
* from each. Here two of those registers are spent on the pmap pointers, so
* instead the body does both ascending outputs, reloads the pair for the
* descending end, and does both of those. That reordering is legal only
* because nothing is read back: in the in-place version the descending
* output has to be read before the ascending one overwrites it, which forces
* the two ends to interleave.
*
* No iteration counter either: the two pmap pointers walk towards each other
* and the loop ends when they cross. Four outputs an iteration and N4 a
* multiple of four -- 60, 120, 240, 480 -- means that lands exactly.
*
* The scalar path below is the fallback for a trig table the packed loads
* cannot address as words, or an N4 the unroll does not divide. No CELT
* mode reaches it; build with -DMDCT_ARMV5E_FORCE_SCALAR to route everything
* through it and confirm it still agrees.
* ------------------------------------------------------------------------ */
.align 2
.global mdct_postrot_pfa_armv5e
.type mdct_postrot_pfa_armv5e, %function
mdct_postrot_pfa_armv5e:
push {r4-r11, lr}
ldr r5, [sp, #36] @ pmap
ldr r12, [sp, #40] @ N4
mov r7, r12, lsl #1 @ byte offset from t[k] to t[N4+k]
sub r8, r12, #1
add r4, r3, r8, lsl #1 @ tb = &t[N4-1], walks down
add r6, r5, r8, lsl #1 @ pd = &pmap[N4-1], walks down
/* The packed path needs the trig table word-addressable and N4 a
multiple of four; the second gives the first its stride and lets the
crossing test land exactly. */
orr r8, r3, r12
#ifdef MDCT_ARMV5E_FORCE_SCALAR
orr r8, r8, #1
#endif
tst r8, #3
bne .Lpost5_pfa_scalar
cmp r5, r6
bhi .Lpost5_pfa_done
.Lpost5_pfa_pairs:
ldrsh r12, [r5], #2 @ index of X[i]; the twiddle loads below
@ stand between it and the address add
ldr r10, [r3] @ t[i] : t[i+1]
ldr r11, [r3, r7] @ t[N4+i] : t[N4+i+1]
add r3, r3, #4
PR5_PFA b, r5, 2, r1, 0, r2, 4 @ X[i]
PR5_PFA t, r6, -2, r1, 8, r2, -4 @ X[i+1]
add lr, r4, r7
ldr r10, [r4, #-2] @ t[N4-i-2] : t[N4-i-1]
ldr r11, [lr, #-2] @ t[N2-i-2] : t[N2-i-1]
sub r4, r4, #4
PR5_PFA t, r6, -2, r2, 0, r1, 4 @ X[N4-1-i]
PR5_PFA b, none, 0, r2, -8, r1, 12 @ X[N4-2-i]
add r1, r1, #16
sub r2, r2, #16
cmp r5, r6
blo .Lpost5_pfa_pairs
b .Lpost5_pfa_done
/* One output per end per iteration, twiddles a halfword at a time. The
crossing test is bls rather than blo so an odd N4 still computes its
middle pair, which is what the C loop to (N4+1)>>1 does. */
.Lpost5_pfa_scalar:
cmp r5, r6
bhi .Lpost5_pfa_done
.Lpost5_pfa_sloop:
ldrsh r11, [r3, r7] @ t[N4+i]
ldrsh r10, [r3], #2 @ t[i]
PR5_PFA_S b, r5, 2, r1, 0, r2, 4
ldrsh r11, [r4, r7] @ t[N2-i-1]
ldrsh r10, [r4], #-2 @ t[N4-i-1]
PR5_PFA_S b, r6, -2, r2, 0, r1, 4
add r1, r1, #8
sub r2, r2, #8
cmp r5, r6
bls .Lpost5_pfa_sloop
.Lpost5_pfa_done:
pop {r4-r11, pc}
.size mdct_postrot_pfa_armv5e, .-mdct_postrot_pfa_armv5e
#endif
.section .note.GNU-stack,"",%progbits

30
lib/rbcodec/codecs/libopus/celt/arm/pfa_arm.h Normal file
View file

@ -0,0 +1,30 @@
/* ARM hooks for the Good-Thomas 15-point kernel.
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/kiss_fft.c are met.
*/
#ifndef PFA_ARM_H
#define PFA_ARM_H
#if defined(FIXED_POINT) && defined(OPUS_PFA)
#if defined(OPUS_ARM_INLINE_EDSP)
#define OVERRIDE_PFA_FFT15
void pfa_fft15_armv5e(const kiss_fft_cpx *in, kiss_fft_cpx *out, int ostride);
#define PFA_FFT15(in, out, ostride) pfa_fft15_armv5e(in, out, ostride)
#elif defined(OPUS_ARM_INLINE_ASM)
#define OVERRIDE_PFA_FFT15
void pfa_fft15_armv4(const kiss_fft_cpx *in, kiss_fft_cpx *out, int ostride);
#define PFA_FFT15(in, out, ostride) pfa_fft15_armv4(in, out, ostride)
#endif
#endif /* FIXED_POINT && OPUS_PFA */
#endif /* PFA_ARM_H */

286
lib/rbcodec/codecs/libopus/celt/arm/pfa_armv4_asm.S Normal file
View file

@ -0,0 +1,286 @@
/* ARMv4 Good-Thomas 15-point DFT for the CELT backward MDCT.
*
* 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/kiss_fft.c are met.
*
* Why this exists
* ---------------
* 15 = 3*5 and the two are co-prime, so this is Good-Thomas inside
* Good-Thomas: five 3-point DFTs, then three 5-point DFTs, and nothing
* between them. A 15-point complex DFT for 40 real multiplies, where the
* mixed-radix chain it replaces spends 16 twiddle multiplies per point on
* the radix-5 pass alone.
*
* The arithmetic is the arithmetic the existing butterflies already use, so
* nothing new is being claimed about accuracy: epi3 is the radix-3 constant
* from kf_bfly3, and the radix-5 half keeps the sqrt(5)/4 collapse from
* kf_bfly5, where cos(2pi/5)+cos(4pi/5) is exactly -1/2 -- and the Q15
* constants honour it exactly, 10126 - 26510 == -16384 -- so four cosine
* products become one multiply by yc. Q15 narrowing is the exact
* smull/shift/orr form the other ARMv4 kernels use, which keeps all 32 bits
* of the product where MULT16_32_Q15_armv4 drops the low one.
*
* Register pressure is the whole difficulty. A 5-point DFT needs eight live
* words across its eight multiplies plus two scratch, and ARM has fourteen
* registers to hold that, the output pointer, its stride and two constants.
* It fits only because the sum and difference pairs s5 and s11 are parked on
* the frame across the multiplies -- sixteen bytes, to the same cache line
* every call -- and because yc borrows the register ya.i will need later,
* reloading it afterwards for the price of one ldr per 5-point DFT.
*
* The output permutation costs nothing. Output k2 of the k1'th 5-point DFT
* belongs at 15-point index (10*k1 + 6*k2) mod 15; walking that from k1 in
* steps of 3 visits every one of the five, so the destinations are a single
* post-incremented pointer and the permutation is only a choice of which
* register to store.
*/
#if defined(__thumb__) || defined(__thumb2__)
#error "pfa_armv4_asm.S must be assembled in ARM mode"
#endif
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
.text
.align 2
/* Frame: v[15] complex at 0..119, the s5/s11 park at 120..135, then the
saved output pointer, the 3-element output step and the 5-element stride
the three 5-point bases are derived from. */
#define VOFF 0
#define PARK 120
#define SV_OUT 136
#define SV_STEP3 140
#define SV_FIVE 144
#define FRAME 152
/* Q15 narrowing that keeps every bit of the product: see kiss_fft_armv4.h. */
.macro QNARROW dst, lo, hi
mov \lo, \lo, lsr #15
orr \dst, \lo, \hi, lsl #17
.endm
/* Advance the output pointer by d groups of three elements. The five
destinations of one 5-point DFT are (10*k1 + 6*k2) mod 15, which walking
k2 in the order 0,1,4,2,3 visits as steps of 3 with one wrap, so every
delta is a small multiple of r3 and none needs a multiply. */
.macro PFA_ADV d
.if \d == 2
add r0, r0, r3, lsl #1
.elseif \d == 1
add r0, r0, r3
.elseif \d == -3
sub r0, r0, r3, lsl #1
sub r0, r0, r3
.elseif \d == -4
sub r0, r0, r3, lsl #2
.else
.error "unexpected PFA store delta"
.endif
.endm
/* One 5-point DFT over v[voff..voff+4], writing through r0.
*
* r0 = output, already at this DFT's base r1 = ya.i r2 = yb.i
* r3 = three output elements, in bytes
*
* Ten multiplies: two for the sqrt(5)/4 term that stands in for the four
* cosine products, eight for the sines. The four sine products of each
* half share their operands between s6 and s12, so both come out of one
* pass over s9 and s10 and neither has to be kept alive across the other.
*/
.macro PFA_DFT5 voff, d1, d2, d3, d4
add r12, sp, #(\voff+8)
ldmia r12, {r4,r5,r6,r7,r8,r9} @ x1, x2, x3
add r12, sp, #(\voff+32)
ldm r12, {r10,r11} @ x4
add r12, r4, r10 @ s7.r = x1.r + x4.r
sub r4, r4, r10 @ s10.r = x1.r - x4.r
add r10, r5, r11 @ s7.i
sub r5, r5, r11 @ s10.i
add r11, r6, r8 @ s8.r = x2.r + x3.r
sub r6, r6, r8 @ s9.r = x2.r - x3.r
add r8, r7, r9 @ s8.i
sub r7, r7, r9 @ s9.i
/* s4 *= sqrt(5)/4. yc borrows ya.i's register; one ldr puts it back,
which is cheaper than spending a register on a constant used twice.
It loads ahead of the four sums below, which is what keeps its I
cycle off the multiply. */
ldr r1, .Lpfa_yc
add r9, r12, r11 @ s3.r = s7.r + s8.r
sub r12, r12, r11 @ s4.r = s7.r - s8.r
add r11, r10, r8 @ s3.i
sub r10, r10, r8 @ s4.i
smull r8, lr, r12, r1
QNARROW r12, r8, lr @ s4.r
smull r8, lr, r10, r1
QNARROW r10, r8, lr @ s4.i
ldr lr, [sp, #(\voff+0)] @ x0.r
ldr r1, .Lpfa_yai @ and this covers its I cycle
/* X0 = x0 + s3 goes out immediately, then s3 becomes x0 - s3/4. x0.r
is loaded a line early so the constant reload above covers it, and
x0.i goes into the scratch that storing X0.r frees, which puts three
instructions between that load and its use. */
add r8, lr, r9 @ X0.r
str r8, [r0]
ldr r8, [sp, #(\voff+4)] @ x0.i
add r9, r9, #2
sub r9, lr, r9, asr #2 @ QUARTER_OF rounds, as kf_bfly5 does
add lr, r8, r11 @ X0.i
str lr, [r0, #4]
add r11, r11, #2
sub r11, r8, r11, asr #2
/* s5 and s11 are not touched by the eight multiplies below and there is
no register to keep them in, so they wait on the frame -- sixteen
bytes, the same line every call. */
add r8, r9, r12 @ s5.r
sub r9, r9, r12 @ s11.r
add lr, r11, r10 @ s5.i
sub r11, r11, r10 @ s11.i
mov r10, lr @ stm stores in register order, not
add r12, sp, #PARK @ list order, so s5.i must sit in r10
stmia r12, {r8, r9, r10, r11}
/* s6.r = s10.i*ya.i + s9.i*yb.i s12.r = s9.i*ya.i - s10.i*yb.i */
smull r10, r12, r5, r1
QNARROW r8, r10, r12
smull r10, r12, r7, r2
QNARROW r9, r10, r12
smull r10, r12, r7, r1
QNARROW r11, r10, r12
smull r10, r12, r5, r2
QNARROW lr, r10, r12
add r5, r8, r9 @ s6.r
sub r7, r11, lr @ s12.r
/* s6.i = -(s10.r*ya.i + s9.r*yb.i) s12.i = s10.r*yb.i - s9.r*ya.i */
smull r10, r12, r4, r1
QNARROW r8, r10, r12
smull r10, r12, r6, r2
QNARROW r9, r10, r12
smull r10, r12, r4, r2
QNARROW r11, r10, r12
smull r10, r12, r6, r1
QNARROW lr, r10, r12
add r4, r8, r9
rsb r4, r4, #0 @ s6.i
sub r6, r11, lr @ s12.i
add r12, sp, #PARK
ldmia r12, {r8, r9, r10, r11} @ s5.r, s11.r, s5.i, s11.i
sub r12, r8, r5 @ X1 = s5 - s6
sub lr, r10, r4
PFA_ADV \d1
stmia r0, {r12, lr}
add r12, r8, r5 @ X4 = s5 + s6
add lr, r10, r4
PFA_ADV \d2
stmia r0, {r12, lr}
add r12, r9, r7 @ X2 = s11 + s12
add lr, r11, r6
PFA_ADV \d3
stmia r0, {r12, lr}
sub r12, r9, r7 @ X3 = s11 - s12
sub lr, r11, r6
PFA_ADV \d4
stmia r0, {r12, lr}
.endm
/* ------------------------------------------------------------------------
* void pfa_fft15_armv4(const kiss_fft_cpx *in, kiss_fft_cpx *out,
* int ostride)
*
* r0 = in (15 contiguous complex) r1 = out r2 = ostride, in elements
* ------------------------------------------------------------------------ */
.global pfa_fft15_armv4
.type pfa_fft15_armv4, %function
pfa_fft15_armv4:
push {r4-r11, lr}
sub sp, sp, #FRAME
mov r2, r2, lsl #3 @ ostride in bytes
str r1, [sp, #SV_OUT]
add r3, r2, r2, lsl #1 @ 3 elements
str r3, [sp, #SV_STEP3]
add r3, r2, r2, lsl #2 @ 5 elements
str r3, [sp, #SV_FIVE]
/* ---- Five 3-point DFTs -------------------------------------------------
* Group g reads in[3g..3g+2] and writes v[g], v[5+g], v[10+g], so the
* 5-point DFTs below each read a contiguous run of five.
*/
add r1, sp, #VOFF @ &v[0]
add r2, sp, #(VOFF+40) @ &v[5]
add r3, sp, #(VOFF+80) @ &v[10]
ldr lr, .Lpfa_epi3
mov r12, #5
.Lpfa_b3:
ldmia r0!, {r4,r5,r6,r7,r8,r9} @ x0, x1, x2
add r10, r6, r8 @ s3.r = x1.r + x2.r
sub r6, r6, r8 @ s0.r = x1.r - x2.r
add r11, r7, r9 @ s3.i
sub r7, r7, r9 @ s0.i
sub r8, r4, r10, asr #1 @ xm.r = x0.r - HALF_OF(s3.r)
sub r9, r5, r11, asr #1 @ xm.i
add r4, r4, r10 @ X0.r = x0.r + s3.r
add r5, r5, r11 @ X0.i
smull r10, r11, r6, lr @ s0.r * epi3
QNARROW r6, r10, r11
smull r10, r11, r7, lr @ s0.i * epi3
QNARROW r7, r10, r11
stmia r1!, {r4,r5} @ v[g]
sub r10, r8, r7
add r11, r9, r6
stmia r2!, {r10,r11} @ v[5+g]
add r10, r8, r7
sub r11, r9, r6
stmia r3!, {r10,r11} @ v[10+g]
subs r12, r12, #1
bne .Lpfa_b3
/* ---- Three 5-point DFTs ------------------------------------------------ */
ldr r1, .Lpfa_yai
ldr r2, .Lpfa_ybi
ldr r3, [sp, #SV_STEP3]
ldr r0, [sp, #SV_OUT] @ k1 = 0: base is out itself
PFA_DFT5 (VOFF+0), 2, 1, 1, -3
ldr r0, [sp, #SV_OUT]
ldr r12, [sp, #SV_FIVE]
add r0, r0, r12, lsl #1 @ k1 = 1: base is out + 10 elements
PFA_DFT5 (VOFF+40), -3, 1, 1, 2
ldr r0, [sp, #SV_OUT]
ldr r12, [sp, #SV_FIVE]
add r0, r0, r12 @ k1 = 2: base is out + 5 elements
PFA_DFT5 (VOFF+80), 2, 1, -4, 2
add sp, sp, #FRAME
pop {r4-r11, pc}
.size pfa_fft15_armv4, .-pfa_fft15_armv4
.align 2
.Lpfa_epi3:
.word -28378 @ -sin(2*pi/3), as kf_bfly3 uses
.Lpfa_yai:
.word -31164 @ ya.i
.Lpfa_ybi:
.word -19261 @ yb.i
.Lpfa_yc:
.word 18318 @ sqrt(5)/4

227
lib/rbcodec/codecs/libopus/celt/arm/pfa_armv5e_asm.S Normal file
View file

@ -0,0 +1,227 @@
/* ARMv5E Good-Thomas 15-point DFT for the CELT backward MDCT.
*
* 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/kiss_fft.c are met.
*
* Why this exists
* ---------------
* The same transform as pfa_armv4_asm.S -- five 3-point DFTs then three
* 5-point DFTs, co-prime so nothing sits between them -- but this core
* changes what the arithmetic costs, and with it the register budget.
*
* A Q15 narrowing is one smulwb rather than smull/shift/orr, so each of the
* forty multiplies drops from three instructions to one and needs one
* scratch register rather than a pair. Where the C adds two products,
* smlawb folds the add in and is bit-exact with it, because both truncate
* at Q16 before the sum and the shift distributes. Where it subtracts them
* the two stay separate: smlaw only adds, so fusing would mean negating a
* twiddle, and truncation is not symmetric about zero.
*
* The Q15 doubling is deferred. MULT16_32_Q15 on this core is smulwb
* followed by a shift left, and every consumer of a product here is an add
* or a subtract, so the shift folds into it as a free immediate. That is
* exact rather than approximate: a +/- (b << 1) is the same value as
* a +/- SHL32(b, 1) in all 32 bits.
*
* Between them those two savings free enough registers that the 5-point DFT
* needs no frame at all -- where the ARMv4 version has to park its sum and
* difference pair across the multiplies, this one keeps everything live.
* ya.i and yb.i share one word, picked apart by the b and t halves of
* smulwb and smulwt, which is the same packing the radix-5 pass in
* kiss_fft_armv5e_asm.S uses.
*/
#if defined(__thumb__) || defined(__thumb2__)
#error "pfa_armv5e_asm.S must be assembled in ARM mode"
#endif
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
.text
.align 2
#define VOFF 0
#define SV_OUT 120
#define SV_STEP3 124
#define SV_FIVE 128
#define FRAME 136
/* Advance the output pointer by d groups of three elements; see the ARMv4
kernel for why every delta is a small multiple of one stride. */
.macro PFA5_ADV d
.if \d == 2
add r0, r0, r3, lsl #1
.elseif \d == 1
add r0, r0, r3
.elseif \d == -3
sub r0, r0, r3, lsl #1
sub r0, r0, r3
.elseif \d == -4
sub r0, r0, r3, lsl #2
.else
.error "unexpected PFA store delta"
.endif
.endm
/* One 5-point DFT over v[voff..voff+4].
* r0 = output, at this DFT's base r1 = ya.i packed with yb.i
* r2 = sqrt(5)/4 r3 = three output elements, in bytes
*/
.macro PFA5_DFT5 voff, d1, d2, d3, d4
add r12, sp, #(\voff+8)
ldmia r12, {r4,r5,r6,r7,r8,r9} @ x1, x2, x3
add r12, sp, #(\voff+32)
ldm r12, {r10,r11} @ x4
add r12, r4, r10 @ s7.r = x1.r + x4.r
sub r4, r4, r10 @ s10.r = x1.r - x4.r
add r10, r5, r11 @ s7.i
sub r5, r5, r11 @ s10.i
add r11, r6, r8 @ s8.r = x2.r + x3.r
sub r6, r6, r8 @ s9.r = x2.r - x3.r
add r8, r7, r9 @ s8.i
sub r7, r7, r9 @ s9.i
add r9, r12, r11 @ s3.r = s7.r + s8.r
sub r12, r12, r11 @ s4.r = s7.r - s8.r
add r11, r10, r8 @ s3.i
sub r10, r10, r8 @ s4.i
ldr lr, [sp, #(\voff+0)] @ x0.r, ahead of the two multiplies
smulwb r12, r12, r2 @ s4.r * sqrt(5)/4, left at Q16
smulwb r10, r10, r2 @ s4.i * sqrt(5)/4
add r8, lr, r9 @ X0.r
str r8, [r0]
ldr r8, [sp, #(\voff+4)] @ x0.i, into the scratch that frees
add r9, r9, #2
sub r9, lr, r9, asr #2 @ QUARTER_OF rounds, as kf_bfly5 does
add lr, r8, r11 @ X0.i
str lr, [r0, #4]
add r11, r11, #2
sub r11, r8, r11, asr #2
add r8, r9, r12, lsl #1 @ s5.r (the deferred doubling)
sub r9, r9, r12, lsl #1 @ s11.r
add lr, r11, r10, lsl #1 @ s5.i
sub r11, r11, r10, lsl #1 @ s11.i
/* s6.r = s10.i*ya + s9.i*yb, so smlawt carries the add. */
smulwb r10, r5, r1
smlawt r10, r7, r1, r10
/* s12.r = s9.i*ya - s10.i*yb: a difference, so the products stay apart. */
smulwb r12, r7, r1
smulwt r5, r5, r1
sub r12, r12, r5
/* s6.i is the negation of this sum; the sign folds into the consumers. */
smulwb r5, r4, r1
smlawt r5, r6, r1, r5
/* s12.i = s10.r*yb - s9.r*ya */
smulwt r4, r4, r1
smulwb r6, r6, r1
sub r4, r4, r6
sub r6, r8, r10, lsl #1 @ X1.r = s5.r - s6.r
add r7, lr, r5, lsl #1 @ X1.i = s5.i - s6.i
PFA5_ADV \d1
stmia r0, {r6, r7}
add r6, r8, r10, lsl #1 @ X4.r
sub r7, lr, r5, lsl #1 @ X4.i
PFA5_ADV \d2
stmia r0, {r6, r7}
add r6, r9, r12, lsl #1 @ X2.r = s11.r + s12.r
add r7, r11, r4, lsl #1 @ X2.i
PFA5_ADV \d3
stmia r0, {r6, r7}
sub r6, r9, r12, lsl #1 @ X3.r
sub r7, r11, r4, lsl #1 @ X3.i
PFA5_ADV \d4
stmia r0, {r6, r7}
.endm
/* ------------------------------------------------------------------------
* void pfa_fft15_armv5e(const kiss_fft_cpx *in, kiss_fft_cpx *out,
* int ostride)
* ------------------------------------------------------------------------ */
.global pfa_fft15_armv5e
.type pfa_fft15_armv5e, %function
pfa_fft15_armv5e:
push {r4-r11, lr}
sub sp, sp, #FRAME
mov r2, r2, lsl #3 @ ostride in bytes
str r1, [sp, #SV_OUT]
add r3, r2, r2, lsl #1
str r3, [sp, #SV_STEP3]
add r3, r2, r2, lsl #2
str r3, [sp, #SV_FIVE]
/* ---- Five 3-point DFTs ------------------------------------------------ */
add r1, sp, #VOFF
add r2, sp, #(VOFF+40)
add r3, sp, #(VOFF+80)
ldr lr, .Lpfa5_epi3
mov r12, #5
.Lpfa5_b3:
ldmia r0!, {r4,r5,r6,r7,r8,r9} @ x0, x1, x2
add r10, r6, r8 @ s3.r
sub r6, r6, r8 @ s0.r
add r11, r7, r9 @ s3.i
sub r7, r7, r9 @ s0.i
sub r8, r4, r10, asr #1 @ xm.r = x0.r - HALF_OF(s3.r)
sub r9, r5, r11, asr #1 @ xm.i
add r4, r4, r10 @ X0.r
add r5, r5, r11 @ X0.i
smulwb r6, r6, lr @ s0.r * epi3, left at Q16
smulwb r7, r7, lr @ s0.i * epi3
stmia r1!, {r4,r5} @ v[g]
sub r10, r8, r7, lsl #1
add r11, r9, r6, lsl #1
stmia r2!, {r10,r11} @ v[5+g]
add r10, r8, r7, lsl #1
sub r11, r9, r6, lsl #1
stmia r3!, {r10,r11} @ v[10+g]
subs r12, r12, #1
bne .Lpfa5_b3
/* ---- Three 5-point DFTs ----------------------------------------------- */
ldr r1, .Lpfa5_yab
ldr r2, .Lpfa5_yc
ldr r3, [sp, #SV_STEP3]
ldr r0, [sp, #SV_OUT]
PFA5_DFT5 (VOFF+0), 2, 1, 1, -3
ldr r0, [sp, #SV_OUT]
ldr r12, [sp, #SV_FIVE]
add r0, r0, r12, lsl #1
PFA5_DFT5 (VOFF+40), -3, 1, 1, 2
ldr r0, [sp, #SV_OUT]
ldr r12, [sp, #SV_FIVE]
add r0, r0, r12
PFA5_DFT5 (VOFF+80), 2, 1, -4, 2
add sp, sp, #FRAME
pop {r4-r11, pc}
.size pfa_fft15_armv5e, .-pfa_fft15_armv5e
.align 2
.Lpfa5_epi3:
.word -28378 @ -sin(2*pi/3), in the low half
.Lpfa5_yab:
.word 0xb4c38644 @ ya.i = -31164 low, yb.i = -19261 top
.Lpfa5_yc:
.word 18318 @ sqrt(5)/4
.section .note.GNU-stack,"",%progbits

58
lib/rbcodec/codecs/libopus/celt/kiss_fft.c
View file

@ -40,6 +40,10 @@
#include "os_support.h"
#include "mathops.h"
#include "stack_alloc.h"
#ifdef OPUS_PFA
#include "pfa.h"
#include "pfa_tables.h"
#endif
/* The guts header contains all the multiplication and addition macros that are defined for
complex numbers. It also delares the kf_ internal functions.
@ -619,6 +623,60 @@ void opus_fft_impl(const kiss_fft_state *st,kiss_fft_cpx *fout)
}
}
#ifdef OPUS_PFA
/* kf_bfly4 reads nothing from the state but the twiddle table, and the prime
factor sub-transforms want their own: one 32-entry table serves all four
sizes, read with a stride, exactly as kiss_fft shares a single table
between the four transform lengths. */
static const kiss_fft_state pfa_st32 = { 0, 0, 0, 0, {0}, 0, pfa_tw32, 0 };
void opus_pfa_impl(const kiss_fft_cpx *fin,
kiss_fft_cpx *fout, int nfft)
{
const opus_int16 *brev;
int M = nfft/15;
int n2;
switch (nfft)
{
case 60: brev = pfa_brev_4; break;
case 120: brev = pfa_brev_8; break;
case 240: brev = pfa_brev_16; break;
default: brev = pfa_brev_32; break;
}
/* Pass 1: M fifteen-point DFTs. The bit-reversal the radix-4 chain below
expects is folded into this scatter, so pass 2 needs no permutation of
its own. */
for (n2=0;n2<M;n2++)
PFA_FFT15(fin + 15*n2, fout + brev[n2], M);
/* Pass 2: fifteen M-point FFTs, contiguous and identical, so one call per
stage covers all fifteen with N counting groups across the lot. The
factorisation is the one kf_factor would pick for M, which is why the
existing butterflies serve unchanged. */
switch (M)
{
case 4:
kf_bfly4(fout, 1, &pfa_st32, 1, 15, 4);
break;
case 8:
kf_bfly4(fout, 1, &pfa_st32, 1, 2*15, 4);
kf_bfly2(fout, 4, 15);
break;
case 16:
kf_bfly4(fout, 1, &pfa_st32, 1, 4*15, 4);
kf_bfly4(fout, 2, &pfa_st32, 4, 15, 16);
break;
default:
kf_bfly4(fout, 1, &pfa_st32, 1, 8*15, 4);
kf_bfly2(fout, 4, 4*15);
kf_bfly4(fout, 1, &pfa_st32, 8, 15, 32);
break;
}
}
#endif /* OPUS_PFA */
void opus_fft_c(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
{
int i;

107
lib/rbcodec/codecs/libopus/celt/mdct.c
View file

@ -51,6 +51,41 @@
#include <math.h>
#include "os_support.h"
#include "mathops.h"
#ifdef OPUS_PFA
#include "pfa.h"
#include "pfa_tables.h"
/* The pre-rotation already scatters through a table; for the prime factor
transform it reads this one instead of the bit-reversal, which is why the
gather costs nothing. */
static const opus_int16 *pfa_gather(int n)
{
switch (n)
{
case 60: return pfa_in_60;
case 120: return pfa_in_120;
case 240: return pfa_in_240;
default: return pfa_in_480;
}
}
/* Where the post-rotation finds X[k], as a byte offset. */
static const opus_int16 *pfa_postmap(int n)
{
switch (n)
{
case 60: return pfa_post_60;
case 120: return pfa_post_120;
case 240: return pfa_post_240;
default: return pfa_post_480;
}
}
/* Walk p(k) = M*(k mod 15) + (k mod M) by one step. Both residues are
running counters, so the post-rotation needs no index table. */
#define PFA_ADV_UP(p, k1, k2, M) do { (p) += (M)+1; if (++(k1) == 15) { (k1) = 0; (p) -= 15*(M); } if (++(k2) == (M)) { (k2) = 0; (p) -= (M); } } while (0)
#define PFA_ADV_DN(p, k1, k2, M) do { (p) -= (M)+1; if ((k1)-- == 0) { (k1) = 14; (p) += 15*(M); } if ((k2)-- == 0) { (k2) = (M)-1; (p) += (M); } } while (0)
#endif
#if defined(OPUS_ARM_ASM)
#include "arm/mdct_armv4.h"
#include "arm/mdct_armv5e.h"
@ -249,6 +284,12 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
int i;
int N, N2, N4;
const kiss_twiddle_scalar *trig;
#ifdef OPUS_PFA
int use_pfa;
kiss_fft_cpx *pfa_out = NULL;
VARDECL(kiss_fft_cpx, pfa_tmp);
#endif
SAVE_STACK;
(void) arch;
N = l->n;
@ -261,6 +302,18 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
N2 = N>>1;
N4 = N>>2;
#ifdef OPUS_PFA
/* Pass 1 of the prime factor transform reads fifteen contiguous points
and writes them strided, so it cannot work in place. The input buffer
is exactly the right size and is documented as destroyed here, so for
the usual stride==1 call it costs nothing; a strided call cannot reuse
it, but those are always the short transform, N4 == 60. */
use_pfa = OPUS_PFA_SIZE(N4) && (stride == 1 || N4 == 60);
ALLOC(pfa_tmp, (use_pfa && stride != 1) ? N4 : ALLOC_NONE, kiss_fft_cpx);
if (use_pfa)
pfa_out = (stride == 1) ? (kiss_fft_cpx *)in : pfa_tmp;
#endif
/* Pre-rotate */
{
/* Temp pointers to make it really clear to the compiler what we're doing */
@ -269,6 +322,9 @@ 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 OPUS_PFA
if (use_pfa) bitrev = pfa_gather(N4);
#endif
#ifdef OVERRIDE_MDCT_PREROT
mdct_prerot_opt(xp1, xp2, t, bitrev, yp, N4,
2*stride*(int)sizeof(kiss_fft_scalar));
@ -290,6 +346,11 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
#endif
}
#ifdef OPUS_PFA
if (use_pfa)
opus_pfa_impl((const kiss_fft_cpx*)(out+(overlap>>1)), pfa_out, N4);
else
#endif
opus_fft_impl(l->kfft[shift], (kiss_fft_cpx*)(out+(overlap>>1)));
/* Post-rotate and de-shuffle from both ends of the buffer at once to make
@ -300,6 +361,51 @@ 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 OPUS_PFA
if (use_pfa)
{
/* Same rotation, but the transform is in Good-Thomas order, so the
two ends are gathered through p(k) rather than read in place.
That is also why this writes a different buffer than it reads. */
const kiss_fft_scalar *S = (const kiss_fft_scalar *)pfa_out;
int M = N4/15;
int pa = 0, ka = 0, qa = 0;
int kd = (N4-1)%15, qd = (N4-1)&(M-1);
int pd = M*kd + qd;
#ifdef OVERRIDE_MDCT_POSTROT_PFA
mdct_postrot_pfa_opt(S, yp0, yp1, t, pfa_postmap(N4), N4);
(void)pa; (void)ka; (void)qa; (void)pd; (void)kd; (void)qd; (void)M;
#else
for(i=0;i<(N4+1)>>1;i++)
{
kiss_fft_scalar re, im, yr, yi;
kiss_twiddle_scalar t0, t1;
re = S[2*pa+1];
im = S[2*pa];
t0 = t[i];
t1 = t[N4+i];
yr = ADD32_ovflw(S_MUL(re,t0), S_MUL(im,t1));
yi = SUB32_ovflw(S_MUL(re,t1), S_MUL(im,t0));
re = S[2*pd+1];
im = S[2*pd];
yp0[0] = yr;
yp1[1] = yi;
t0 = t[(N4-i-1)];
t1 = t[(N2-i-1)];
yr = ADD32_ovflw(S_MUL(re,t0), S_MUL(im,t1));
yi = SUB32_ovflw(S_MUL(re,t1), S_MUL(im,t0));
yp1[0] = yr;
yp0[1] = yi;
yp0 += 2;
yp1 -= 2;
PFA_ADV_UP(pa, ka, qa, M);
PFA_ADV_DN(pd, kd, qd, M);
}
#endif
}
else
#endif
#ifdef OVERRIDE_MDCT_POSTROT
mdct_postrot_opt(yp0, yp1, t, N4, (N4+1)>>1);
#else
@ -356,5 +462,6 @@ void clt_mdct_backward_c(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_sca
}
#endif
}
RESTORE_STACK;
}
#endif /* OVERRIDE_clt_mdct_backward */

148
lib/rbcodec/codecs/libopus/celt/pfa.c Normal file
View file

@ -0,0 +1,148 @@
/* Good-Thomas (prime factor) FFT for the CELT backward MDCT.
*
* Copyright (c) 2007-2008 CSIRO
* Copyright (c) 2007-2008 Xiph.Org Foundation
* 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/kiss_fft.c are met.
*
* Why this exists
* ---------------
* Every FFT length a 48 kHz Opus stream uses is 15 times a power of two:
* 60, 120, 240, 480. 15 and 2^k are co-prime, so the prime factor algorithm
* applies, and its whole point is that co-prime factors need no twiddle
* factors between them -- only an index permutation.
*
* That matters because in the mixed-radix decomposition kiss_fft picks
* (5, 3, 4, 2, 4 for N=480) the inter-stage twiddles are 4,020 of the 5,540
* real multiplies, 73% of the total. Splitting 480 as 15 x 32 deletes that
* class of work rather than making it cheaper: 2,780 multiplies, half as
* many. Multiplies are 28.4% of ARMv4 decode, so this is the largest single
* item the profile still offers on that CPU.
*
* Structure, for N = 15*M:
*
* gather x[n] -> u[15*n2 + 3*m2 + m1], a table the pre-rotation reads
* in place of the bit-reversal it read before, so it is free
* pass 1 M independent 15-point DFTs, contiguous in, strided M out
* pass 2 15 independent M-point DFTs, contiguous, which is exactly
* what the existing kf_bfly4/kf_bfly2 kernels already do
* scatter X[k] = B[k mod 15][k mod M], absorbed into the post-rotation
* as two running counters, no table
*
* Buffers. Pass 1 cannot be in place (it reads 15 contiguous and writes 15
* strided), so it needs somewhere to land. clt_mdct_backward already
* destroys its input -- celt_synthesis copies the spectrum out first for
* exactly that reason -- and that buffer is N2 scalars, which is N4 complex,
* exactly the size wanted. So the common stride==1 case costs no memory at
* all. Short blocks arrive with stride>1 and a strided input that cannot be
* reused, but there N4 is only 60, so a 480-byte local covers it.
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "kiss_fft.h"
#include "_kiss_fft_guts.h"
#include "mathops.h"
#include "stack_alloc.h"
#include "pfa.h"
#ifdef OPUS_PFA
/* The same constants the mixed-radix butterflies use, so the arithmetic
below is the arithmetic that was already shipping, only reassociated.
epi3 is -sin(2pi/3); ya/yb are the radix-5 sines; yc is sqrt(5)/4, which
collapses the four radix-5 cosine products into one multiply because
cos(2pi/5)+cos(4pi/5) is exactly -1/2 and the Q15 constants honour that
identity exactly (10126 - 26510 == -16384). */
#define PFA_EPI3 (-28378)
#define PFA_YAI (-31164)
#define PFA_YBI (-19261)
#define PFA_YC ( 18318)
/* PFA_KEEP_C keeps the reference compiled alongside an override so a
harness can score one against the other on target. */
#if !defined(OVERRIDE_PFA_FFT15) || defined(PFA_KEEP_C)
/* One 15-point DFT, itself Good-Thomas over 3 and 5 so it too carries no
twiddles: five 3-point DFTs, then three 5-point DFTs. 40 real multiplies.
Input is five contiguous triples, laid out by the gather table. Output is
scattered by ostride complex elements; within one 5-point DFT the
destinations run at a stride of 3, which is what lets the kernel walk them
with a single post-indexed store. */
void pfa_fft15_c(const kiss_fft_cpx *in, kiss_fft_cpx *out, int ostride)
{
kiss_fft_cpx v[15];
int g, k1;
/* Five 3-point DFTs. Results land strided by 5 so each 5-point DFT
below reads a contiguous run. */
for (g=0;g<5;g++)
{
kiss_fft_cpx x0, s3, s0, xm;
x0 = in[3*g];
s3.r = ADD32_ovflw(in[3*g+1].r, in[3*g+2].r);
s3.i = ADD32_ovflw(in[3*g+1].i, in[3*g+2].i);
s0.r = SUB32_ovflw(in[3*g+1].r, in[3*g+2].r);
s0.i = SUB32_ovflw(in[3*g+1].i, in[3*g+2].i);
xm.r = SUB32_ovflw(x0.r, HALF_OF(s3.r));
xm.i = SUB32_ovflw(x0.i, HALF_OF(s3.i));
s0.r = S_MUL(s0.r, PFA_EPI3);
s0.i = S_MUL(s0.i, PFA_EPI3);
v[g].r = ADD32_ovflw(x0.r, s3.r);
v[g].i = ADD32_ovflw(x0.i, s3.i);
v[5+g].r = SUB32_ovflw(xm.r, s0.i);
v[5+g].i = ADD32_ovflw(xm.i, s0.r);
v[10+g].r = ADD32_ovflw(xm.r, s0.i);
v[10+g].i = SUB32_ovflw(xm.i, s0.r);
}
/* Three 5-point DFTs. Output k2 of DFT k1 belongs at 15-point index
j = (10*k1 + 6*k2) mod 15; walking j from k1 in steps of 3 visits them
in the order k2 = (k1 + 3*t) mod 5, which is the permutation applied
when storing. */
for (k1=0;k1<3;k1++)
{
const kiss_fft_cpx *x = v + 5*k1;
kiss_fft_cpx X[5], s7, s8, s9, s10, s3, s4, s5, s11, s6, s12;
int t;
s7.r = ADD32_ovflw(x[1].r, x[4].r); s7.i = ADD32_ovflw(x[1].i, x[4].i);
s10.r = SUB32_ovflw(x[1].r, x[4].r); s10.i = SUB32_ovflw(x[1].i, x[4].i);
s8.r = ADD32_ovflw(x[2].r, x[3].r); s8.i = ADD32_ovflw(x[2].i, x[3].i);
s9.r = SUB32_ovflw(x[2].r, x[3].r); s9.i = SUB32_ovflw(x[2].i, x[3].i);
s3.r = ADD32_ovflw(s7.r, s8.r); s3.i = ADD32_ovflw(s7.i, s8.i);
s4.r = SUB32_ovflw(s7.r, s8.r); s4.i = SUB32_ovflw(s7.i, s8.i);
X[0].r = ADD32_ovflw(x[0].r, s3.r);
X[0].i = ADD32_ovflw(x[0].i, s3.i);
s3.r = SUB32_ovflw(x[0].r, QUARTER_OF(s3.r));
s3.i = SUB32_ovflw(x[0].i, QUARTER_OF(s3.i));
s4.r = S_MUL(s4.r, PFA_YC);
s4.i = S_MUL(s4.i, PFA_YC);
s5.r = ADD32_ovflw(s3.r, s4.r); s5.i = ADD32_ovflw(s3.i, s4.i);
s11.r = SUB32_ovflw(s3.r, s4.r); s11.i = SUB32_ovflw(s3.i, s4.i);
s6.r = ADD32_ovflw(S_MUL(s10.i, PFA_YAI), S_MUL(s9.i, PFA_YBI));
s6.i = NEG32_ovflw(ADD32_ovflw(S_MUL(s10.r, PFA_YAI), S_MUL(s9.r, PFA_YBI)));
X[1].r = SUB32_ovflw(s5.r, s6.r); X[1].i = SUB32_ovflw(s5.i, s6.i);
X[4].r = ADD32_ovflw(s5.r, s6.r); X[4].i = ADD32_ovflw(s5.i, s6.i);
s12.r = SUB32_ovflw(S_MUL(s9.i, PFA_YAI), S_MUL(s10.i, PFA_YBI));
s12.i = SUB32_ovflw(S_MUL(s10.r, PFA_YBI), S_MUL(s9.r, PFA_YAI));
X[2].r = ADD32_ovflw(s11.r, s12.r); X[2].i = ADD32_ovflw(s11.i, s12.i);
X[3].r = SUB32_ovflw(s11.r, s12.r); X[3].i = SUB32_ovflw(s11.i, s12.i);
for (t=0;t<5;t++)
out[(k1 + 3*t)*ostride] = X[(k1 + 3*t) % 5];
}
}
#endif /* !OVERRIDE_PFA_FFT15 || PFA_KEEP_C */
#endif /* OPUS_PFA */

43
lib/rbcodec/codecs/libopus/celt/pfa.h Normal file
View file

@ -0,0 +1,43 @@
/* Good-Thomas (prime factor) FFT for the CELT backward MDCT.
See pfa.c for what it is and why.
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/kiss_fft.c are met.
*/
#ifndef PFA_H
#define PFA_H
#include "kiss_fft.h"
#ifdef OPUS_PFA
void pfa_fft15_c(const kiss_fft_cpx *in, kiss_fft_cpx *out, int ostride);
#if defined(OPUS_ARM_ASM) && !defined(OPUS_ARM_NO_PFA_ASM)
# include "arm/pfa_arm.h"
#endif
#ifndef OVERRIDE_PFA_FFT15
# define PFA_FFT15(in, out, ostride) pfa_fft15_c(in, out, ostride)
#endif
/* Pass 1 and pass 2. fin holds the gathered spectrum and fout receives the
transform; they must be different buffers, because pass 1 reads fifteen
contiguous points and writes them strided. nfft is 60, 120, 240 or 480.
The result is left in Good-Thomas order: X[k] sits at
fout[(nfft/15)*(k % 15) + (k % (nfft/15))]. The post-rotation reads it
through that index, which costs two running counters and no table. */
void opus_pfa_impl(const kiss_fft_cpx *fin, kiss_fft_cpx *fout, int nfft);
/* True for the lengths opus_pfa_impl handles. Every 48 kHz CELT transform
qualifies; the test exists so a custom mode falls back to kiss_fft. */
#define OPUS_PFA_SIZE(n) \
((n) == 60 || (n) == 120 || (n) == 240 || (n) == 480)
#endif /* OPUS_PFA */
#endif /* PFA_H */

281
lib/rbcodec/codecs/libopus/celt/pfa_tables.h Normal file
View file

@ -0,0 +1,281 @@
/* 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/kiss_fft.c are met.
*/
/* Good-Thomas index tables for the CELT backward MDCT.
Generated by an offline script. Do not edit.
N = 15*M with M a power of two, so 15 and M are co-prime and the prime
factor algorithm applies: the 15-point and M-point passes need no twiddle
factors between them at all, which is where the multiply saving comes
from. All that is left is the index algebra, and all of it is either a
table the pre-rotation was already reading or a running counter.
*/
#ifndef PFA_TABLES_H
#define PFA_TABLES_H
#include "opus_types.h"
#include "arch.h"
/* n -> 15*n2 + 3*m2 + m1, for N=60, M=4. */
static const opus_int16 pfa_in_60[60] ICONST_ATTR = {
0, 56, 34, 27, 8, 46, 39, 20, 13, 51, 32, 25,
3, 59, 37, 15, 11, 49, 42, 23, 1, 54, 35, 28,
6, 47, 40, 18, 14, 52, 30, 26, 4, 57, 38, 16,
9, 50, 43, 21, 2, 55, 33, 29, 7, 45, 41, 19,
12, 53, 31, 24, 5, 58, 36, 17, 10, 48, 44, 22,
};
/* n -> 15*n2 + 3*m2 + m1, for N=120, M=8. */
static const opus_int16 pfa_in_120[120] ICONST_ATTR = {
0, 118, 101, 81, 64, 47, 42, 25, 8, 108, 91, 89,
69, 52, 35, 15, 13, 116, 96, 79, 62, 57, 40, 23,
3, 106, 104, 84, 67, 50, 30, 28, 11, 111, 94, 77,
72, 55, 38, 18, 1, 119, 99, 82, 65, 45, 43, 26,
6, 109, 92, 87, 70, 53, 33, 16, 14, 114, 97, 80,
60, 58, 41, 21, 4, 107, 102, 85, 68, 48, 31, 29,
9, 112, 95, 75, 73, 56, 36, 19, 2, 117, 100, 83,
63, 46, 44, 24, 7, 110, 90, 88, 71, 51, 34, 17,
12, 115, 98, 78, 61, 59, 39, 22, 5, 105, 103, 86,
66, 49, 32, 27, 10, 113, 93, 76, 74, 54, 37, 20,
};
/* n -> 15*n2 + 3*m2 + m1, for N=240, M=16. */
static const opus_int16 pfa_in_240[240] ICONST_ATTR = {
0, 233, 223, 198, 191, 166, 156, 149, 124, 114, 92, 82,
72, 50, 40, 15, 8, 238, 213, 206, 181, 171, 164, 139,
129, 107, 97, 87, 65, 55, 30, 23, 13, 228, 221, 196,
186, 179, 154, 144, 122, 112, 102, 80, 70, 45, 38, 28,
3, 236, 211, 201, 194, 169, 159, 137, 127, 117, 95, 85,
60, 53, 43, 18, 11, 226, 216, 209, 184, 174, 152, 142,
132, 110, 100, 75, 68, 58, 33, 26, 1, 231, 224, 199,
189, 167, 157, 147, 125, 115, 90, 83, 73, 48, 41, 16,
6, 239, 214, 204, 182, 172, 162, 140, 130, 105, 98, 88,
63, 56, 31, 21, 14, 229, 219, 197, 187, 177, 155, 145,
120, 113, 103, 78, 71, 46, 36, 29, 4, 234, 212, 202,
192, 170, 160, 135, 128, 118, 93, 86, 61, 51, 44, 19,
9, 227, 217, 207, 185, 175, 150, 143, 133, 108, 101, 76,
66, 59, 34, 24, 2, 232, 222, 200, 190, 165, 158, 148,
123, 116, 91, 81, 74, 49, 39, 17, 7, 237, 215, 205,
180, 173, 163, 138, 131, 106, 96, 89, 64, 54, 32, 22,
12, 230, 220, 195, 188, 178, 153, 146, 121, 111, 104, 79,
69, 47, 37, 27, 5, 235, 210, 203, 193, 168, 161, 136,
126, 119, 94, 84, 62, 52, 42, 20, 10, 225, 218, 208,
183, 176, 151, 141, 134, 109, 99, 77, 67, 57, 35, 25,
};
/* n -> 15*n2 + 3*m2 + m1, for N=480, M=32. */
static const opus_int16 pfa_in_480[480] ICONST_ATTR = {
0, 229, 458, 204, 433, 167, 393, 142, 371, 117, 331, 80,
306, 55, 284, 15, 244, 473, 219, 448, 182, 408, 157, 386,
132, 346, 95, 321, 70, 299, 30, 259, 8, 234, 463, 197,
423, 172, 401, 147, 361, 110, 336, 85, 314, 45, 274, 23,
249, 478, 212, 438, 187, 416, 162, 376, 125, 351, 100, 329,
60, 289, 38, 264, 13, 227, 453, 202, 431, 177, 391, 140,
366, 115, 344, 75, 304, 53, 279, 28, 242, 468, 217, 446,
192, 406, 155, 381, 130, 359, 90, 319, 68, 294, 43, 257,
3, 232, 461, 207, 421, 170, 396, 145, 374, 105, 334, 83,
309, 58, 272, 18, 247, 476, 222, 436, 185, 411, 160, 389,
120, 349, 98, 324, 73, 287, 33, 262, 11, 237, 451, 200,
426, 175, 404, 135, 364, 113, 339, 88, 302, 48, 277, 26,
252, 466, 215, 441, 190, 419, 150, 379, 128, 354, 103, 317,
63, 292, 41, 267, 1, 230, 456, 205, 434, 165, 394, 143,
369, 118, 332, 78, 307, 56, 282, 16, 245, 471, 220, 449,
180, 409, 158, 384, 133, 347, 93, 322, 71, 297, 31, 260,
6, 235, 464, 195, 424, 173, 399, 148, 362, 108, 337, 86,
312, 46, 275, 21, 250, 479, 210, 439, 188, 414, 163, 377,
123, 352, 101, 327, 61, 290, 36, 265, 14, 225, 454, 203,
429, 178, 392, 138, 367, 116, 342, 76, 305, 51, 280, 29,
240, 469, 218, 444, 193, 407, 153, 382, 131, 357, 91, 320,
66, 295, 44, 255, 4, 233, 459, 208, 422, 168, 397, 146,
372, 106, 335, 81, 310, 59, 270, 19, 248, 474, 223, 437,
183, 412, 161, 387, 121, 350, 96, 325, 74, 285, 34, 263,
9, 238, 452, 198, 427, 176, 402, 136, 365, 111, 340, 89,
300, 49, 278, 24, 253, 467, 213, 442, 191, 417, 151, 380,
126, 355, 104, 315, 64, 293, 39, 268, 2, 228, 457, 206,
432, 166, 395, 141, 370, 119, 330, 79, 308, 54, 283, 17,
243, 472, 221, 447, 181, 410, 156, 385, 134, 345, 94, 323,
69, 298, 32, 258, 7, 236, 462, 196, 425, 171, 400, 149,
360, 109, 338, 84, 313, 47, 273, 22, 251, 477, 211, 440,
186, 415, 164, 375, 124, 353, 99, 328, 62, 288, 37, 266,
12, 226, 455, 201, 430, 179, 390, 139, 368, 114, 343, 77,
303, 52, 281, 27, 241, 470, 216, 445, 194, 405, 154, 383,
129, 358, 92, 318, 67, 296, 42, 256, 5, 231, 460, 209,
420, 169, 398, 144, 373, 107, 333, 82, 311, 57, 271, 20,
246, 475, 224, 435, 184, 413, 159, 388, 122, 348, 97, 326,
72, 286, 35, 261, 10, 239, 450, 199, 428, 174, 403, 137,
363, 112, 341, 87, 301, 50, 276, 25, 254, 465, 214, 443,
189, 418, 152, 378, 127, 356, 102, 316, 65, 291, 40, 269,
};
/* Bit-reversal of the M-point sub-FFT, folded into the 15-point
kernel's output scatter so the radix-4 chain sees natural order. */
static const opus_int16 pfa_brev_4[4] ICONST_ATTR = {
0, 1, 2, 3,
};
static const opus_int16 pfa_brev_8[8] ICONST_ATTR = {
0, 4, 1, 5, 2, 6, 3, 7,
};
static const opus_int16 pfa_brev_16[16] ICONST_ATTR = {
0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15,
};
static const opus_int16 pfa_brev_32[32] ICONST_ATTR = {
0, 8, 16, 24, 4, 12, 20, 28, 1, 9, 17, 25, 5, 13, 21, 29,
2, 10, 18, 26, 6, 14, 22, 30, 3, 11, 19, 27, 7, 15, 23, 31,
};
/* Where the post-rotation finds X[k]. The transform is left in
Good-Thomas order, X[k] at M*(k mod 15) + (k mod M). Both residues are
running counters, but keeping four of them plus their wrap constants
costs more registers than the kernel has to spare, so the composed index
is tabulated instead, pre-scaled to a byte offset. The post-rotation
walks it forwards for the ascending end and backwards for the descending
one, so a single table serves both. */
static const opus_int16 pfa_post_60[60] ICONST_ATTR = {
0, 40, 80, 120, 128, 168, 208, 248, 256, 296,
336, 376, 384, 424, 464, 24, 32, 72, 112, 152,
160, 200, 240, 280, 288, 328, 368, 408, 416, 456,
16, 56, 64, 104, 144, 184, 192, 232, 272, 312,
320, 360, 400, 440, 448, 8, 48, 88, 96, 136,
176, 216, 224, 264, 304, 344, 352, 392, 432, 472,
};
static const opus_int16 pfa_post_120[120] ICONST_ATTR = {
0, 72, 144, 216, 288, 360, 432, 504, 512, 584,
656, 728, 800, 872, 944, 56, 64, 136, 208, 280,
352, 424, 496, 568, 576, 648, 720, 792, 864, 936,
48, 120, 128, 200, 272, 344, 416, 488, 560, 632,
640, 712, 784, 856, 928, 40, 112, 184, 192, 264,
336, 408, 480, 552, 624, 696, 704, 776, 848, 920,
32, 104, 176, 248, 256, 328, 400, 472, 544, 616,
688, 760, 768, 840, 912, 24, 96, 168, 240, 312,
320, 392, 464, 536, 608, 680, 752, 824, 832, 904,
16, 88, 160, 232, 304, 376, 384, 456, 528, 600,
672, 744, 816, 888, 896, 8, 80, 152, 224, 296,
368, 440, 448, 520, 592, 664, 736, 808, 880, 952,
};
static const opus_int16 pfa_post_240[240] ICONST_ATTR = {
0, 136, 272, 408, 544, 680, 816, 952, 1088, 1224,
1360, 1496, 1632, 1768, 1904, 120, 128, 264, 400, 536,
672, 808, 944, 1080, 1216, 1352, 1488, 1624, 1760, 1896,
112, 248, 256, 392, 528, 664, 800, 936, 1072, 1208,
1344, 1480, 1616, 1752, 1888, 104, 240, 376, 384, 520,
656, 792, 928, 1064, 1200, 1336, 1472, 1608, 1744, 1880,
96, 232, 368, 504, 512, 648, 784, 920, 1056, 1192,
1328, 1464, 1600, 1736, 1872, 88, 224, 360, 496, 632,
640, 776, 912, 1048, 1184, 1320, 1456, 1592, 1728, 1864,
80, 216, 352, 488, 624, 760, 768, 904, 1040, 1176,
1312, 1448, 1584, 1720, 1856, 72, 208, 344, 480, 616,
752, 888, 896, 1032, 1168, 1304, 1440, 1576, 1712, 1848,
64, 200, 336, 472, 608, 744, 880, 1016, 1024, 1160,
1296, 1432, 1568, 1704, 1840, 56, 192, 328, 464, 600,
736, 872, 1008, 1144, 1152, 1288, 1424, 1560, 1696, 1832,
48, 184, 320, 456, 592, 728, 864, 1000, 1136, 1272,
1280, 1416, 1552, 1688, 1824, 40, 176, 312, 448, 584,
720, 856, 992, 1128, 1264, 1400, 1408, 1544, 1680, 1816,
32, 168, 304, 440, 576, 712, 848, 984, 1120, 1256,
1392, 1528, 1536, 1672, 1808, 24, 160, 296, 432, 568,
704, 840, 976, 1112, 1248, 1384, 1520, 1656, 1664, 1800,
16, 152, 288, 424, 560, 696, 832, 968, 1104, 1240,
1376, 1512, 1648, 1784, 1792, 8, 144, 280, 416, 552,
688, 824, 960, 1096, 1232, 1368, 1504, 1640, 1776, 1912,
};
static const opus_int16 pfa_post_480[480] ICONST_ATTR = {
0, 264, 528, 792, 1056, 1320, 1584, 1848, 2112, 2376,
2640, 2904, 3168, 3432, 3696, 120, 384, 648, 912, 1176,
1440, 1704, 1968, 2232, 2496, 2760, 3024, 3288, 3552, 3816,
240, 504, 512, 776, 1040, 1304, 1568, 1832, 2096, 2360,
2624, 2888, 3152, 3416, 3680, 104, 368, 632, 896, 1160,
1424, 1688, 1952, 2216, 2480, 2744, 3008, 3272, 3536, 3800,
224, 488, 752, 1016, 1024, 1288, 1552, 1816, 2080, 2344,
2608, 2872, 3136, 3400, 3664, 88, 352, 616, 880, 1144,
1408, 1672, 1936, 2200, 2464, 2728, 2992, 3256, 3520, 3784,
208, 472, 736, 1000, 1264, 1528, 1536, 1800, 2064, 2328,
2592, 2856, 3120, 3384, 3648, 72, 336, 600, 864, 1128,
1392, 1656, 1920, 2184, 2448, 2712, 2976, 3240, 3504, 3768,
192, 456, 720, 984, 1248, 1512, 1776, 2040, 2048, 2312,
2576, 2840, 3104, 3368, 3632, 56, 320, 584, 848, 1112,
1376, 1640, 1904, 2168, 2432, 2696, 2960, 3224, 3488, 3752,
176, 440, 704, 968, 1232, 1496, 1760, 2024, 2288, 2552,
2560, 2824, 3088, 3352, 3616, 40, 304, 568, 832, 1096,
1360, 1624, 1888, 2152, 2416, 2680, 2944, 3208, 3472, 3736,
160, 424, 688, 952, 1216, 1480, 1744, 2008, 2272, 2536,
2800, 3064, 3072, 3336, 3600, 24, 288, 552, 816, 1080,
1344, 1608, 1872, 2136, 2400, 2664, 2928, 3192, 3456, 3720,
144, 408, 672, 936, 1200, 1464, 1728, 1992, 2256, 2520,
2784, 3048, 3312, 3576, 3584, 8, 272, 536, 800, 1064,
1328, 1592, 1856, 2120, 2384, 2648, 2912, 3176, 3440, 3704,
128, 392, 656, 920, 1184, 1448, 1712, 1976, 2240, 2504,
2768, 3032, 3296, 3560, 3824, 248, 256, 520, 784, 1048,
1312, 1576, 1840, 2104, 2368, 2632, 2896, 3160, 3424, 3688,
112, 376, 640, 904, 1168, 1432, 1696, 1960, 2224, 2488,
2752, 3016, 3280, 3544, 3808, 232, 496, 760, 768, 1032,
1296, 1560, 1824, 2088, 2352, 2616, 2880, 3144, 3408, 3672,
96, 360, 624, 888, 1152, 1416, 1680, 1944, 2208, 2472,
2736, 3000, 3264, 3528, 3792, 216, 480, 744, 1008, 1272,
1280, 1544, 1808, 2072, 2336, 2600, 2864, 3128, 3392, 3656,
80, 344, 608, 872, 1136, 1400, 1664, 1928, 2192, 2456,
2720, 2984, 3248, 3512, 3776, 200, 464, 728, 992, 1256,
1520, 1784, 1792, 2056, 2320, 2584, 2848, 3112, 3376, 3640,
64, 328, 592, 856, 1120, 1384, 1648, 1912, 2176, 2440,
2704, 2968, 3232, 3496, 3760, 184, 448, 712, 976, 1240,
1504, 1768, 2032, 2296, 2304, 2568, 2832, 3096, 3360, 3624,
48, 312, 576, 840, 1104, 1368, 1632, 1896, 2160, 2424,
2688, 2952, 3216, 3480, 3744, 168, 432, 696, 960, 1224,
1488, 1752, 2016, 2280, 2544, 2808, 2816, 3080, 3344, 3608,
32, 296, 560, 824, 1088, 1352, 1616, 1880, 2144, 2408,
2672, 2936, 3200, 3464, 3728, 152, 416, 680, 944, 1208,
1472, 1736, 2000, 2264, 2528, 2792, 3056, 3320, 3328, 3592,
16, 280, 544, 808, 1072, 1336, 1600, 1864, 2128, 2392,
2656, 2920, 3184, 3448, 3712, 136, 400, 664, 928, 1192,
1456, 1720, 1984, 2248, 2512, 2776, 3040, 3304, 3568, 3832,
};
/* Twiddles for the M-point half. One 32-entry table serves every size:
M=16 reads every second entry, M=8 every fourth, M=4 every eighth, which
is exactly the shift mechanism kiss_fft already uses to share one table
between the four transform lengths. */
static const kiss_twiddle_cpx pfa_tw32[32] ICONST_ATTR = {
{ 32767, 0},
{ 32137, -6393},
{ 30273, -12539},
{ 27245, -18204},
{ 23170, -23170},
{ 18204, -27245},
{ 12539, -30273},
{ 6393, -32137},
{ 0, -32767},
{ -6393, -32137},
{-12539, -30273},
{-18204, -27245},
{-23170, -23170},
{-27245, -18204},
{-30273, -12539},
{-32137, -6393},
{-32767, 0},
{-32137, 6393},
{-30273, 12539},
{-27245, 18204},
{-23170, 23170},
{-18204, 27245},
{-12539, 30273},
{ -6393, 32137},
{ 0, 32767},
{ 6393, 32137},
{ 12539, 30273},
{ 18204, 27245},
{ 23170, 23170},
{ 27245, 18204},
{ 30273, 12539},
{ 32137, 6393},
};
#endif /* PFA_TABLES_H */

8
lib/rbcodec/codecs/libopus/config.h
View file

@ -82,6 +82,14 @@
#endif
/* Good-Thomas FFT for the backward MDCT. Every 48 kHz CELT transform length
is 15 times a power of two, so the prime factor algorithm applies and the
inter-stage twiddles -- 73% of the FFT multiplies -- disappear. Define
OPUS_NO_PFA to fall back to the mixed-radix chain. */
#ifndef OPUS_NO_PFA
#define OPUS_PFA
#endif
#if defined(CPU_COLDFIRE)
#define OPUS_CF_INLINE_ASM
#endif