Removes the separate 1/8N rotation in the (I)MDCT and unmerges the MDCT sizes
[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 #ifdef FIXED_POINT
123    /* FIXME: This should eventually just go in the state. */
124    opus_val16 scale;
125    int scale_shift;
126    scale_shift = celt_ilog2(st->nfft);
127    if (st->nfft == 1<<scale_shift)
128       scale = Q15ONE;
129    else
130       scale = (1073741824+st->nfft/2)/st->nfft>>(15-scale_shift);
131 #endif
132    SAVE_STACK;
133
134    N = l->n;
135    trig = l->trig;
136    for (i=0;i<shift;i++)
137    {
138       N >>= 1;
139       trig += N;
140    }
141    N2 = N>>1;
142    N4 = N>>2;
143
144    ALLOC(f, N2, kiss_fft_scalar);
145    ALLOC(f2, N2, kiss_fft_cpx);
146
147    /* Consider the input to be composed of four blocks: [a, b, c, d] */
148    /* Window, shuffle, fold */
149    {
150       /* Temp pointers to make it really clear to the compiler what we're doing */
151       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in+(overlap>>1);
152       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+N2-1+(overlap>>1);
153       kiss_fft_scalar * OPUS_RESTRICT yp = f;
154       const opus_val16 * OPUS_RESTRICT wp1 = window+(overlap>>1);
155       const opus_val16 * OPUS_RESTRICT wp2 = window+(overlap>>1)-1;
156       for(i=0;i<((overlap+3)>>2);i++)
157       {
158          /* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
159          *yp++ = MULT16_32_Q15(*wp2, xp1[N2]) + MULT16_32_Q15(*wp1,*xp2);
160          *yp++ = MULT16_32_Q15(*wp1, *xp1)    - MULT16_32_Q15(*wp2, xp2[-N2]);
161          xp1+=2;
162          xp2-=2;
163          wp1+=2;
164          wp2-=2;
165       }
166       wp1 = window;
167       wp2 = window+overlap-1;
168       for(;i<N4-((overlap+3)>>2);i++)
169       {
170          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
171          *yp++ = *xp2;
172          *yp++ = *xp1;
173          xp1+=2;
174          xp2-=2;
175       }
176       for(;i<N4;i++)
177       {
178          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
179          *yp++ =  -MULT16_32_Q15(*wp1, xp1[-N2]) + MULT16_32_Q15(*wp2, *xp2);
180          *yp++ = MULT16_32_Q15(*wp2, *xp1)     + MULT16_32_Q15(*wp1, xp2[N2]);
181          xp1+=2;
182          xp2-=2;
183          wp1+=2;
184          wp2-=2;
185       }
186    }
187    /* Pre-rotation */
188    {
189       kiss_fft_scalar * OPUS_RESTRICT yp = f;
190       const kiss_twiddle_scalar *t = &trig[0];
191       for(i=0;i<N4;i++)
192       {
193          kiss_fft_cpx yc;
194          kiss_twiddle_scalar t0, t1;
195          kiss_fft_scalar re, im, yr, yi;
196          t0 = t[i];
197          t1 = t[N4+i];
198 #ifdef FIXED_POINT
199          t0 = MULT16_16_P15(t0, scale);
200          t1 = MULT16_16_P15(t1, scale);
201 #else
202          t0 *= st->scale;
203          t1 *= st->scale;
204 #endif
205          re = *yp++;
206          im = *yp++;
207          yr = -S_MUL(re,t0)  +  S_MUL(im,t1);
208          yi = -S_MUL(im,t0)  -  S_MUL(re,t1);
209          yc.r = yr;
210          yc.i = yi;
211 #ifdef FIXED_POINT
212          yc.r = SHR32(yc.r, scale_shift);
213          yc.i = SHR32(yc.i, scale_shift);
214 #endif
215          f2[st->bitrev[i]] = yc;
216       }
217    }
218
219    /* N/4 complex FFT, down-scales by 4/N */
220    opus_fft_impl(st, f2);
221
222    /* Post-rotate */
223    {
224       /* Temp pointers to make it really clear to the compiler what we're doing */
225       const kiss_fft_cpx * OPUS_RESTRICT fp = f2;
226       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
227       kiss_fft_scalar * OPUS_RESTRICT yp2 = out+stride*(N2-1);
228       const kiss_twiddle_scalar *t = &trig[0];
229       /* Temp pointers to make it really clear to the compiler what we're doing */
230       for(i=0;i<N4;i++)
231       {
232          kiss_fft_scalar yr, yi;
233          yr = -S_MUL(fp->i,t[N4+i]) + S_MUL(fp->r,t[i]);
234          yi = -S_MUL(fp->r,t[N4+i]) - S_MUL(fp->i,t[i]);
235          *yp1 = yr;
236          *yp2 = yi;
237          fp++;
238          yp1 += 2*stride;
239          yp2 -= 2*stride;
240       }
241    }
242    RESTORE_STACK;
243 }
244
245 void clt_mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * OPUS_RESTRICT out,
246       const opus_val16 * OPUS_RESTRICT window, int overlap, int shift, int stride)
247 {
248    int i;
249    int N, N2, N4;
250    const kiss_twiddle_scalar *trig;
251
252    N = l->n;
253    trig = l->trig;
254    for (i=0;i<shift;i++)
255    {
256       N >>= 1;
257       trig += N;
258    }
259    N2 = N>>1;
260    N4 = N>>2;
261
262    /* Pre-rotate */
263    {
264       /* Temp pointers to make it really clear to the compiler what we're doing */
265       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in;
266       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+stride*(N2-1);
267       kiss_fft_scalar * OPUS_RESTRICT yp = out+(overlap>>1);
268       const kiss_twiddle_scalar * OPUS_RESTRICT t = &trig[0];
269       const opus_int16 * OPUS_RESTRICT bitrev = l->kfft[shift]->bitrev;
270       for(i=0;i<N4;i++)
271       {
272          int rev;
273          kiss_fft_scalar yr, yi;
274          rev = *bitrev++;
275          yr = -S_MUL(*xp2, t[i]) - S_MUL(*xp1,t[N4+i]);
276          yi =  S_MUL(*xp2, t[N4+i]) - S_MUL(*xp1,t[i]);
277          /* We swap real and imag because we use an FFT instead of an IFFT. */
278          yp[2*rev+1] = yr;
279          yp[2*rev] = yi;
280          /* Storing the pre-rotation directly in the bitrev order. */
281          xp1+=2*stride;
282          xp2-=2*stride;
283       }
284    }
285
286    opus_fft_impl(l->kfft[shift], (kiss_fft_cpx*)(out+(overlap>>1)));
287
288    /* Post-rotate and de-shuffle from both ends of the buffer at once to make
289       it in-place. */
290    {
291       kiss_fft_scalar * OPUS_RESTRICT yp0 = out+(overlap>>1);
292       kiss_fft_scalar * OPUS_RESTRICT yp1 = out+(overlap>>1)+N2-2;
293       const kiss_twiddle_scalar *t = &trig[0];
294       /* Loop to (N4+1)>>1 to handle odd N4. When N4 is odd, the
295          middle pair will be computed twice. */
296       for(i=0;i<(N4+1)>>1;i++)
297       {
298          kiss_fft_scalar re, im, yr, yi;
299          kiss_twiddle_scalar t0, t1;
300          /* We swap real and imag because we're using an FFT instead of an IFFT. */
301          re = yp0[1];
302          im = yp0[0];
303          t0 = t[i];
304          t1 = t[N4+i];
305          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
306          yr = S_MUL(re,t0) + S_MUL(im,t1);
307          yi = S_MUL(im,t0) - S_MUL(re,t1);
308          /* We swap real and imag because we're using an FFT instead of an IFFT. */
309          re = yp1[1];
310          im = yp1[0];
311          yp0[0] = -yr;
312          yp1[1] = yi;
313
314          t0 = t[(N4-i-1)];
315          t1 = t[(N2-i-1)];
316          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
317          yr = S_MUL(re,t0) + S_MUL(im,t1);
318          yi = S_MUL(im,t0) - S_MUL(re,t1);
319          yp1[0] = -yr;
320          yp0[1] = yi;
321          yp0 += 2;
322          yp1 -= 2;
323       }
324    }
325
326    /* Mirror on both sides for TDAC */
327    {
328       kiss_fft_scalar * OPUS_RESTRICT xp1 = out+overlap-1;
329       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
330       const opus_val16 * OPUS_RESTRICT wp1 = window;
331       const opus_val16 * OPUS_RESTRICT wp2 = window+overlap-1;
332
333       for(i = 0; i < overlap/2; i++)
334       {
335          kiss_fft_scalar x1, x2;
336          x1 = *xp1;
337          x2 = *yp1;
338          *yp1++ = MULT16_32_Q15(*wp2, x2) - MULT16_32_Q15(*wp1, x1);
339          *xp1-- = MULT16_32_Q15(*wp1, x2) + MULT16_32_Q15(*wp2, x1);
340          wp1++;
341          wp2--;
342       }
343    }
344 }