Getting rid of some negations
[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    /* FIXME: This should eventually just go in the state. */
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    /* Allows us to scale with MULT16_32_Q16(), which is faster than
132       MULT16_32_Q15() on ARM. */
133    scale_shift--;
134 #else
135    scale = st->scale;
136 #endif
137    SAVE_STACK;
138
139    N = l->n;
140    trig = l->trig;
141    for (i=0;i<shift;i++)
142    {
143       N >>= 1;
144       trig += N;
145    }
146    N2 = N>>1;
147    N4 = N>>2;
148
149    ALLOC(f, N2, kiss_fft_scalar);
150    ALLOC(f2, N4, kiss_fft_cpx);
151
152    /* Consider the input to be composed of four blocks: [a, b, c, d] */
153    /* Window, shuffle, fold */
154    {
155       /* Temp pointers to make it really clear to the compiler what we're doing */
156       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in+(overlap>>1);
157       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+N2-1+(overlap>>1);
158       kiss_fft_scalar * OPUS_RESTRICT yp = f;
159       const opus_val16 * OPUS_RESTRICT wp1 = window+(overlap>>1);
160       const opus_val16 * OPUS_RESTRICT wp2 = window+(overlap>>1)-1;
161       for(i=0;i<((overlap+3)>>2);i++)
162       {
163          /* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
164          *yp++ = MULT16_32_Q15(*wp2, xp1[N2]) + MULT16_32_Q15(*wp1,*xp2);
165          *yp++ = MULT16_32_Q15(*wp1, *xp1)    - MULT16_32_Q15(*wp2, xp2[-N2]);
166          xp1+=2;
167          xp2-=2;
168          wp1+=2;
169          wp2-=2;
170       }
171       wp1 = window;
172       wp2 = window+overlap-1;
173       for(;i<N4-((overlap+3)>>2);i++)
174       {
175          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
176          *yp++ = *xp2;
177          *yp++ = *xp1;
178          xp1+=2;
179          xp2-=2;
180       }
181       for(;i<N4;i++)
182       {
183          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
184          *yp++ =  -MULT16_32_Q15(*wp1, xp1[-N2]) + MULT16_32_Q15(*wp2, *xp2);
185          *yp++ = MULT16_32_Q15(*wp2, *xp1)     + MULT16_32_Q15(*wp1, xp2[N2]);
186          xp1+=2;
187          xp2-=2;
188          wp1+=2;
189          wp2-=2;
190       }
191    }
192    /* Pre-rotation */
193    {
194       kiss_fft_scalar * OPUS_RESTRICT yp = f;
195       const kiss_twiddle_scalar *t = &trig[0];
196       for(i=0;i<N4;i++)
197       {
198          kiss_fft_cpx yc;
199          kiss_twiddle_scalar t0, t1;
200          kiss_fft_scalar re, im, yr, yi;
201          t0 = t[i];
202          t1 = t[N4+i];
203          re = *yp++;
204          im = *yp++;
205          yr = S_MUL(re,t0)  -  S_MUL(im,t1);
206          yi = S_MUL(im,t0)  +  S_MUL(re,t1);
207          yc.r = yr;
208          yc.i = yi;
209          yc.r = PSHR32(MULT16_32_Q16(scale, yc.r), scale_shift);
210          yc.i = PSHR32(MULT16_32_Q16(scale, yc.i), scale_shift);
211          f2[st->bitrev[i]] = yc;
212       }
213    }
214
215    /* N/4 complex FFT, does not downscale anymore */
216    opus_fft_impl(st, f2);
217
218    /* Post-rotate */
219    {
220       /* Temp pointers to make it really clear to the compiler what we're doing */
221       const kiss_fft_cpx * OPUS_RESTRICT fp = f2;
222       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
223       kiss_fft_scalar * OPUS_RESTRICT yp2 = out+stride*(N2-1);
224       const kiss_twiddle_scalar *t = &trig[0];
225       /* Temp pointers to make it really clear to the compiler what we're doing */
226       for(i=0;i<N4;i++)
227       {
228          kiss_fft_scalar yr, yi;
229          yr = S_MUL(fp->i,t[N4+i]) - S_MUL(fp->r,t[i]);
230          yi = S_MUL(fp->r,t[N4+i]) + S_MUL(fp->i,t[i]);
231          *yp1 = yr;
232          *yp2 = yi;
233          fp++;
234          yp1 += 2*stride;
235          yp2 -= 2*stride;
236       }
237    }
238    RESTORE_STACK;
239 }
240
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 * OPUS_RESTRICT yp0 = out+(overlap>>1);
288       kiss_fft_scalar * OPUS_RESTRICT 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 }