Merge branch 'exp_mips_opt'
[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 #if defined(MIPSr1_ASM)
57 #include "mips/mdct_mipsr1.h"
58 #endif
59
60
61 #ifdef CUSTOM_MODES
62
63 int clt_mdct_init(mdct_lookup *l,int N, int maxshift)
64 {
65    int i;
66    kiss_twiddle_scalar *trig;
67    int shift;
68    int N2=N>>1;
69    l->n = N;
70    l->maxshift = maxshift;
71    for (i=0;i<=maxshift;i++)
72    {
73       if (i==0)
74          l->kfft[i] = opus_fft_alloc(N>>2>>i, 0, 0);
75       else
76          l->kfft[i] = opus_fft_alloc_twiddles(N>>2>>i, 0, 0, l->kfft[0]);
77 #ifndef ENABLE_TI_DSPLIB55
78       if (l->kfft[i]==NULL)
79          return 0;
80 #endif
81    }
82    l->trig = trig = (kiss_twiddle_scalar*)opus_alloc((N-(N2>>maxshift))*sizeof(kiss_twiddle_scalar));
83    if (l->trig==NULL)
84      return 0;
85    for (shift=0;shift<=maxshift;shift++)
86    {
87       /* We have enough points that sine isn't necessary */
88 #if defined(FIXED_POINT)
89 #if 1
90       for (i=0;i<N2;i++)
91          trig[i] = TRIG_UPSCALE*celt_cos_norm(DIV32(ADD32(SHL32(EXTEND32(i),17),N2+16384),N));
92 #else
93       for (i=0;i<N2;i++)
94          trig[i] = (kiss_twiddle_scalar)MAX32(-32767,MIN32(32767,floor(.5+32768*cos(2*M_PI*(i+.125)/N))));
95 #endif
96 #else
97       for (i=0;i<N2;i++)
98          trig[i] = (kiss_twiddle_scalar)cos(2*PI*(i+.125)/N);
99 #endif
100       trig += N2;
101       N2 >>= 1;
102       N >>= 1;
103    }
104    return 1;
105 }
106
107 void clt_mdct_clear(mdct_lookup *l)
108 {
109    int i;
110    for (i=0;i<=l->maxshift;i++)
111       opus_fft_free(l->kfft[i]);
112    opus_free((kiss_twiddle_scalar*)l->trig);
113 }
114
115 #endif /* CUSTOM_MODES */
116
117 /* Forward MDCT trashes the input array */
118 #ifndef OVERRIDE_clt_mdct_forward
119 void clt_mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * OPUS_RESTRICT out,
120       const opus_val16 *window, int overlap, int shift, int stride)
121 {
122    int i;
123    int N, N2, N4;
124    VARDECL(kiss_fft_scalar, f);
125    VARDECL(kiss_fft_cpx, f2);
126    const kiss_fft_state *st = l->kfft[shift];
127    const kiss_twiddle_scalar *trig;
128    opus_val16 scale;
129 #ifdef FIXED_POINT
130    /* Allows us to scale with MULT16_32_Q16(), which is faster than
131       MULT16_32_Q15() on ARM. */
132    int scale_shift = st->scale_shift-1;
133 #endif
134    SAVE_STACK;
135    scale = st->scale;
136
137    N = l->n;
138    trig = l->trig;
139    for (i=0;i<shift;i++)
140    {
141       N >>= 1;
142       trig += N;
143    }
144    N2 = N>>1;
145    N4 = N>>2;
146
147    ALLOC(f, N2, kiss_fft_scalar);
148    ALLOC(f2, N4, kiss_fft_cpx);
149
150    /* Consider the input to be composed of four blocks: [a, b, c, d] */
151    /* Window, shuffle, fold */
152    {
153       /* Temp pointers to make it really clear to the compiler what we're doing */
154       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in+(overlap>>1);
155       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+N2-1+(overlap>>1);
156       kiss_fft_scalar * OPUS_RESTRICT yp = f;
157       const opus_val16 * OPUS_RESTRICT wp1 = window+(overlap>>1);
158       const opus_val16 * OPUS_RESTRICT wp2 = window+(overlap>>1)-1;
159       for(i=0;i<((overlap+3)>>2);i++)
160       {
161          /* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
162          *yp++ = MULT16_32_Q15(*wp2, xp1[N2]) + MULT16_32_Q15(*wp1,*xp2);
163          *yp++ = MULT16_32_Q15(*wp1, *xp1)    - MULT16_32_Q15(*wp2, xp2[-N2]);
164          xp1+=2;
165          xp2-=2;
166          wp1+=2;
167          wp2-=2;
168       }
169       wp1 = window;
170       wp2 = window+overlap-1;
171       for(;i<N4-((overlap+3)>>2);i++)
172       {
173          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
174          *yp++ = *xp2;
175          *yp++ = *xp1;
176          xp1+=2;
177          xp2-=2;
178       }
179       for(;i<N4;i++)
180       {
181          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
182          *yp++ =  -MULT16_32_Q15(*wp1, xp1[-N2]) + MULT16_32_Q15(*wp2, *xp2);
183          *yp++ = MULT16_32_Q15(*wp2, *xp1)     + MULT16_32_Q15(*wp1, xp2[N2]);
184          xp1+=2;
185          xp2-=2;
186          wp1+=2;
187          wp2-=2;
188       }
189    }
190    /* Pre-rotation */
191    {
192       kiss_fft_scalar * OPUS_RESTRICT yp = f;
193       const kiss_twiddle_scalar *t = &trig[0];
194       for(i=0;i<N4;i++)
195       {
196          kiss_fft_cpx yc;
197          kiss_twiddle_scalar t0, t1;
198          kiss_fft_scalar re, im, yr, yi;
199          t0 = t[i];
200          t1 = t[N4+i];
201          re = *yp++;
202          im = *yp++;
203          yr = S_MUL(re,t0)  -  S_MUL(im,t1);
204          yi = S_MUL(im,t0)  +  S_MUL(re,t1);
205          yc.r = yr;
206          yc.i = yi;
207          yc.r = PSHR32(MULT16_32_Q16(scale, yc.r), scale_shift);
208          yc.i = PSHR32(MULT16_32_Q16(scale, yc.i), scale_shift);
209          f2[st->bitrev[i]] = yc;
210       }
211    }
212
213    /* N/4 complex FFT, does not downscale anymore */
214    opus_fft_impl(st, f2);
215
216    /* Post-rotate */
217    {
218       /* Temp pointers to make it really clear to the compiler what we're doing */
219       const kiss_fft_cpx * OPUS_RESTRICT fp = f2;
220       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
221       kiss_fft_scalar * OPUS_RESTRICT yp2 = out+stride*(N2-1);
222       const kiss_twiddle_scalar *t = &trig[0];
223       /* Temp pointers to make it really clear to the compiler what we're doing */
224       for(i=0;i<N4;i++)
225       {
226          kiss_fft_scalar yr, yi;
227          yr = S_MUL(fp->i,t[N4+i]) - S_MUL(fp->r,t[i]);
228          yi = S_MUL(fp->r,t[N4+i]) + S_MUL(fp->i,t[i]);
229          *yp1 = yr;
230          *yp2 = yi;
231          fp++;
232          yp1 += 2*stride;
233          yp2 -= 2*stride;
234       }
235    }
236    RESTORE_STACK;
237 }
238 #endif /* OVERRIDE_clt_mdct_forward */
239
240 #ifndef OVERRIDE_clt_mdct_backward
241 void clt_mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * OPUS_RESTRICT out,
242       const opus_val16 * OPUS_RESTRICT window, int overlap, int shift, int stride)
243 {
244    int i;
245    int N, N2, N4;
246    const kiss_twiddle_scalar *trig;
247
248    N = l->n;
249    trig = l->trig;
250    for (i=0;i<shift;i++)
251    {
252       N >>= 1;
253       trig += N;
254    }
255    N2 = N>>1;
256    N4 = N>>2;
257
258    /* Pre-rotate */
259    {
260       /* Temp pointers to make it really clear to the compiler what we're doing */
261       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in;
262       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+stride*(N2-1);
263       kiss_fft_scalar * OPUS_RESTRICT yp = out+(overlap>>1);
264       const kiss_twiddle_scalar * OPUS_RESTRICT t = &trig[0];
265       const opus_int16 * OPUS_RESTRICT bitrev = l->kfft[shift]->bitrev;
266       for(i=0;i<N4;i++)
267       {
268          int rev;
269          kiss_fft_scalar yr, yi;
270          rev = *bitrev++;
271          yr = S_MUL(*xp2, t[i]) + S_MUL(*xp1, t[N4+i]);
272          yi = S_MUL(*xp1, t[i]) - S_MUL(*xp2, t[N4+i]);
273          /* We swap real and imag because we use an FFT instead of an IFFT. */
274          yp[2*rev+1] = yr;
275          yp[2*rev] = yi;
276          /* Storing the pre-rotation directly in the bitrev order. */
277          xp1+=2*stride;
278          xp2-=2*stride;
279       }
280    }
281
282    opus_fft_impl(l->kfft[shift], (kiss_fft_cpx*)(out+(overlap>>1)));
283
284    /* Post-rotate and de-shuffle from both ends of the buffer at once to make
285       it in-place. */
286    {
287       kiss_fft_scalar * yp0 = out+(overlap>>1);
288       kiss_fft_scalar * yp1 = out+(overlap>>1)+N2-2;
289       const kiss_twiddle_scalar *t = &trig[0];
290       /* Loop to (N4+1)>>1 to handle odd N4. When N4 is odd, the
291          middle pair will be computed twice. */
292       for(i=0;i<(N4+1)>>1;i++)
293       {
294          kiss_fft_scalar re, im, yr, yi;
295          kiss_twiddle_scalar t0, t1;
296          /* We swap real and imag because we're using an FFT instead of an IFFT. */
297          re = yp0[1];
298          im = yp0[0];
299          t0 = t[i];
300          t1 = t[N4+i];
301          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
302          yr = S_MUL(re,t0) + S_MUL(im,t1);
303          yi = S_MUL(re,t1) - S_MUL(im,t0);
304          /* We swap real and imag because we're using an FFT instead of an IFFT. */
305          re = yp1[1];
306          im = yp1[0];
307          yp0[0] = yr;
308          yp1[1] = yi;
309
310          t0 = t[(N4-i-1)];
311          t1 = t[(N2-i-1)];
312          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
313          yr = S_MUL(re,t0) + S_MUL(im,t1);
314          yi = S_MUL(re,t1) - S_MUL(im,t0);
315          yp1[0] = yr;
316          yp0[1] = yi;
317          yp0 += 2;
318          yp1 -= 2;
319       }
320    }
321
322    /* Mirror on both sides for TDAC */
323    {
324       kiss_fft_scalar * OPUS_RESTRICT xp1 = out+overlap-1;
325       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
326       const opus_val16 * OPUS_RESTRICT wp1 = window;
327       const opus_val16 * OPUS_RESTRICT wp2 = window+overlap-1;
328
329       for(i = 0; i < overlap/2; i++)
330       {
331          kiss_fft_scalar x1, x2;
332          x1 = *xp1;
333          x2 = *yp1;
334          *yp1++ = MULT16_32_Q15(*wp2, x2) - MULT16_32_Q15(*wp1, x1);
335          *xp1-- = MULT16_32_Q15(*wp1, x2) + MULT16_32_Q15(*wp2, x1);
336          wp1++;
337          wp2--;
338       }
339    }
340 }
341 #endif /* OVERRIDE_clt_mdct_backward */