← Back to branch summary

~ubuntu-branches/ubuntu/quantal/linphone/quantal

« back to all changes in this revision

Viewing changes to ffmpeg/libavcodec/libac3/imdct.c

  • Committer: Bazaar Package Importer
  • Author(s): Samuel Mimram
  • Date: 2004-06-30 13:58:16 UTC
  • Revision ID: james.westby@ubuntu.com-20040630135816-wwx75gdlodkqbabb
Tags: upstream-0.12.2
Import upstream version 0.12.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
ffmpeg/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
}