b3516a91e09d812af39a5f1bce7e3b7a788c2582
[ffmpeg.git] / libavcodec / libac3 / imdct.c
1 /* 
2  *  imdct.c
3  *
4  *      Copyright (C) Aaron Holtzman - May 1999
5  *
6  *  This file is part of ac3dec, a free Dolby AC-3 stream decoder.
7  *      
8  *  ac3dec is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2, or (at your option)
11  *  any later version.
12  *   
13  *  ac3dec is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *   
18  *  You should have received a copy of the GNU General Public License
19  *  along with GNU Make; see the file COPYING.  If not, write to
20  *  the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
21  *
22  *
23  */
24
25 #include "../common.h"
26 #include <math.h>
27 #include "ac3.h"
28 #include "ac3_internal.h"
29
30 void (* imdct_256) (float data[], float delay[]);
31 void (* imdct_512) (float data[], float delay[]);
32
33 typedef struct complex_s
34 {
35     float real;
36     float imag;
37 } complex_t;
38
39
40 /* 128 point bit-reverse LUT */
41 static uint8_t bit_reverse_512[] = {
42         0x00, 0x40, 0x20, 0x60, 0x10, 0x50, 0x30, 0x70, 
43         0x08, 0x48, 0x28, 0x68, 0x18, 0x58, 0x38, 0x78, 
44         0x04, 0x44, 0x24, 0x64, 0x14, 0x54, 0x34, 0x74, 
45         0x0c, 0x4c, 0x2c, 0x6c, 0x1c, 0x5c, 0x3c, 0x7c, 
46         0x02, 0x42, 0x22, 0x62, 0x12, 0x52, 0x32, 0x72, 
47         0x0a, 0x4a, 0x2a, 0x6a, 0x1a, 0x5a, 0x3a, 0x7a, 
48         0x06, 0x46, 0x26, 0x66, 0x16, 0x56, 0x36, 0x76, 
49         0x0e, 0x4e, 0x2e, 0x6e, 0x1e, 0x5e, 0x3e, 0x7e, 
50         0x01, 0x41, 0x21, 0x61, 0x11, 0x51, 0x31, 0x71, 
51         0x09, 0x49, 0x29, 0x69, 0x19, 0x59, 0x39, 0x79, 
52         0x05, 0x45, 0x25, 0x65, 0x15, 0x55, 0x35, 0x75, 
53         0x0d, 0x4d, 0x2d, 0x6d, 0x1d, 0x5d, 0x3d, 0x7d, 
54         0x03, 0x43, 0x23, 0x63, 0x13, 0x53, 0x33, 0x73, 
55         0x0b, 0x4b, 0x2b, 0x6b, 0x1b, 0x5b, 0x3b, 0x7b, 
56         0x07, 0x47, 0x27, 0x67, 0x17, 0x57, 0x37, 0x77, 
57         0x0f, 0x4f, 0x2f, 0x6f, 0x1f, 0x5f, 0x3f, 0x7f};
58
59 static uint8_t bit_reverse_256[] = {
60         0x00, 0x20, 0x10, 0x30, 0x08, 0x28, 0x18, 0x38, 
61         0x04, 0x24, 0x14, 0x34, 0x0c, 0x2c, 0x1c, 0x3c, 
62         0x02, 0x22, 0x12, 0x32, 0x0a, 0x2a, 0x1a, 0x3a, 
63         0x06, 0x26, 0x16, 0x36, 0x0e, 0x2e, 0x1e, 0x3e, 
64         0x01, 0x21, 0x11, 0x31, 0x09, 0x29, 0x19, 0x39, 
65         0x05, 0x25, 0x15, 0x35, 0x0d, 0x2d, 0x1d, 0x3d, 
66         0x03, 0x23, 0x13, 0x33, 0x0b, 0x2b, 0x1b, 0x3b, 
67         0x07, 0x27, 0x17, 0x37, 0x0f, 0x2f, 0x1f, 0x3f};
68
69 static complex_t buf[128];
70
71 /* Twiddle factor LUT */
72 static complex_t w_1[1];
73 static complex_t w_2[2];
74 static complex_t w_4[4];
75 static complex_t w_8[8];
76 static complex_t w_16[16];
77 static complex_t w_32[32];
78 static complex_t w_64[64];
79 static complex_t * w[7] = {w_1, w_2, w_4, w_8, w_16, w_32, w_64};
80
81 /* Twiddle factors for IMDCT */
82 static float xcos1[128];
83 static float xsin1[128];
84 static float xcos2[64];
85 static float xsin2[64];
86
87 /* Windowing function for Modified DCT - Thank you acroread */
88 float imdct_window[] = {
89         0.00014, 0.00024, 0.00037, 0.00051, 0.00067, 0.00086, 0.00107, 0.00130,
90         0.00157, 0.00187, 0.00220, 0.00256, 0.00297, 0.00341, 0.00390, 0.00443,
91         0.00501, 0.00564, 0.00632, 0.00706, 0.00785, 0.00871, 0.00962, 0.01061,
92         0.01166, 0.01279, 0.01399, 0.01526, 0.01662, 0.01806, 0.01959, 0.02121,
93         0.02292, 0.02472, 0.02662, 0.02863, 0.03073, 0.03294, 0.03527, 0.03770,
94         0.04025, 0.04292, 0.04571, 0.04862, 0.05165, 0.05481, 0.05810, 0.06153,
95         0.06508, 0.06878, 0.07261, 0.07658, 0.08069, 0.08495, 0.08935, 0.09389,
96         0.09859, 0.10343, 0.10842, 0.11356, 0.11885, 0.12429, 0.12988, 0.13563,
97         0.14152, 0.14757, 0.15376, 0.16011, 0.16661, 0.17325, 0.18005, 0.18699,
98         0.19407, 0.20130, 0.20867, 0.21618, 0.22382, 0.23161, 0.23952, 0.24757,
99         0.25574, 0.26404, 0.27246, 0.28100, 0.28965, 0.29841, 0.30729, 0.31626,
100         0.32533, 0.33450, 0.34376, 0.35311, 0.36253, 0.37204, 0.38161, 0.39126,
101         0.40096, 0.41072, 0.42054, 0.43040, 0.44030, 0.45023, 0.46020, 0.47019,
102         0.48020, 0.49022, 0.50025, 0.51028, 0.52031, 0.53033, 0.54033, 0.55031,
103         0.56026, 0.57019, 0.58007, 0.58991, 0.59970, 0.60944, 0.61912, 0.62873,
104         0.63827, 0.64774, 0.65713, 0.66643, 0.67564, 0.68476, 0.69377, 0.70269,
105         0.71150, 0.72019, 0.72877, 0.73723, 0.74557, 0.75378, 0.76186, 0.76981,
106         0.77762, 0.78530, 0.79283, 0.80022, 0.80747, 0.81457, 0.82151, 0.82831,
107         0.83496, 0.84145, 0.84779, 0.85398, 0.86001, 0.86588, 0.87160, 0.87716,
108         0.88257, 0.88782, 0.89291, 0.89785, 0.90264, 0.90728, 0.91176, 0.91610,
109         0.92028, 0.92432, 0.92822, 0.93197, 0.93558, 0.93906, 0.94240, 0.94560,
110         0.94867, 0.95162, 0.95444, 0.95713, 0.95971, 0.96217, 0.96451, 0.96674,
111         0.96887, 0.97089, 0.97281, 0.97463, 0.97635, 0.97799, 0.97953, 0.98099,
112         0.98236, 0.98366, 0.98488, 0.98602, 0.98710, 0.98811, 0.98905, 0.98994,
113         0.99076, 0.99153, 0.99225, 0.99291, 0.99353, 0.99411, 0.99464, 0.99513,
114         0.99558, 0.99600, 0.99639, 0.99674, 0.99706, 0.99736, 0.99763, 0.99788,
115         0.99811, 0.99831, 0.99850, 0.99867, 0.99882, 0.99895, 0.99908, 0.99919,
116         0.99929, 0.99938, 0.99946, 0.99953, 0.99959, 0.99965, 0.99969, 0.99974,
117         0.99978, 0.99981, 0.99984, 0.99986, 0.99988, 0.99990, 0.99992, 0.99993,
118         0.99994, 0.99995, 0.99996, 0.99997, 0.99998, 0.99998, 0.99998, 0.99999,
119         0.99999, 0.99999, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
120         1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000 };
121
122
123 static inline void swap_cmplx(complex_t *a, complex_t *b)
124 {
125     complex_t tmp;
126
127     tmp = *a;
128     *a = *b;
129     *b = tmp;
130 }
131
132
133
134 static inline complex_t cmplx_mult(complex_t a, complex_t b)
135 {
136     complex_t ret;
137
138     ret.real = a.real * b.real - a.imag * b.imag;
139     ret.imag = a.real * b.imag + a.imag * b.real;
140
141     return ret;
142 }
143
144 void
145 imdct_do_512(float data[],float delay[])
146 {
147     int i,k;
148     int p,q;
149     int m;
150     int two_m;
151     int two_m_plus_one;
152
153     float tmp_a_i;
154     float tmp_a_r;
155     float tmp_b_i;
156     float tmp_b_r;
157
158     float *data_ptr;
159     float *delay_ptr;
160     float *window_ptr;
161         
162     //
163     // 512 IMDCT with source and dest data in 'data'
164     //
165         
166     // Pre IFFT complex multiply plus IFFT cmplx conjugate 
167     for( i=0; i < 128; i++) {
168         /* z[i] = (X[256-2*i-1] + j * X[2*i]) * (xcos1[i] + j * xsin1[i]) ; */ 
169         buf[i].real =         (data[256-2*i-1] * xcos1[i])  -  (data[2*i]       * xsin1[i]);
170         buf[i].imag = -1.0 * ((data[2*i]       * xcos1[i])  +  (data[256-2*i-1] * xsin1[i]));
171     }
172
173     //Bit reversed shuffling
174     for(i=0; i<128; i++) {
175         k = bit_reverse_512[i];
176         if (k < i)
177             swap_cmplx(&buf[i],&buf[k]);
178     }
179
180     /* FFT Merge */
181     for (m=0; m < 7; m++) {
182         if(m)
183             two_m = (1 << m);
184         else
185             two_m = 1;
186
187         two_m_plus_one = (1 << (m+1));
188
189         for(k = 0; k < two_m; k++) {
190             for(i = 0; i < 128; i += two_m_plus_one) {
191                 p = k + i;
192                 q = p + two_m;
193                 tmp_a_r = buf[p].real;
194                 tmp_a_i = buf[p].imag;
195                 tmp_b_r = buf[q].real * w[m][k].real - buf[q].imag * w[m][k].imag;
196                 tmp_b_i = buf[q].imag * w[m][k].real + buf[q].real * w[m][k].imag;
197                 buf[p].real = tmp_a_r + tmp_b_r;
198                 buf[p].imag =  tmp_a_i + tmp_b_i;
199                 buf[q].real = tmp_a_r - tmp_b_r;
200                 buf[q].imag =  tmp_a_i - tmp_b_i;
201             }
202         }
203     }
204
205     /* Post IFFT complex multiply  plus IFFT complex conjugate*/
206     for( i=0; i < 128; i++) {
207         /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ; */
208         tmp_a_r =        buf[i].real;
209         tmp_a_i = -1.0 * buf[i].imag;
210         buf[i].real =(tmp_a_r * xcos1[i])  -  (tmp_a_i  * xsin1[i]);
211         buf[i].imag =(tmp_a_r * xsin1[i])  +  (tmp_a_i  * xcos1[i]);
212     }
213         
214     data_ptr = data;
215     delay_ptr = delay;
216     window_ptr = imdct_window;
217
218     /* Window and convert to real valued signal */
219     for(i=0; i< 64; i++) { 
220         *data_ptr++   = 2.0f * (-buf[64+i].imag   * *window_ptr++ + *delay_ptr++); 
221         *data_ptr++   = 2.0f * ( buf[64-i-1].real * *window_ptr++ + *delay_ptr++); 
222     }
223
224     for(i=0; i< 64; i++) { 
225         *data_ptr++  = 2.0f * (-buf[i].real       * *window_ptr++ + *delay_ptr++); 
226         *data_ptr++  = 2.0f * ( buf[128-i-1].imag * *window_ptr++ + *delay_ptr++); 
227     }
228
229     /* The trailing edge of the window goes into the delay line */
230     delay_ptr = delay;
231
232     for(i=0; i< 64; i++) { 
233         *delay_ptr++  = -buf[64+i].real   * *--window_ptr; 
234         *delay_ptr++  =  buf[64-i-1].imag * *--window_ptr; 
235     }
236
237     for(i=0; i<64; i++) {
238         *delay_ptr++  =  buf[i].imag       * *--window_ptr; 
239         *delay_ptr++  = -buf[128-i-1].real * *--window_ptr; 
240     }
241 }
242
243 void
244 imdct_do_256(float data[],float delay[])
245 {
246     int i,k;
247     int p,q;
248     int m;
249     int two_m;
250     int two_m_plus_one;
251
252     float tmp_a_i;
253     float tmp_a_r;
254     float tmp_b_i;
255     float tmp_b_r;
256
257     float *data_ptr;
258     float *delay_ptr;
259     float *window_ptr;
260
261     complex_t *buf_1, *buf_2;
262
263     buf_1 = &buf[0];
264     buf_2 = &buf[64];
265
266     /* Pre IFFT complex multiply plus IFFT cmplx conjugate */
267     for(k=0; k<64; k++) { 
268         /* X1[k] = X[2*k]  */
269         /* X2[k] = X[2*k+1]     */
270
271         p = 2 * (128-2*k-1);
272         q = 2 * (2 * k);
273
274         /* Z1[k] = (X1[128-2*k-1] + j * X1[2*k]) * (xcos2[k] + j * xsin2[k]); */ 
275         buf_1[k].real =         data[p] * xcos2[k] - data[q] * xsin2[k];
276         buf_1[k].imag = -1.0f * (data[q] * xcos2[k] + data[p] * xsin2[k]); 
277         /* Z2[k] = (X2[128-2*k-1] + j * X2[2*k]) * (xcos2[k] + j * xsin2[k]); */ 
278         buf_2[k].real =          data[p + 1] * xcos2[k] - data[q + 1] * xsin2[k];
279         buf_2[k].imag = -1.0f * ( data[q + 1] * xcos2[k] + data[p + 1] * xsin2[k]); 
280     }
281
282     //IFFT Bit reversed shuffling
283     for(i=0; i<64; i++) { 
284         k = bit_reverse_256[i];
285         if (k < i) {
286             swap_cmplx(&buf_1[i],&buf_1[k]);
287             swap_cmplx(&buf_2[i],&buf_2[k]);
288         }
289     }
290
291     /* FFT Merge */
292     for (m=0; m < 6; m++) {
293         two_m = (1 << m);
294         two_m_plus_one = (1 << (m+1));
295
296         //FIXME
297         if(m)
298             two_m = (1 << m);
299         else
300             two_m = 1;
301
302         for(k = 0; k < two_m; k++) {
303             for(i = 0; i < 64; i += two_m_plus_one) {
304                 p = k + i;
305                 q = p + two_m;
306                 //Do block 1
307                 tmp_a_r = buf_1[p].real;
308                 tmp_a_i = buf_1[p].imag;
309                 tmp_b_r = buf_1[q].real * w[m][k].real - buf_1[q].imag * w[m][k].imag;
310                 tmp_b_i = buf_1[q].imag * w[m][k].real + buf_1[q].real * w[m][k].imag;
311                 buf_1[p].real = tmp_a_r + tmp_b_r;
312                 buf_1[p].imag =  tmp_a_i + tmp_b_i;
313                 buf_1[q].real = tmp_a_r - tmp_b_r;
314                 buf_1[q].imag =  tmp_a_i - tmp_b_i;
315
316                 //Do block 2
317                 tmp_a_r = buf_2[p].real;
318                 tmp_a_i = buf_2[p].imag;
319                 tmp_b_r = buf_2[q].real * w[m][k].real - buf_2[q].imag * w[m][k].imag;
320                 tmp_b_i = buf_2[q].imag * w[m][k].real + buf_2[q].real * w[m][k].imag;
321                 buf_2[p].real = tmp_a_r + tmp_b_r;
322                 buf_2[p].imag =  tmp_a_i + tmp_b_i;
323                 buf_2[q].real = tmp_a_r - tmp_b_r;
324                 buf_2[q].imag =  tmp_a_i - tmp_b_i;
325             }
326         }
327     }
328
329     /* Post IFFT complex multiply */
330     for( i=0; i < 64; i++) {
331         /* y1[n] = z1[n] * (xcos2[n] + j * xs in2[n]) ; */ 
332         tmp_a_r =  buf_1[i].real;
333         tmp_a_i = -buf_1[i].imag;
334         buf_1[i].real =(tmp_a_r * xcos2[i])  -  (tmp_a_i  * xsin2[i]);
335         buf_1[i].imag =(tmp_a_r * xsin2[i])  +  (tmp_a_i  * xcos2[i]);
336         /* y2[n] = z2[n] * (xcos2[n] + j * xsin2[n]) ; */ 
337         tmp_a_r =  buf_2[i].real;
338         tmp_a_i = -buf_2[i].imag;
339         buf_2[i].real =(tmp_a_r * xcos2[i])  -  (tmp_a_i  * xsin2[i]);
340         buf_2[i].imag =(tmp_a_r * xsin2[i])  +  (tmp_a_i  * xcos2[i]);
341     }
342         
343     data_ptr = data;
344     delay_ptr = delay;
345     window_ptr = imdct_window;
346
347     /* Window and convert to real valued signal */
348     for(i=0; i< 64; i++) { 
349         *data_ptr++  = 2.0f * (-buf_1[i].imag      * *window_ptr++ + *delay_ptr++);
350         *data_ptr++  = 2.0f * ( buf_1[64-i-1].real * *window_ptr++ + *delay_ptr++);
351     }
352
353     for(i=0; i< 64; i++) {
354         *data_ptr++  = 2.0f * (-buf_1[i].real      * *window_ptr++ + *delay_ptr++);
355         *data_ptr++  = 2.0f * ( buf_1[64-i-1].imag * *window_ptr++ + *delay_ptr++);
356     }
357         
358     delay_ptr = delay;
359
360     for(i=0; i< 64; i++) {
361         *delay_ptr++ = -buf_2[i].real      * *--window_ptr;
362         *delay_ptr++ =  buf_2[64-i-1].imag * *--window_ptr;
363     }
364
365     for(i=0; i< 64; i++) {
366         *delay_ptr++ =  buf_2[i].imag      * *--window_ptr;
367         *delay_ptr++ = -buf_2[64-i-1].real * *--window_ptr;
368     }
369 }
370
371 void imdct_init (void)
372 {
373 #ifdef LIBAC3_MLIB
374     void imdct_do_256_mlib(float data[],float delay[]);
375     void imdct_do_512_mlib(float data[],float delay[]);
376
377     imdct_512 = imdct_do_512_mlib;
378     imdct_256 = imdct_do_256_mlib;
379 #else
380     int i, j, k;
381
382     /* Twiddle factors to turn IFFT into IMDCT */
383     for (i = 0; i < 128; i++) {
384         xcos1[i] = -cos ((M_PI / 2048) * (8 * i + 1));
385         xsin1[i] = -sin ((M_PI / 2048) * (8 * i + 1));
386     }
387
388     /* More twiddle factors to turn IFFT into IMDCT */
389     for (i = 0; i < 64; i++) {
390         xcos2[i] = -cos ((M_PI / 1024) * (8 * i + 1));
391         xsin2[i] = -sin ((M_PI / 1024) * (8 * i + 1));
392     }
393
394     for (i = 0; i < 7; i++) {
395         j = 1 << i;
396         for (k = 0; k < j; k++) {
397             w[i][k].real = cos (-M_PI * k / j);
398             w[i][k].imag = sin (-M_PI * k / j);
399         }
400     }
401     imdct_512 = imdct_do_512;
402     imdct_256 = imdct_do_256;
403 #endif
404 }