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