ref: 38e535a1e7177c4cede1b11cc79089d900d04c81
parent: 3116093d34f794238acdc8c509dea6ace3aa1b95
author: Jean-Marc Valin <jeanmarcv@google.com>
date: Thu Jul 4 21:28:30 EDT 2024
Improve MDCT scaling to maximize accuracy
--- a/celt/mdct.c
+++ b/celt/mdct.c
@@ -133,6 +133,7 @@
/* Allows us to scale with MULT16_32_Q16(), which is faster than
MULT16_32_Q15() on ARM. */
int scale_shift = st->scale_shift-1;
+ int headroom;
#endif
SAVE_STACK;
(void)arch;
@@ -195,6 +196,9 @@
{
kiss_fft_scalar * OPUS_RESTRICT yp = f;
const kiss_twiddle_scalar *t = &trig[0];
+#ifdef FIXED_POINT
+ opus_val32 maxval=1;
+#endif
for(i=0;i<N4;i++)
{
kiss_fft_cpx yc;
@@ -208,10 +212,21 @@
yi = S_MUL(im,t0) + S_MUL(re,t1);
yc.r = yr;
yc.i = yi;
- yc.r = PSHR32(S_MUL2(yc.r, scale), scale_shift);
- yc.i = PSHR32(S_MUL2(yc.i, scale), scale_shift);
+ yc.r = S_MUL2(yc.r, scale);
+ yc.i = S_MUL2(yc.i, scale);
+#ifdef FIXED_POINT
+ maxval = MAX32(maxval, MAX32(yc.r, yc.i));
+#endif
f2[st->bitrev[i]] = yc;
}
+#ifdef FIXED_POINT
+ headroom = IMAX(0, IMIN(scale_shift-1, 28-celt_ilog2(maxval)));
+ for(i=0;i<N4;i++)
+ {
+ f2[i].r = PSHR32(f2[i].r, scale_shift-headroom);
+ f2[i].i = PSHR32(f2[i].i, scale_shift-headroom);
+ }
+#endif
}
/* N/4 complex FFT, does not downscale anymore */
@@ -228,8 +243,8 @@
for(i=0;i<N4;i++)
{
kiss_fft_scalar yr, yi;
- yr = S_MUL(fp->i,t[N4+i]) - S_MUL(fp->r,t[i]);
- yi = S_MUL(fp->r,t[N4+i]) + S_MUL(fp->i,t[i]);
+ yr = PSHR32(S_MUL(fp->i,t[N4+i]) - S_MUL(fp->r,t[i]), headroom);
+ yi = PSHR32(S_MUL(fp->r,t[N4+i]) + S_MUL(fp->i,t[i]), headroom);
*yp1 = yr;
*yp2 = yi;
fp++;
@@ -272,9 +287,12 @@
{
int rev;
kiss_fft_scalar yr, yi;
+ opus_val32 x1, x2;
rev = *bitrev++;
- yr = ADD32_ovflw(S_MUL(*xp2, t[i]), S_MUL(*xp1, t[N4+i]));
- yi = SUB32_ovflw(S_MUL(*xp1, t[i]), S_MUL(*xp2, t[N4+i]));
+ x1 = SHL32(*xp1, IMDCT_HEADROOM);
+ x2 = SHL32(*xp2, IMDCT_HEADROOM);
+ yr = ADD32_ovflw(S_MUL(x2, t[i]), S_MUL(x1, t[N4+i]));
+ yi = SUB32_ovflw(S_MUL(x1, t[i]), S_MUL(x2, t[N4+i]));
/* We swap real and imag because we use an FFT instead of an IFFT. */
yp[2*rev+1] = yr;
yp[2*rev] = yi;
@@ -304,8 +322,8 @@
t0 = t[i];
t1 = t[N4+i];
/* We'd scale up by 2 here, but instead it's done when mixing the windows */
- yr = ADD32_ovflw(S_MUL(re,t0), S_MUL(im,t1));
- yi = SUB32_ovflw(S_MUL(re,t1), S_MUL(im,t0));
+ yr = PSHR32(ADD32_ovflw(S_MUL(re,t0), S_MUL(im,t1)), IMDCT_HEADROOM);
+ yi = PSHR32(SUB32_ovflw(S_MUL(re,t1), S_MUL(im,t0)), IMDCT_HEADROOM);
/* We swap real and imag because we're using an FFT instead of an IFFT. */
re = yp1[1];
im = yp1[0];
@@ -315,8 +333,8 @@
t0 = t[(N4-i-1)];
t1 = t[(N2-i-1)];
/* We'd scale up by 2 here, but instead it's done when mixing the windows */
- yr = ADD32_ovflw(S_MUL(re,t0), S_MUL(im,t1));
- yi = SUB32_ovflw(S_MUL(re,t1), S_MUL(im,t0));
+ yr = PSHR32(ADD32_ovflw(S_MUL(re,t0), S_MUL(im,t1)), IMDCT_HEADROOM);
+ yi = PSHR32(SUB32_ovflw(S_MUL(re,t1), S_MUL(im,t0)), IMDCT_HEADROOM);
yp1[0] = yr;
yp0[1] = yi;
yp0 += 2;
--- a/celt/mdct.h
+++ b/celt/mdct.h
@@ -57,6 +57,9 @@
#include "arm/mdct_arm.h"
#endif
+/* There should be 2 bits of headroom in the IMDCT which we can take
+ advantage of to maximize accuracy. */
+#define IMDCT_HEADROOM 2
int clt_mdct_init(mdct_lookup *l,int N, int maxshift, int arch);
void clt_mdct_clear(mdct_lookup *l, int arch);
--
⑨