Mercurial > mplayer.hg

comparison libfaad2/mdct.c @ 10725:e989150f8216

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
libfaad2 v2.0rc1 imported
author arpi
date Sat, 30 Aug 2003 22:30:28 +0000
parents
children 3185f64f6350
comparison
equal deleted inserted replaced
10724:adf5697b9d83 10725:e989150f8216
1 /*
2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
3 ** Copyright (C) 2003 M. Bakker, Ahead Software AG, http://www.nero.com
4 **
5 ** This program is free software; you can redistribute it and/or modify
6 ** it under the terms of the GNU General Public License as published by
7 ** the Free Software Foundation; either version 2 of the License, or
8 ** (at your option) any later version.
9 **
10 ** This program is distributed in the hope that it will be useful,
11 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
12 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 ** GNU General Public License for more details.
14 **
15 ** You should have received a copy of the GNU General Public License
16 ** along with this program; if not, write to the Free Software
17 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18 **
19 ** Any non-GPL usage of this software or parts of this software is strictly
20 ** forbidden.
21 **
22 ** Commercial non-GPL licensing of this software is possible.
23 ** For more info contact Ahead Software through Mpeg4AAClicense@nero.com.
24 **
25 ** $Id: mdct.c,v 1.26 2003/07/29 08:20:12 menno Exp $
26 **/
27
28 /*
29 * Fast (I)MDCT Implementation using (I)FFT ((Inverse) Fast Fourier Transform)
30 * and consists of three steps: pre-(I)FFT complex multiplication, complex
31 * (I)FFT, post-(I)FFT complex multiplication,
32 *
33 * As described in:
34 * P. Duhamel, Y. Mahieux, and J.P. Petit, "A Fast Algorithm for the
35 * Implementation of Filter Banks Based on 'Time Domain Aliasing
36 * Cancellation’," IEEE Proc. on ICASSP‘91, 1991, pp. 2209-2212.
37 *
38 *
39 * As of April 6th 2002 completely rewritten.
40 * This (I)MDCT can now be used for any data size n, where n is divisible by 8.
41 *
42 */
43
44 #include "common.h"
45 #include "structs.h"
46
47 #include <stdlib.h>
48 #ifdef _WIN32_WCE
49 #define assert(x)
50 #else
51 #include <assert.h>
52 #endif
53
54 #include "cfft.h"
55 #include "mdct.h"
56
57 /* const_tab[]:
58 0: sqrt(2 / N)
59 1: cos(2 * PI / N)
60 2: sin(2 * PI / N)
61 3: cos(2 * PI * (1/8) / N)
62 4: sin(2 * PI * (1/8) / N)
63 */
64 #ifndef FIXED_POINT
65 #ifdef _MSC_VER
66 #pragma warning(disable:4305)
67 #pragma warning(disable:4244)
68 #endif
69 real_t const_tab[][5] =
70 {
71 { COEF_CONST(0.0312500000), COEF_CONST(0.9999952938), COEF_CONST(0.0030679568),
72 COEF_CONST(0.9999999265), COEF_CONST(0.0003834952) }, /* 2048 */
73 { COEF_CONST(0.0322748612), COEF_CONST(0.9999946356), COEF_CONST(0.0032724866),
74 COEF_CONST(0.9999999404), COEF_CONST(0.0004090615) }, /* 1920 */
75 { COEF_CONST(0.0441941738), COEF_CONST(0.9999811649), COEF_CONST(0.0061358847),
76 COEF_CONST(0.9999997020), COEF_CONST(0.0007669903) }, /* 1024 */
77 { COEF_CONST(0.0456435465), COEF_CONST(0.9999786019), COEF_CONST(0.0065449383),
78 COEF_CONST(0.9999996424), COEF_CONST(0.0008181230) }, /* 960 */
79 { COEF_CONST(0.0883883476), COEF_CONST(0.9996988177), COEF_CONST(0.0245412290),
80 COEF_CONST(0.9999952912), COEF_CONST(0.0030679568) }, /* 256 */
81 { COEF_CONST(0.0912870929), COEF_CONST(0.9996573329), COEF_CONST(0.0261769500),
82 COEF_CONST(0.9999946356), COEF_CONST(0.0032724866) } /* 240 */
83 #ifdef SSR_DEC
84 ,{ COEF_CONST(0.062500000), COEF_CONST(0.999924702), COEF_CONST(0.012271538),
85 COEF_CONST(0.999998823), COEF_CONST(0.00153398) }, /* 512 */
86 { COEF_CONST(0.176776695), COEF_CONST(0.995184727), COEF_CONST(0.09801714),
87 COEF_CONST(0.999924702), COEF_CONST(0.012271538) } /* 64 */
88 #endif
89 };
90 #else
91 real_t const_tab[][5] =
92 {
93 { COEF_CONST(1), COEF_CONST(0.9999952938), COEF_CONST(0.0030679568),
94 COEF_CONST(0.9999999265), COEF_CONST(0.0003834952) }, /* 2048 */
95 { COEF_CONST(/* sqrt(1024/960) */ 1.03279556), COEF_CONST(0.9999946356), COEF_CONST(0.0032724866),
96 COEF_CONST(0), COEF_CONST(0.0004090615) }, /* 1920 */
97 { COEF_CONST(1), COEF_CONST(0.9999811649), COEF_CONST(0.0061358847),
98 COEF_CONST(0.9999997020), COEF_CONST(0.0007669903) }, /* 1024 */
99 { COEF_CONST(/* sqrt(512/480) */ 1.03279556), COEF_CONST(0.9999786019), COEF_CONST(0.0065449383),
100 COEF_CONST(0.9999996424), COEF_CONST(0.0008181230) }, /* 960 */
101 { COEF_CONST(1), COEF_CONST(0.9996988177), COEF_CONST(0.0245412290),
102 COEF_CONST(0.9999952912), COEF_CONST(0.0030679568) }, /* 256 */
103 { COEF_CONST(/* sqrt(256/240) */ 1.03279556), COEF_CONST(0.9996573329), COEF_CONST(0.0261769500),
104 COEF_CONST(0.9999946356), COEF_CONST(0.0032724866) } /* 240 */
105 #ifdef SSR_DEC
106 ,{ COEF_CONST(0), COEF_CONST(0.999924702), COEF_CONST(0.012271538),
107 COEF_CONST(0.999998823), COEF_CONST(0.00153398) }, /* 512 */
108 { COEF_CONST(0), COEF_CONST(0.995184727), COEF_CONST(0.09801714),
109 COEF_CONST(0.999924702), COEF_CONST(0.012271538) } /* 64 */
110 #endif
111 };
112 #endif
113
114 uint8_t map_N_to_idx(uint16_t N)
115 {
116 switch(N)
117 {
118 case 2048: return 0;
119 case 1920: return 1;
120 case 1024: return 2;
121 case 960: return 3;
122 case 256: return 4;
123 case 240: return 5;
124 #ifdef SSR_DEC
125 case 512: return 6;
126 case 64: return 7;
127 #endif
128 }
129 return 0;
130 }
131
132 mdct_info *faad_mdct_init(uint16_t N)
133 {
134 uint16_t k, N_idx;
135 real_t cangle, sangle, c, s, cold;
136 real_t scale;
137
138 mdct_info *mdct = (mdct_info*)malloc(sizeof(mdct_info));
139
140 assert(N % 8 == 0);
141
142 mdct->N = N;
143 mdct->sincos = (complex_t*)malloc(N/4*sizeof(complex_t));
144 mdct->Z1 = (complex_t*)malloc(N/4*sizeof(complex_t));
145
146 N_idx = map_N_to_idx(N);
147
148 scale = const_tab[N_idx][0];
149 cangle = const_tab[N_idx][1];
150 sangle = const_tab[N_idx][2];
151 c = const_tab[N_idx][3];
152 s = const_tab[N_idx][4];
153
154 for (k = 0; k < N/4; k++)
155 {
156 RE(mdct->sincos[k]) = -1*MUL_C_C(c,scale);
157 IM(mdct->sincos[k]) = -1*MUL_C_C(s,scale);
158
159 cold = c;
160 c = MUL_C_C(c,cangle) - MUL_C_C(s,sangle);
161 s = MUL_C_C(s,cangle) + MUL_C_C(cold,sangle);
162 }
163
164 /* initialise fft */
165 mdct->cfft = cffti(N/4);
166
167 return mdct;
168 }
169
170 void faad_mdct_end(mdct_info *mdct)
171 {
172 if (mdct != NULL)
173 {
174 cfftu(mdct->cfft);
175
176 if (mdct->Z1) free(mdct->Z1);
177 if (mdct->sincos) free(mdct->sincos);
178
179 free(mdct);
180 }
181 }
182
183 void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
184 {
185 uint16_t k;
186
187 complex_t x;
188 complex_t *Z1 = mdct->Z1;
189 complex_t *sincos = mdct->sincos;
190
191 uint16_t N = mdct->N;
192 uint16_t N2 = N >> 1;
193 uint16_t N4 = N >> 2;
194 uint16_t N8 = N >> 3;
195
196 /* pre-IFFT complex multiplication */
197 for (k = 0; k < N4; k++)
198 {
199 uint16_t n = k << 1;
200 RE(x) = X_in[ n];
201 IM(x) = X_in[N2 - 1 - n];
202 RE(Z1[k]) = MUL_R_C(IM(x), RE(sincos[k])) - MUL_R_C(RE(x), IM(sincos[k]));
203 IM(Z1[k]) = MUL_R_C(RE(x), RE(sincos[k])) + MUL_R_C(IM(x), IM(sincos[k]));
204 }
205
206 /* complex IFFT */
207 cfftb(mdct->cfft, Z1);
208
209 /* post-IFFT complex multiplication */
210 for (k = 0; k < N4; k++)
211 {
212 uint16_t n = k << 1;
213 RE(x) = RE(Z1[k]);
214 IM(x) = IM(Z1[k]);
215
216 RE(Z1[k]) = MUL_R_C(RE(x), RE(sincos[k])) - MUL_R_C(IM(x), IM(sincos[k]));
217 IM(Z1[k]) = MUL_R_C(IM(x), RE(sincos[k])) + MUL_R_C(RE(x), IM(sincos[k]));
218 }
219
220 /* reordering */
221 for (k = 0; k < N8; k++)
222 {
223 uint16_t n = k << 1;
224 X_out[ n] = IM(Z1[N8 + k]);
225 X_out[ 1 + n] = -RE(Z1[N8 - 1 - k]);
226 X_out[N4 + n] = RE(Z1[ k]);
227 X_out[N4 + 1 + n] = -IM(Z1[N4 - 1 - k]);
228 X_out[N2 + n] = RE(Z1[N8 + k]);
229 X_out[N2 + 1 + n] = -IM(Z1[N8 - 1 - k]);
230 X_out[N2 + N4 + n] = -IM(Z1[ k]);
231 X_out[N2 + N4 + 1 + n] = RE(Z1[N4 - 1 - k]);
232 }
233 }
234
235 #ifdef LTP_DEC
236 void faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
237 {
238 uint16_t k;
239
240 complex_t x;
241 complex_t *Z1 = mdct->Z1;
242 complex_t *sincos = mdct->sincos;
243
244 uint16_t N = mdct->N;
245 uint16_t N2 = N >> 1;
246 uint16_t N4 = N >> 2;
247 uint16_t N8 = N >> 3;
248
249 real_t scale = REAL_CONST(N);
250
251 /* pre-FFT complex multiplication */
252 for (k = 0; k < N8; k++)
253 {
254 uint16_t n = k << 1;
255 RE(x) = X_in[N - N4 - 1 - n] + X_in[N - N4 + n];
256 IM(x) = X_in[ N4 + n] - X_in[ N4 - 1 - n];
257
258 RE(Z1[k]) = -MUL_R_C(RE(x), RE(sincos[k])) - MUL_R_C(IM(x), IM(sincos[k]));
259 IM(Z1[k]) = -MUL_R_C(IM(x), RE(sincos[k])) + MUL_R_C(RE(x), IM(sincos[k]));
260
261 RE(x) = X_in[N2 - 1 - n] - X_in[ n];
262 IM(x) = X_in[N2 + n] + X_in[N - 1 - n];
263
264 RE(Z1[k + N8]) = -MUL_R_C(RE(x), RE(sincos[k + N8])) - MUL_R_C(IM(x), IM(sincos[k + N8]));
265 IM(Z1[k + N8]) = -MUL_R_C(IM(x), RE(sincos[k + N8])) + MUL_R_C(RE(x), IM(sincos[k + N8]));
266 }
267
268 /* complex FFT */
269 cfftf(mdct->cfft, Z1);
270
271 /* post-FFT complex multiplication */
272 for (k = 0; k < N4; k++)
273 {
274 uint16_t n = k << 1;
275 RE(x) = MUL(MUL_R_C(RE(Z1[k]), RE(sincos[k])) + MUL_R_C(IM(Z1[k]), IM(sincos[k])), scale);
276 IM(x) = MUL(MUL_R_C(IM(Z1[k]), RE(sincos[k])) - MUL_R_C(RE(Z1[k]), IM(sincos[k])), scale);
277
278 X_out[ n] = RE(x);
279 X_out[N2 - 1 - n] = -IM(x);
280 X_out[N2 + n] = IM(x);
281 X_out[N - 1 - n] = -RE(x);
282 }
283 }
284 #endif