Fixes an aliasing bug in the MDCT when the frame size isn't a multiple of 4.
[opus.git] / celt / mdct.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2008 Xiph.Org Foundation
3    Written by Jean-Marc Valin */
4 /*
5    Redistribution and use in source and binary forms, with or without
6    modification, are permitted provided that the following conditions
7    are met:
8
9    - Redistributions of source code must retain the above copyright
10    notice, this list of conditions and the following disclaimer.
11
12    - Redistributions in binary form must reproduce the above copyright
13    notice, this list of conditions and the following disclaimer in the
14    documentation and/or other materials provided with the distribution.
15
16    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
20    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
21    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
22    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 /* This is a simple MDCT implementation that uses a N/4 complex FFT
30    to do most of the work. It should be relatively straightforward to
31    plug in pretty much and FFT here.
32
33    This replaces the Vorbis FFT (and uses the exact same API), which
34    was a bit too messy and that was ending up duplicating code
35    (might as well use the same FFT everywhere).
36
37    The algorithm is similar to (and inspired from) Fabrice Bellard's
38    MDCT implementation in FFMPEG, but has differences in signs, ordering
39    and scaling in many places.
40 */
41
42 #ifndef SKIP_CONFIG_H
43 #ifdef HAVE_CONFIG_H
44 #include "config.h"
45 #endif
46 #endif
47
48 #include "mdct.h"
49 #include "kiss_fft.h"
50 #include "_kiss_fft_guts.h"
51 #include <math.h>
52 #include "os_support.h"
53 #include "mathops.h"
54 #include "stack_alloc.h"
55
56 #ifdef CUSTOM_MODES
57
58 int clt_mdct_init(mdct_lookup *l,int N, int maxshift)
59 {
60    int i;
61    kiss_twiddle_scalar *trig;
62    int shift;
63    int N2=N>>1;
64    l->n = N;
65    l->maxshift = maxshift;
66    for (i=0;i<=maxshift;i++)
67    {
68       if (i==0)
69          l->kfft[i] = opus_fft_alloc(N>>2>>i, 0, 0);
70       else
71          l->kfft[i] = opus_fft_alloc_twiddles(N>>2>>i, 0, 0, l->kfft[0]);
72 #ifndef ENABLE_TI_DSPLIB55
73       if (l->kfft[i]==NULL)
74          return 0;
75 #endif
76    }
77    l->trig = trig = (kiss_twiddle_scalar*)opus_alloc((N-(N2>>maxshift))*sizeof(kiss_twiddle_scalar));
78    if (l->trig==NULL)
79      return 0;
80    for (shift=0;shift<=maxshift;shift++)
81    {
82       /* We have enough points that sine isn't necessary */
83 #if defined(FIXED_POINT)
84 #if 1
85       for (i=0;i<N2;i++)
86          trig[i] = TRIG_UPSCALE*celt_cos_norm(DIV32(ADD32(SHL32(EXTEND32(i),17),N2+16384),N));
87 #else
88       for (i=0;i<N2;i++)
89          trig[i] = (kiss_twiddle_scalar)MAX32(-32767,MIN32(32767,floor(.5+32768*cos(2*M_PI*(i+.125)/N))));
90 #endif
91 #else
92       for (i=0;i<N2;i++)
93          trig[i] = (kiss_twiddle_scalar)cos(2*PI*(i+.125)/N);
94 #endif
95       trig += N2;
96       N2 >>= 1;
97       N >>= 1;
98    }
99    return 1;
100 }
101
102 void clt_mdct_clear(mdct_lookup *l)
103 {
104    int i;
105    for (i=0;i<=l->maxshift;i++)
106       opus_fft_free(l->kfft[i]);
107    opus_free((kiss_twiddle_scalar*)l->trig);
108 }
109
110 #endif /* CUSTOM_MODES */
111
112 /* Forward MDCT trashes the input array */
113 void clt_mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * OPUS_RESTRICT out,
114       const opus_val16 *window, int overlap, int shift, int stride)
115 {
116    int i;
117    int N, N2, N4;
118    VARDECL(kiss_fft_scalar, f);
119    VARDECL(kiss_fft_cpx, f2);
120    const kiss_fft_state *st = l->kfft[shift];
121    const kiss_twiddle_scalar *trig;
122    opus_val16 scale;
123 #ifdef FIXED_POINT
124    /* Allows us to scale with MULT16_32_Q16(), which is faster than
125       MULT16_32_Q15() on ARM. */
126    int scale_shift = st->scale_shift-1;
127 #endif
128    SAVE_STACK;
129    scale = st->scale;
130
131    N = l->n;
132    trig = l->trig;
133    for (i=0;i<shift;i++)
134    {
135       N >>= 1;
136       trig += N;
137    }
138    N2 = N>>1;
139    N4 = N>>2;
140
141    ALLOC(f, N2, kiss_fft_scalar);
142    ALLOC(f2, N4, kiss_fft_cpx);
143
144    /* Consider the input to be composed of four blocks: [a, b, c, d] */
145    /* Window, shuffle, fold */
146    {
147       /* Temp pointers to make it really clear to the compiler what we're doing */
148       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in+(overlap>>1);
149       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+N2-1+(overlap>>1);
150       kiss_fft_scalar * OPUS_RESTRICT yp = f;
151       const opus_val16 * OPUS_RESTRICT wp1 = window+(overlap>>1);
152       const opus_val16 * OPUS_RESTRICT wp2 = window+(overlap>>1)-1;
153       for(i=0;i<((overlap+3)>>2);i++)
154       {
155          /* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
156          *yp++ = MULT16_32_Q15(*wp2, xp1[N2]) + MULT16_32_Q15(*wp1,*xp2);
157          *yp++ = MULT16_32_Q15(*wp1, *xp1)    - MULT16_32_Q15(*wp2, xp2[-N2]);
158          xp1+=2;
159          xp2-=2;
160          wp1+=2;
161          wp2-=2;
162       }
163       wp1 = window;
164       wp2 = window+overlap-1;
165       for(;i<N4-((overlap+3)>>2);i++)
166       {
167          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
168          *yp++ = *xp2;
169          *yp++ = *xp1;
170          xp1+=2;
171          xp2-=2;
172       }
173       for(;i<N4;i++)
174       {
175          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
176          *yp++ =  -MULT16_32_Q15(*wp1, xp1[-N2]) + MULT16_32_Q15(*wp2, *xp2);
177          *yp++ = MULT16_32_Q15(*wp2, *xp1)     + MULT16_32_Q15(*wp1, xp2[N2]);
178          xp1+=2;
179          xp2-=2;
180          wp1+=2;
181          wp2-=2;
182       }
183    }
184    /* Pre-rotation */
185    {
186       kiss_fft_scalar * OPUS_RESTRICT yp = f;
187       const kiss_twiddle_scalar *t = &trig[0];
188       for(i=0;i<N4;i++)
189       {
190          kiss_fft_cpx yc;
191          kiss_twiddle_scalar t0, t1;
192          kiss_fft_scalar re, im, yr, yi;
193          t0 = t[i];
194          t1 = t[N4+i];
195          re = *yp++;
196          im = *yp++;
197          yr = S_MUL(re,t0)  -  S_MUL(im,t1);
198          yi = S_MUL(im,t0)  +  S_MUL(re,t1);
199          yc.r = yr;
200          yc.i = yi;
201          yc.r = PSHR32(MULT16_32_Q16(scale, yc.r), scale_shift);
202          yc.i = PSHR32(MULT16_32_Q16(scale, yc.i), scale_shift);
203          f2[st->bitrev[i]] = yc;
204       }
205    }
206
207    /* N/4 complex FFT, does not downscale anymore */
208    opus_fft_impl(st, f2);
209
210    /* Post-rotate */
211    {
212       /* Temp pointers to make it really clear to the compiler what we're doing */
213       const kiss_fft_cpx * OPUS_RESTRICT fp = f2;
214       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
215       kiss_fft_scalar * OPUS_RESTRICT yp2 = out+stride*(N2-1);
216       const kiss_twiddle_scalar *t = &trig[0];
217       /* Temp pointers to make it really clear to the compiler what we're doing */
218       for(i=0;i<N4;i++)
219       {
220          kiss_fft_scalar yr, yi;
221          yr = S_MUL(fp->i,t[N4+i]) - S_MUL(fp->r,t[i]);
222          yi = S_MUL(fp->r,t[N4+i]) + S_MUL(fp->i,t[i]);
223          *yp1 = yr;
224          *yp2 = yi;
225          fp++;
226          yp1 += 2*stride;
227          yp2 -= 2*stride;
228       }
229    }
230    RESTORE_STACK;
231 }
232
233 void clt_mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * OPUS_RESTRICT out,
234       const opus_val16 * OPUS_RESTRICT window, int overlap, int shift, int stride)
235 {
236    int i;
237    int N, N2, N4;
238    const kiss_twiddle_scalar *trig;
239
240    N = l->n;
241    trig = l->trig;
242    for (i=0;i<shift;i++)
243    {
244       N >>= 1;
245       trig += N;
246    }
247    N2 = N>>1;
248    N4 = N>>2;
249
250    /* Pre-rotate */
251    {
252       /* Temp pointers to make it really clear to the compiler what we're doing */
253       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in;
254       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+stride*(N2-1);
255       kiss_fft_scalar * OPUS_RESTRICT yp = out+(overlap>>1);
256       const kiss_twiddle_scalar * OPUS_RESTRICT t = &trig[0];
257       const opus_int16 * OPUS_RESTRICT bitrev = l->kfft[shift]->bitrev;
258       for(i=0;i<N4;i++)
259       {
260          int rev;
261          kiss_fft_scalar yr, yi;
262          rev = *bitrev++;
263          yr = S_MUL(*xp2, t[i]) + S_MUL(*xp1, t[N4+i]);
264          yi = S_MUL(*xp1, t[i]) - S_MUL(*xp2, t[N4+i]);
265          /* We swap real and imag because we use an FFT instead of an IFFT. */
266          yp[2*rev+1] = yr;
267          yp[2*rev] = yi;
268          /* Storing the pre-rotation directly in the bitrev order. */
269          xp1+=2*stride;
270          xp2-=2*stride;
271       }
272    }
273
274    opus_fft_impl(l->kfft[shift], (kiss_fft_cpx*)(out+(overlap>>1)));
275
276    /* Post-rotate and de-shuffle from both ends of the buffer at once to make
277       it in-place. */
278    {
279       kiss_fft_scalar * yp0 = out+(overlap>>1);
280       kiss_fft_scalar * yp1 = out+(overlap>>1)+N2-2;
281       const kiss_twiddle_scalar *t = &trig[0];
282       /* Loop to (N4+1)>>1 to handle odd N4. When N4 is odd, the
283          middle pair will be computed twice. */
284       for(i=0;i<(N4+1)>>1;i++)
285       {
286          kiss_fft_scalar re, im, yr, yi;
287          kiss_twiddle_scalar t0, t1;
288          /* We swap real and imag because we're using an FFT instead of an IFFT. */
289          re = yp0[1];
290          im = yp0[0];
291          t0 = t[i];
292          t1 = t[N4+i];
293          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
294          yr = S_MUL(re,t0) + S_MUL(im,t1);
295          yi = S_MUL(re,t1) - S_MUL(im,t0);
296          /* We swap real and imag because we're using an FFT instead of an IFFT. */
297          re = yp1[1];
298          im = yp1[0];
299          yp0[0] = yr;
300          yp1[1] = yi;
301
302          t0 = t[(N4-i-1)];
303          t1 = t[(N2-i-1)];
304          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
305          yr = S_MUL(re,t0) + S_MUL(im,t1);
306          yi = S_MUL(re,t1) - S_MUL(im,t0);
307          yp1[0] = yr;
308          yp0[1] = yi;
309          yp0 += 2;
310          yp1 -= 2;
311       }
312    }
313
314    /* Mirror on both sides for TDAC */
315    {
316       kiss_fft_scalar * OPUS_RESTRICT xp1 = out+overlap-1;
317       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
318       const opus_val16 * OPUS_RESTRICT wp1 = window;
319       const opus_val16 * OPUS_RESTRICT wp2 = window+overlap-1;
320
321       for(i = 0; i < overlap/2; i++)
322       {
323          kiss_fft_scalar x1, x2;
324          x1 = *xp1;
325          x2 = *yp1;
326          *yp1++ = MULT16_32_Q15(*wp2, x2) - MULT16_32_Q15(*wp1, x1);
327          *xp1-- = MULT16_32_Q15(*wp1, x2) + MULT16_32_Q15(*wp2, x1);
328          wp1++;
329          wp2--;
330       }
331    }
332 }