MIPS optimizations
[opus.git] / celt / mips / mdct_mipsr1.h
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 #ifndef __MDCT_MIPSR1_H__
42 #define __MDCT_MIPSR1_H__
43
44 #ifndef SKIP_CONFIG_H
45 #ifdef HAVE_CONFIG_H
46 #include "config.h"
47 #endif
48 #endif
49
50 #include "mdct.h"
51 #include "kiss_fft.h"
52 #include "_kiss_fft_guts.h"
53 #include <math.h>
54 #include "os_support.h"
55 #include "mathops.h"
56 #include "stack_alloc.h"
57
58 /* Forward MDCT trashes the input array */
59 #define OVERRIDE_clt_mdct_forward
60 void clt_mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * OPUS_RESTRICT out,
61       const opus_val16 *window, int overlap, int shift, int stride)
62 {
63    int i;
64    int N, N2, N4;
65    VARDECL(kiss_fft_scalar, f);
66    VARDECL(kiss_fft_cpx, f2);
67    const kiss_fft_state *st = l->kfft[shift];
68    const kiss_twiddle_scalar *trig;
69    opus_val16 scale;
70 #ifdef FIXED_POINT
71    /* Allows us to scale with MULT16_32_Q16(), which is faster than
72       MULT16_32_Q15() on ARM. */
73    int scale_shift = st->scale_shift-1;
74 #endif
75    SAVE_STACK;
76    scale = st->scale;
77
78    N = l->n;
79    trig = l->trig;
80    for (i=0;i<shift;i++)
81    {
82       N >>= 1;
83       trig += N;
84    }
85    N2 = N>>1;
86    N4 = N>>2;
87
88    ALLOC(f, N2, kiss_fft_scalar);
89    ALLOC(f2, N4, kiss_fft_cpx);
90
91    /* Consider the input to be composed of four blocks: [a, b, c, d] */
92    /* Window, shuffle, fold */
93    {
94       /* Temp pointers to make it really clear to the compiler what we're doing */
95       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in+(overlap>>1);
96       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+N2-1+(overlap>>1);
97       kiss_fft_scalar * OPUS_RESTRICT yp = f;
98       const opus_val16 * OPUS_RESTRICT wp1 = window+(overlap>>1);
99       const opus_val16 * OPUS_RESTRICT wp2 = window+(overlap>>1)-1;
100       for(i=0;i<((overlap+3)>>2);i++)
101       {
102          /* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
103           *yp++ = S_MUL_ADD(*wp2, xp1[N2],*wp1,*xp2);
104           *yp++ = S_MUL_SUB(*wp1, *xp1,*wp2, xp2[-N2]);
105          xp1+=2;
106          xp2-=2;
107          wp1+=2;
108          wp2-=2;
109       }
110       wp1 = window;
111       wp2 = window+overlap-1;
112       for(;i<N4-((overlap+3)>>2);i++)
113       {
114          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
115          *yp++ = *xp2;
116          *yp++ = *xp1;
117          xp1+=2;
118          xp2-=2;
119       }
120       for(;i<N4;i++)
121       {
122          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
123           *yp++ =  S_MUL_SUB(*wp2, *xp2, *wp1, xp1[-N2]);
124           *yp++ = S_MUL_ADD(*wp2, *xp1, *wp1, xp2[N2]);
125          xp1+=2;
126          xp2-=2;
127          wp1+=2;
128          wp2-=2;
129       }
130    }
131    /* Pre-rotation */
132    {
133       kiss_fft_scalar * OPUS_RESTRICT yp = f;
134       const kiss_twiddle_scalar *t = &trig[0];
135       for(i=0;i<N4;i++)
136       {
137          kiss_fft_cpx yc;
138          kiss_twiddle_scalar t0, t1;
139          kiss_fft_scalar re, im, yr, yi;
140          t0 = t[i];
141          t1 = t[N4+i];
142          re = *yp++;
143          im = *yp++;
144
145          yr = S_MUL_SUB(re,t0,im,t1);
146          yi = S_MUL_ADD(im,t0,re,t1);
147
148          yc.r = yr;
149          yc.i = yi;
150          yc.r = PSHR32(MULT16_32_Q16(scale, yc.r), scale_shift);
151          yc.i = PSHR32(MULT16_32_Q16(scale, yc.i), scale_shift);
152          f2[st->bitrev[i]] = yc;
153       }
154    }
155
156    /* N/4 complex FFT, does not downscale anymore */
157    opus_fft_impl(st, f2);
158
159    /* Post-rotate */
160    {
161       /* Temp pointers to make it really clear to the compiler what we're doing */
162       const kiss_fft_cpx * OPUS_RESTRICT fp = f2;
163       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
164       kiss_fft_scalar * OPUS_RESTRICT yp2 = out+stride*(N2-1);
165       const kiss_twiddle_scalar *t = &trig[0];
166       /* Temp pointers to make it really clear to the compiler what we're doing */
167       for(i=0;i<N4;i++)
168       {
169          kiss_fft_scalar yr, yi;
170          yr = S_MUL_SUB(fp->i,t[N4+i] , fp->r,t[i]);
171          yi = S_MUL_ADD(fp->r,t[N4+i] ,fp->i,t[i]);
172          *yp1 = yr;
173          *yp2 = yi;
174          fp++;
175          yp1 += 2*stride;
176          yp2 -= 2*stride;
177       }
178    }
179    RESTORE_STACK;
180 }
181
182 #define OVERRIDE_clt_mdct_backward
183 void clt_mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * OPUS_RESTRICT out,
184       const opus_val16 * OPUS_RESTRICT window, int overlap, int shift, int stride)
185 {
186    int i;
187    int N, N2, N4;
188    const kiss_twiddle_scalar *trig;
189
190    N = l->n;
191    trig = l->trig;
192    for (i=0;i<shift;i++)
193    {
194       N >>= 1;
195       trig += N;
196    }
197    N2 = N>>1;
198    N4 = N>>2;
199
200    /* Pre-rotate */
201    {
202       /* Temp pointers to make it really clear to the compiler what we're doing */
203       const kiss_fft_scalar * OPUS_RESTRICT xp1 = in;
204       const kiss_fft_scalar * OPUS_RESTRICT xp2 = in+stride*(N2-1);
205       kiss_fft_scalar * OPUS_RESTRICT yp = out+(overlap>>1);
206       const kiss_twiddle_scalar * OPUS_RESTRICT t = &trig[0];
207       const opus_int16 * OPUS_RESTRICT bitrev = l->kfft[shift]->bitrev;
208       for(i=0;i<N4;i++)
209       {
210          int rev;
211          kiss_fft_scalar yr, yi;
212          rev = *bitrev++;
213          yr = S_MUL_ADD(*xp2, t[i] , *xp1, t[N4+i]);
214          yi = S_MUL_SUB(*xp1, t[i] , *xp2, t[N4+i]);
215          /* We swap real and imag because we use an FFT instead of an IFFT. */
216          yp[2*rev+1] = yr;
217          yp[2*rev] = yi;
218          /* Storing the pre-rotation directly in the bitrev order. */
219          xp1+=2*stride;
220          xp2-=2*stride;
221       }
222    }
223
224    opus_fft_impl(l->kfft[shift], (kiss_fft_cpx*)(out+(overlap>>1)));
225
226    /* Post-rotate and de-shuffle from both ends of the buffer at once to make
227       it in-place. */
228    {
229       kiss_fft_scalar * OPUS_RESTRICT yp0 = out+(overlap>>1);
230       kiss_fft_scalar * OPUS_RESTRICT yp1 = out+(overlap>>1)+N2-2;
231       const kiss_twiddle_scalar *t = &trig[0];
232       /* Loop to (N4+1)>>1 to handle odd N4. When N4 is odd, the
233          middle pair will be computed twice. */
234       for(i=0;i<(N4+1)>>1;i++)
235       {
236          kiss_fft_scalar re, im, yr, yi;
237          kiss_twiddle_scalar t0, t1;
238          /* We swap real and imag because we're using an FFT instead of an IFFT. */
239          re = yp0[1];
240          im = yp0[0];
241          t0 = t[i];
242          t1 = t[N4+i];
243          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
244          yr = S_MUL_ADD(re,t0 , im,t1);
245          yi = S_MUL_SUB(re,t1 , im,t0);
246          /* We swap real and imag because we're using an FFT instead of an IFFT. */
247          re = yp1[1];
248          im = yp1[0];
249          yp0[0] = yr;
250          yp1[1] = yi;
251
252          t0 = t[(N4-i-1)];
253          t1 = t[(N2-i-1)];
254          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
255          yr = S_MUL_ADD(re,t0,im,t1);
256          yi = S_MUL_SUB(re,t1,im,t0);
257          yp1[0] = yr;
258          yp0[1] = yi;
259          yp0 += 2;
260          yp1 -= 2;
261       }
262    }
263
264    /* Mirror on both sides for TDAC */
265    {
266       kiss_fft_scalar * OPUS_RESTRICT xp1 = out+overlap-1;
267       kiss_fft_scalar * OPUS_RESTRICT yp1 = out;
268       const opus_val16 * OPUS_RESTRICT wp1 = window;
269       const opus_val16 * OPUS_RESTRICT wp2 = window+overlap-1;
270
271       for(i = 0; i < overlap/2; i++)
272       {
273          kiss_fft_scalar x1, x2;
274          x1 = *xp1;
275          x2 = *yp1;
276          *yp1++ = MULT16_32_Q15(*wp2, x2) - MULT16_32_Q15(*wp1, x1);
277          *xp1-- = MULT16_32_Q15(*wp1, x2) + MULT16_32_Q15(*wp2, x1);
278          wp1++;
279          wp2--;
280       }
281    }
282 }
283 #endif /* __MDCT_MIPSR1_H__ */