Mercurial > mplayer.hg
annotate libfaad2/mdct.c @ 27428:bea545589195
fixes spotted by diego
author | compn |
---|---|
date | Wed, 13 Aug 2008 18:22:10 +0000 |
parents | 59b6fa5b4201 |
children | e83eef58b30a |
rev | line source |
---|---|
10725 | 1 /* |
2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding | |
12527 | 3 ** Copyright (C) 2003-2004 M. Bakker, Ahead Software AG, http://www.nero.com |
10725 | 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 ** | |
18141 | 25 ** $Id: mdct.c,v 1.43 2004/09/04 14:56:28 menno Exp $ |
10725 | 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" | |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
56 #include "mdct_tab.h" |
10725 | 57 |
58 | |
59 mdct_info *faad_mdct_init(uint16_t N) | |
60 { | |
12527 | 61 mdct_info *mdct = (mdct_info*)faad_malloc(sizeof(mdct_info)); |
10725 | 62 |
63 assert(N % 8 == 0); | |
64 | |
65 mdct->N = N; | |
66 | |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
67 /* NOTE: For "small framelengths" in FIXED_POINT the coefficients need to be |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
68 * scaled by sqrt("(nearest power of 2) > N" / N) */ |
10725 | 69 |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
70 /* RE(mdct->sincos[k]) = scale*(real_t)(cos(2.0*M_PI*(k+1./8.) / (real_t)N)); |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
71 * IM(mdct->sincos[k]) = scale*(real_t)(sin(2.0*M_PI*(k+1./8.) / (real_t)N)); */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
72 /* scale is 1 for fixed point, sqrt(N) for floating point */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
73 switch (N) |
10725 | 74 { |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
75 case 2048: mdct->sincos = (complex_t*)mdct_tab_2048; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
76 case 256: mdct->sincos = (complex_t*)mdct_tab_256; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
77 #ifdef LD_DEC |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
78 case 1024: mdct->sincos = (complex_t*)mdct_tab_1024; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
79 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
80 #ifdef ALLOW_SMALL_FRAMELENGTH |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
81 case 1920: mdct->sincos = (complex_t*)mdct_tab_1920; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
82 case 240: mdct->sincos = (complex_t*)mdct_tab_240; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
83 #ifdef LD_DEC |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
84 case 960: mdct->sincos = (complex_t*)mdct_tab_960; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
85 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
86 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
87 #ifdef SSR_DEC |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
88 case 512: mdct->sincos = (complex_t*)mdct_tab_512; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
89 case 64: mdct->sincos = (complex_t*)mdct_tab_64; break; |
10989 | 90 #endif |
10725 | 91 } |
92 | |
93 /* initialise fft */ | |
94 mdct->cfft = cffti(N/4); | |
95 | |
12527 | 96 #ifdef PROFILE |
97 mdct->cycles = 0; | |
98 mdct->fft_cycles = 0; | |
99 #endif | |
100 | |
10725 | 101 return mdct; |
102 } | |
103 | |
104 void faad_mdct_end(mdct_info *mdct) | |
105 { | |
106 if (mdct != NULL) | |
107 { | |
12527 | 108 #ifdef PROFILE |
109 printf("MDCT[%.4d]: %I64d cycles\n", mdct->N, mdct->cycles); | |
110 printf("CFFT[%.4d]: %I64d cycles\n", mdct->N/4, mdct->fft_cycles); | |
111 #endif | |
112 | |
10725 | 113 cfftu(mdct->cfft); |
114 | |
12527 | 115 faad_free(mdct); |
10725 | 116 } |
117 } | |
118 | |
119 void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out) | |
120 { | |
121 uint16_t k; | |
122 | |
123 complex_t x; | |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
124 #ifdef ALLOW_SMALL_FRAMELENGTH |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
125 #ifdef FIXED_POINT |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
126 real_t scale, b_scale = 0; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
127 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
128 #endif |
12527 | 129 ALIGN complex_t Z1[512]; |
130 complex_t *sincos = mdct->sincos; | |
131 | |
132 uint16_t N = mdct->N; | |
133 uint16_t N2 = N >> 1; | |
134 uint16_t N4 = N >> 2; | |
135 uint16_t N8 = N >> 3; | |
136 | |
137 #ifdef PROFILE | |
138 int64_t count1, count2 = faad_get_ts(); | |
139 #endif | |
140 | |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
141 #ifdef ALLOW_SMALL_FRAMELENGTH |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
142 #ifdef FIXED_POINT |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
143 /* detect non-power of 2 */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
144 if (N & (N-1)) |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
145 { |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
146 /* adjust scale for non-power of 2 MDCT */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
147 /* 2048/1920 */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
148 b_scale = 1; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
149 scale = COEF_CONST(1.0666666666666667); |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
150 } |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
151 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
152 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
153 |
12527 | 154 /* pre-IFFT complex multiplication */ |
155 for (k = 0; k < N4; k++) | |
156 { | |
157 ComplexMult(&IM(Z1[k]), &RE(Z1[k]), | |
158 X_in[2*k], X_in[N2 - 1 - 2*k], RE(sincos[k]), IM(sincos[k])); | |
159 } | |
160 | |
161 #ifdef PROFILE | |
162 count1 = faad_get_ts(); | |
163 #endif | |
164 | |
165 /* complex IFFT, any non-scaling FFT can be used here */ | |
166 cfftb(mdct->cfft, Z1); | |
167 | |
168 #ifdef PROFILE | |
169 count1 = faad_get_ts() - count1; | |
170 #endif | |
171 | |
172 /* post-IFFT complex multiplication */ | |
173 for (k = 0; k < N4; k++) | |
174 { | |
175 RE(x) = RE(Z1[k]); | |
176 IM(x) = IM(Z1[k]); | |
177 ComplexMult(&IM(Z1[k]), &RE(Z1[k]), | |
178 IM(x), RE(x), RE(sincos[k]), IM(sincos[k])); | |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
179 |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
180 #ifdef ALLOW_SMALL_FRAMELENGTH |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
181 #ifdef FIXED_POINT |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
182 /* non-power of 2 MDCT scaling */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
183 if (b_scale) |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
184 { |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
185 RE(Z1[k]) = MUL_C(RE(Z1[k]), scale); |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
186 IM(Z1[k]) = MUL_C(IM(Z1[k]), scale); |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
187 } |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
188 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
189 #endif |
12527 | 190 } |
191 | |
192 /* reordering */ | |
193 for (k = 0; k < N8; k+=2) | |
194 { | |
195 X_out[ 2*k] = IM(Z1[N8 + k]); | |
196 X_out[ 2 + 2*k] = IM(Z1[N8 + 1 + k]); | |
197 | |
198 X_out[ 1 + 2*k] = -RE(Z1[N8 - 1 - k]); | |
199 X_out[ 3 + 2*k] = -RE(Z1[N8 - 2 - k]); | |
200 | |
201 X_out[N4 + 2*k] = RE(Z1[ k]); | |
202 X_out[N4 + + 2 + 2*k] = RE(Z1[ 1 + k]); | |
203 | |
204 X_out[N4 + 1 + 2*k] = -IM(Z1[N4 - 1 - k]); | |
205 X_out[N4 + 3 + 2*k] = -IM(Z1[N4 - 2 - k]); | |
206 | |
207 X_out[N2 + 2*k] = RE(Z1[N8 + k]); | |
208 X_out[N2 + + 2 + 2*k] = RE(Z1[N8 + 1 + k]); | |
209 | |
210 X_out[N2 + 1 + 2*k] = -IM(Z1[N8 - 1 - k]); | |
211 X_out[N2 + 3 + 2*k] = -IM(Z1[N8 - 2 - k]); | |
212 | |
213 X_out[N2 + N4 + 2*k] = -IM(Z1[ k]); | |
214 X_out[N2 + N4 + 2 + 2*k] = -IM(Z1[ 1 + k]); | |
215 | |
216 X_out[N2 + N4 + 1 + 2*k] = RE(Z1[N4 - 1 - k]); | |
217 X_out[N2 + N4 + 3 + 2*k] = RE(Z1[N4 - 2 - k]); | |
218 } | |
219 | |
220 #ifdef PROFILE | |
221 count2 = faad_get_ts() - count2; | |
222 mdct->fft_cycles += count1; | |
223 mdct->cycles += (count2 - count1); | |
224 #endif | |
225 } | |
226 | |
10725 | 227 #ifdef LTP_DEC |
228 void faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out) | |
229 { | |
230 uint16_t k; | |
231 | |
232 complex_t x; | |
12527 | 233 ALIGN complex_t Z1[512]; |
10725 | 234 complex_t *sincos = mdct->sincos; |
235 | |
236 uint16_t N = mdct->N; | |
237 uint16_t N2 = N >> 1; | |
238 uint16_t N4 = N >> 2; | |
239 uint16_t N8 = N >> 3; | |
240 | |
12527 | 241 #ifndef FIXED_POINT |
10725 | 242 real_t scale = REAL_CONST(N); |
12527 | 243 #else |
244 real_t scale = REAL_CONST(4.0/N); | |
245 #endif | |
10725 | 246 |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
247 #ifdef ALLOW_SMALL_FRAMELENGTH |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
248 #ifdef FIXED_POINT |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
249 /* detect non-power of 2 */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
250 if (N & (N-1)) |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
251 { |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
252 /* adjust scale for non-power of 2 MDCT */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
253 /* *= sqrt(2048/1920) */ |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
254 scale = MUL_C(scale, COEF_CONST(1.0327955589886444)); |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
255 } |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
256 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
257 #endif |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
258 |
10725 | 259 /* pre-FFT complex multiplication */ |
260 for (k = 0; k < N8; k++) | |
261 { | |
262 uint16_t n = k << 1; | |
263 RE(x) = X_in[N - N4 - 1 - n] + X_in[N - N4 + n]; | |
264 IM(x) = X_in[ N4 + n] - X_in[ N4 - 1 - n]; | |
265 | |
12527 | 266 ComplexMult(&RE(Z1[k]), &IM(Z1[k]), |
267 RE(x), IM(x), RE(sincos[k]), IM(sincos[k])); | |
268 | |
269 RE(Z1[k]) = MUL_R(RE(Z1[k]), scale); | |
270 IM(Z1[k]) = MUL_R(IM(Z1[k]), scale); | |
10725 | 271 |
272 RE(x) = X_in[N2 - 1 - n] - X_in[ n]; | |
273 IM(x) = X_in[N2 + n] + X_in[N - 1 - n]; | |
274 | |
12527 | 275 ComplexMult(&RE(Z1[k + N8]), &IM(Z1[k + N8]), |
276 RE(x), IM(x), RE(sincos[k + N8]), IM(sincos[k + N8])); | |
277 | |
278 RE(Z1[k + N8]) = MUL_R(RE(Z1[k + N8]), scale); | |
279 IM(Z1[k + N8]) = MUL_R(IM(Z1[k + N8]), scale); | |
10725 | 280 } |
281 | |
10989 | 282 /* complex FFT, any non-scaling FFT can be used here */ |
10725 | 283 cfftf(mdct->cfft, Z1); |
284 | |
285 /* post-FFT complex multiplication */ | |
286 for (k = 0; k < N4; k++) | |
287 { | |
288 uint16_t n = k << 1; | |
12527 | 289 ComplexMult(&RE(x), &IM(x), |
290 RE(Z1[k]), IM(Z1[k]), RE(sincos[k]), IM(sincos[k])); | |
10725 | 291 |
12527 | 292 X_out[ n] = -RE(x); |
293 X_out[N2 - 1 - n] = IM(x); | |
294 X_out[N2 + n] = -IM(x); | |
295 X_out[N - 1 - n] = RE(x); | |
10725 | 296 } |
297 } | |
298 #endif |