2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
3 ** Copyright (C) 2003-2005 M. Bakker, Nero AG, http://www.nero.com
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.
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.
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.
19 ** Any non-GPL usage of this software or parts of this software is strictly
22 ** The "appropriate copyright message" mentioned in section 2c of the GPLv2
23 ** must read: "Code from FAAD2 is copyright (c) Nero AG, www.nero.com"
25 ** Commercial non-GPL licensing of this software is possible.
26 ** For more info contact Nero AG through Mpeg4AAClicense@nero.com.
28 ** $Id: mdct.c,v 1.47 2007/11/01 12:33:31 menno Exp $
32 * Fast (I)MDCT Implementation using (I)FFT ((Inverse) Fast Fourier Transform)
33 * and consists of three steps: pre-(I)FFT complex multiplication, complex
34 * (I)FFT, post-(I)FFT complex multiplication,
37 * P. Duhamel, Y. Mahieux, and J.P. Petit, "A Fast Algorithm for the
38 * Implementation of Filter Banks Based on 'Time Domain Aliasing
39 * Cancellation
\92," IEEE Proc. on ICASSP
\9191, 1991, pp. 2209-2212.
42 * As of April 6th 2002 completely rewritten.
43 * This (I)MDCT can now be used for any data size n, where n is divisible by 8.
62 mdct_info *faad_mdct_init(uint16_t N)
64 mdct_info *mdct = (mdct_info*)faad_malloc(sizeof(mdct_info));
70 /* NOTE: For "small framelengths" in FIXED_POINT the coefficients need to be
71 * scaled by sqrt("(nearest power of 2) > N" / N) */
73 /* RE(mdct->sincos[k]) = scale*(real_t)(cos(2.0*M_PI*(k+1./8.) / (real_t)N));
74 * IM(mdct->sincos[k]) = scale*(real_t)(sin(2.0*M_PI*(k+1./8.) / (real_t)N)); */
75 /* scale is 1 for fixed point, sqrt(N) for floating point */
78 case 2048: mdct->sincos = (complex_t*)mdct_tab_2048; break;
79 case 256: mdct->sincos = (complex_t*)mdct_tab_256; break;
81 case 1024: mdct->sincos = (complex_t*)mdct_tab_1024; break;
83 #ifdef ALLOW_SMALL_FRAMELENGTH
84 case 1920: mdct->sincos = (complex_t*)mdct_tab_1920; break;
85 case 240: mdct->sincos = (complex_t*)mdct_tab_240; break;
87 case 960: mdct->sincos = (complex_t*)mdct_tab_960; break;
91 case 512: mdct->sincos = (complex_t*)mdct_tab_512; break;
92 case 64: mdct->sincos = (complex_t*)mdct_tab_64; break;
97 mdct->cfft = cffti(N/4);
101 mdct->fft_cycles = 0;
107 void faad_mdct_end(mdct_info *mdct)
112 printf("MDCT[%.4d]: %I64d cycles\n", mdct->N, mdct->cycles);
113 printf("CFFT[%.4d]: %I64d cycles\n", mdct->N/4, mdct->fft_cycles);
122 void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
127 #ifdef ALLOW_SMALL_FRAMELENGTH
129 real_t scale, b_scale = 0;
132 ALIGN complex_t Z1[512];
133 complex_t *sincos = mdct->sincos;
135 uint16_t N = mdct->N;
136 uint16_t N2 = N >> 1;
137 uint16_t N4 = N >> 2;
138 uint16_t N8 = N >> 3;
141 int64_t count1, count2 = faad_get_ts();
144 #ifdef ALLOW_SMALL_FRAMELENGTH
146 /* detect non-power of 2 */
149 /* adjust scale for non-power of 2 MDCT */
152 scale = COEF_CONST(1.0666666666666667);
157 /* pre-IFFT complex multiplication */
158 for (k = 0; k < N4; k++)
160 ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
161 X_in[2*k], X_in[N2 - 1 - 2*k], RE(sincos[k]), IM(sincos[k]));
165 count1 = faad_get_ts();
168 /* complex IFFT, any non-scaling FFT can be used here */
169 cfftb(mdct->cfft, Z1);
172 count1 = faad_get_ts() - count1;
175 /* post-IFFT complex multiplication */
176 for (k = 0; k < N4; k++)
180 ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
181 IM(x), RE(x), RE(sincos[k]), IM(sincos[k]));
183 #ifdef ALLOW_SMALL_FRAMELENGTH
185 /* non-power of 2 MDCT scaling */
188 RE(Z1[k]) = MUL_C(RE(Z1[k]), scale);
189 IM(Z1[k]) = MUL_C(IM(Z1[k]), scale);
196 for (k = 0; k < N8; k+=2)
198 X_out[ 2*k] = IM(Z1[N8 + k]);
199 X_out[ 2 + 2*k] = IM(Z1[N8 + 1 + k]);
201 X_out[ 1 + 2*k] = -RE(Z1[N8 - 1 - k]);
202 X_out[ 3 + 2*k] = -RE(Z1[N8 - 2 - k]);
204 X_out[N4 + 2*k] = RE(Z1[ k]);
205 X_out[N4 + + 2 + 2*k] = RE(Z1[ 1 + k]);
207 X_out[N4 + 1 + 2*k] = -IM(Z1[N4 - 1 - k]);
208 X_out[N4 + 3 + 2*k] = -IM(Z1[N4 - 2 - k]);
210 X_out[N2 + 2*k] = RE(Z1[N8 + k]);
211 X_out[N2 + + 2 + 2*k] = RE(Z1[N8 + 1 + k]);
213 X_out[N2 + 1 + 2*k] = -IM(Z1[N8 - 1 - k]);
214 X_out[N2 + 3 + 2*k] = -IM(Z1[N8 - 2 - k]);
216 X_out[N2 + N4 + 2*k] = -IM(Z1[ k]);
217 X_out[N2 + N4 + 2 + 2*k] = -IM(Z1[ 1 + k]);
219 X_out[N2 + N4 + 1 + 2*k] = RE(Z1[N4 - 1 - k]);
220 X_out[N2 + N4 + 3 + 2*k] = RE(Z1[N4 - 2 - k]);
224 count2 = faad_get_ts() - count2;
225 mdct->fft_cycles += count1;
226 mdct->cycles += (count2 - count1);
231 void faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
236 ALIGN complex_t Z1[512];
237 complex_t *sincos = mdct->sincos;
239 uint16_t N = mdct->N;
240 uint16_t N2 = N >> 1;
241 uint16_t N4 = N >> 2;
242 uint16_t N8 = N >> 3;
245 real_t scale = REAL_CONST(N);
247 real_t scale = REAL_CONST(4.0/N);
250 #ifdef ALLOW_SMALL_FRAMELENGTH
252 /* detect non-power of 2 */
255 /* adjust scale for non-power of 2 MDCT */
256 /* *= sqrt(2048/1920) */
257 scale = MUL_C(scale, COEF_CONST(1.0327955589886444));
262 /* pre-FFT complex multiplication */
263 for (k = 0; k < N8; k++)
266 RE(x) = X_in[N - N4 - 1 - n] + X_in[N - N4 + n];
267 IM(x) = X_in[ N4 + n] - X_in[ N4 - 1 - n];
269 ComplexMult(&RE(Z1[k]), &IM(Z1[k]),
270 RE(x), IM(x), RE(sincos[k]), IM(sincos[k]));
272 RE(Z1[k]) = MUL_R(RE(Z1[k]), scale);
273 IM(Z1[k]) = MUL_R(IM(Z1[k]), scale);
275 RE(x) = X_in[N2 - 1 - n] - X_in[ n];
276 IM(x) = X_in[N2 + n] + X_in[N - 1 - n];
278 ComplexMult(&RE(Z1[k + N8]), &IM(Z1[k + N8]),
279 RE(x), IM(x), RE(sincos[k + N8]), IM(sincos[k + N8]));
281 RE(Z1[k + N8]) = MUL_R(RE(Z1[k + N8]), scale);
282 IM(Z1[k + N8]) = MUL_R(IM(Z1[k + N8]), scale);
285 /* complex FFT, any non-scaling FFT can be used here */
286 cfftf(mdct->cfft, Z1);
288 /* post-FFT complex multiplication */
289 for (k = 0; k < N4; k++)
292 ComplexMult(&RE(x), &IM(x),
293 RE(Z1[k]), IM(Z1[k]), RE(sincos[k]), IM(sincos[k]));
296 X_out[N2 - 1 - n] = IM(x);
297 X_out[N2 + n] = -IM(x);
298 X_out[N - 1 - n] = RE(x);