Adds 3rd clause to CELT license
[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    - Neither the name of Internet Society, IETF or IETF Trust, nor the
17    names of specific contributors, may be used to endorse or promote
18    products derived from this software without specific prior written
19    permission.
20
21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
25    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 /* This is a simple MDCT implementation that uses a N/4 complex FFT
35    to do most of the work. It should be relatively straightforward to
36    plug in pretty much and FFT here.
37
38    This replaces the Vorbis FFT (and uses the exact same API), which
39    was a bit too messy and that was ending up duplicating code
40    (might as well use the same FFT everywhere).
41
42    The algorithm is similar to (and inspired from) Fabrice Bellard's
43    MDCT implementation in FFMPEG, but has differences in signs, ordering
44    and scaling in many places.
45 */
46
47 #ifndef SKIP_CONFIG_H
48 #ifdef HAVE_CONFIG_H
49 #include "config.h"
50 #endif
51 #endif
52
53 #include "mdct.h"
54 #include "kiss_fft.h"
55 #include "_kiss_fft_guts.h"
56 #include <math.h>
57 #include "os_support.h"
58 #include "mathops.h"
59 #include "stack_alloc.h"
60
61 #ifdef CUSTOM_MODES
62
63 int clt_mdct_init(mdct_lookup *l,int N, int maxshift)
64 {
65    int i;
66    int N4;
67    kiss_twiddle_scalar *trig;
68 #if defined(FIXED_POINT)
69    int N2=N>>1;
70 #endif
71    l->n = N;
72    N4 = N>>2;
73    l->maxshift = maxshift;
74    for (i=0;i<=maxshift;i++)
75    {
76       if (i==0)
77          l->kfft[i] = opus_fft_alloc(N>>2>>i, 0, 0);
78       else
79          l->kfft[i] = opus_fft_alloc_twiddles(N>>2>>i, 0, 0, l->kfft[0]);
80 #ifndef ENABLE_TI_DSPLIB55
81       if (l->kfft[i]==NULL)
82          return 0;
83 #endif
84    }
85    l->trig = trig = (kiss_twiddle_scalar*)opus_alloc((N4+1)*sizeof(kiss_twiddle_scalar));
86    if (l->trig==NULL)
87      return 0;
88    /* We have enough points that sine isn't necessary */
89 #if defined(FIXED_POINT)
90    for (i=0;i<=N4;i++)
91       trig[i] = TRIG_UPSCALE*celt_cos_norm(DIV32(ADD32(SHL32(EXTEND32(i),17),N2),N));
92 #else
93    for (i=0;i<=N4;i++)
94       trig[i] = (kiss_twiddle_scalar)cos(2*PI*i/N);
95 #endif
96    return 1;
97 }
98
99 void clt_mdct_clear(mdct_lookup *l)
100 {
101    int i;
102    for (i=0;i<=l->maxshift;i++)
103       opus_fft_free(l->kfft[i]);
104    opus_free((kiss_twiddle_scalar*)l->trig);
105 }
106
107 #endif /* CUSTOM_MODES */
108
109 /* Forward MDCT trashes the input array */
110 void clt_mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * restrict out,
111       const opus_val16 *window, int overlap, int shift, int stride)
112 {
113    int i;
114    int N, N2, N4;
115    kiss_twiddle_scalar sine;
116    VARDECL(kiss_fft_scalar, f);
117    SAVE_STACK;
118    N = l->n;
119    N >>= shift;
120    N2 = N>>1;
121    N4 = N>>2;
122    ALLOC(f, N2, kiss_fft_scalar);
123    /* sin(x) ~= x here */
124 #ifdef FIXED_POINT
125    sine = TRIG_UPSCALE*(QCONST16(0.7853981f, 15)+N2)/N;
126 #else
127    sine = (kiss_twiddle_scalar)2*PI*(.125f)/N;
128 #endif
129
130    /* Consider the input to be composed of four blocks: [a, b, c, d] */
131    /* Window, shuffle, fold */
132    {
133       /* Temp pointers to make it really clear to the compiler what we're doing */
134       const kiss_fft_scalar * restrict xp1 = in+(overlap>>1);
135       const kiss_fft_scalar * restrict xp2 = in+N2-1+(overlap>>1);
136       kiss_fft_scalar * restrict yp = f;
137       const opus_val16 * restrict wp1 = window+(overlap>>1);
138       const opus_val16 * restrict wp2 = window+(overlap>>1)-1;
139       for(i=0;i<(overlap>>2);i++)
140       {
141          /* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
142          *yp++ = MULT16_32_Q15(*wp2, xp1[N2]) + MULT16_32_Q15(*wp1,*xp2);
143          *yp++ = MULT16_32_Q15(*wp1, *xp1)    - MULT16_32_Q15(*wp2, xp2[-N2]);
144          xp1+=2;
145          xp2-=2;
146          wp1+=2;
147          wp2-=2;
148       }
149       wp1 = window;
150       wp2 = window+overlap-1;
151       for(;i<N4-(overlap>>2);i++)
152       {
153          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
154          *yp++ = *xp2;
155          *yp++ = *xp1;
156          xp1+=2;
157          xp2-=2;
158       }
159       for(;i<N4;i++)
160       {
161          /* Real part arranged as a-bR, Imag part arranged as -c-dR */
162          *yp++ =  -MULT16_32_Q15(*wp1, xp1[-N2]) + MULT16_32_Q15(*wp2, *xp2);
163          *yp++ = MULT16_32_Q15(*wp2, *xp1)     + MULT16_32_Q15(*wp1, xp2[N2]);
164          xp1+=2;
165          xp2-=2;
166          wp1+=2;
167          wp2-=2;
168       }
169    }
170    /* Pre-rotation */
171    {
172       kiss_fft_scalar * restrict yp = f;
173       const kiss_twiddle_scalar *t = &l->trig[0];
174       for(i=0;i<N4;i++)
175       {
176          kiss_fft_scalar re, im, yr, yi;
177          re = yp[0];
178          im = yp[1];
179          yr = -S_MUL(re,t[i<<shift])  -  S_MUL(im,t[(N4-i)<<shift]);
180          yi = -S_MUL(im,t[i<<shift])  +  S_MUL(re,t[(N4-i)<<shift]);
181          /* works because the cos is nearly one */
182          *yp++ = yr + S_MUL(yi,sine);
183          *yp++ = yi - S_MUL(yr,sine);
184       }
185    }
186
187    /* N/4 complex FFT, down-scales by 4/N */
188    opus_fft(l->kfft[shift], (kiss_fft_cpx *)f, (kiss_fft_cpx *)in);
189
190    /* Post-rotate */
191    {
192       /* Temp pointers to make it really clear to the compiler what we're doing */
193       const kiss_fft_scalar * restrict fp = in;
194       kiss_fft_scalar * restrict yp1 = out;
195       kiss_fft_scalar * restrict yp2 = out+stride*(N2-1);
196       const kiss_twiddle_scalar *t = &l->trig[0];
197       /* Temp pointers to make it really clear to the compiler what we're doing */
198       for(i=0;i<N4;i++)
199       {
200          kiss_fft_scalar yr, yi;
201          yr = S_MUL(fp[1],t[(N4-i)<<shift]) + S_MUL(fp[0],t[i<<shift]);
202          yi = S_MUL(fp[0],t[(N4-i)<<shift]) - S_MUL(fp[1],t[i<<shift]);
203          /* works because the cos is nearly one */
204          *yp1 = yr - S_MUL(yi,sine);
205          *yp2 = yi + S_MUL(yr,sine);;
206          fp += 2;
207          yp1 += 2*stride;
208          yp2 -= 2*stride;
209       }
210    }
211    RESTORE_STACK;
212 }
213
214 void clt_mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar * restrict out,
215       const opus_val16 * restrict window, int overlap, int shift, int stride)
216 {
217    int i;
218    int N, N2, N4;
219    kiss_twiddle_scalar sine;
220    VARDECL(kiss_fft_scalar, f);
221    VARDECL(kiss_fft_scalar, f2);
222    SAVE_STACK;
223    N = l->n;
224    N >>= shift;
225    N2 = N>>1;
226    N4 = N>>2;
227    ALLOC(f, N2, kiss_fft_scalar);
228    ALLOC(f2, N2, kiss_fft_scalar);
229    /* sin(x) ~= x here */
230 #ifdef FIXED_POINT
231    sine = TRIG_UPSCALE*(QCONST16(0.7853981f, 15)+N2)/N;
232 #else
233    sine = (kiss_twiddle_scalar)2*PI*(.125f)/N;
234 #endif
235
236    /* Pre-rotate */
237    {
238       /* Temp pointers to make it really clear to the compiler what we're doing */
239       const kiss_fft_scalar * restrict xp1 = in;
240       const kiss_fft_scalar * restrict xp2 = in+stride*(N2-1);
241       kiss_fft_scalar * restrict yp = f2;
242       const kiss_twiddle_scalar *t = &l->trig[0];
243       for(i=0;i<N4;i++)
244       {
245          kiss_fft_scalar yr, yi;
246          yr = -S_MUL(*xp2, t[i<<shift]) + S_MUL(*xp1,t[(N4-i)<<shift]);
247          yi =  -S_MUL(*xp2, t[(N4-i)<<shift]) - S_MUL(*xp1,t[i<<shift]);
248          /* works because the cos is nearly one */
249          *yp++ = yr - S_MUL(yi,sine);
250          *yp++ = yi + S_MUL(yr,sine);
251          xp1+=2*stride;
252          xp2-=2*stride;
253       }
254    }
255
256    /* Inverse N/4 complex FFT. This one should *not* downscale even in fixed-point */
257    opus_ifft(l->kfft[shift], (kiss_fft_cpx *)f2, (kiss_fft_cpx *)f);
258
259    /* Post-rotate */
260    {
261       kiss_fft_scalar * restrict fp = f;
262       const kiss_twiddle_scalar *t = &l->trig[0];
263
264       for(i=0;i<N4;i++)
265       {
266          kiss_fft_scalar re, im, yr, yi;
267          re = fp[0];
268          im = fp[1];
269          /* We'd scale up by 2 here, but instead it's done when mixing the windows */
270          yr = S_MUL(re,t[i<<shift]) - S_MUL(im,t[(N4-i)<<shift]);
271          yi = S_MUL(im,t[i<<shift]) + S_MUL(re,t[(N4-i)<<shift]);
272          /* works because the cos is nearly one */
273          *fp++ = yr - S_MUL(yi,sine);
274          *fp++ = yi + S_MUL(yr,sine);
275       }
276    }
277    /* De-shuffle the components for the middle of the window only */
278    {
279       const kiss_fft_scalar * restrict fp1 = f;
280       const kiss_fft_scalar * restrict fp2 = f+N2-1;
281       kiss_fft_scalar * restrict yp = f2;
282       for(i = 0; i < N4; i++)
283       {
284          *yp++ =-*fp1;
285          *yp++ = *fp2;
286          fp1 += 2;
287          fp2 -= 2;
288       }
289    }
290    out -= (N2-overlap)>>1;
291    /* Mirror on both sides for TDAC */
292    {
293       kiss_fft_scalar * restrict fp1 = f2+N4-1;
294       kiss_fft_scalar * restrict xp1 = out+N2-1;
295       kiss_fft_scalar * restrict yp1 = out+N4-overlap/2;
296       const opus_val16 * restrict wp1 = window;
297       const opus_val16 * restrict wp2 = window+overlap-1;
298       for(i = 0; i< N4-overlap/2; i++)
299       {
300          *xp1 = *fp1;
301          xp1--;
302          fp1--;
303       }
304       for(; i < N4; i++)
305       {
306          kiss_fft_scalar x1;
307          x1 = *fp1--;
308          *yp1++ +=-MULT16_32_Q15(*wp1, x1);
309          *xp1-- += MULT16_32_Q15(*wp2, x1);
310          wp1++;
311          wp2--;
312       }
313    }
314    {
315       kiss_fft_scalar * restrict fp2 = f2+N4;
316       kiss_fft_scalar * restrict xp2 = out+N2;
317       kiss_fft_scalar * restrict yp2 = out+N-1-(N4-overlap/2);
318       const opus_val16 * restrict wp1 = window;
319       const opus_val16 * restrict wp2 = window+overlap-1;
320       for(i = 0; i< N4-overlap/2; i++)
321       {
322          *xp2 = *fp2;
323          xp2++;
324          fp2++;
325       }
326       for(; i < N4; i++)
327       {
328          kiss_fft_scalar x2;
329          x2 = *fp2++;
330          *yp2--  = MULT16_32_Q15(*wp1, x2);
331          *xp2++  = MULT16_32_Q15(*wp2, x2);
332          wp1++;
333          wp2--;
334       }
335    }
336    RESTORE_STACK;
337 }