Mercurial > mplayer.hg
comparison libfaad2/mdct.c @ 10725:e989150f8216
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 |