summaryrefslogtreecommitdiff
path: root/audio_codec/libfaad/mdct.c (plain)
blob: 71273dabf4f9fd45eebe2596919aed9472f2cb81
1/*
2** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
3** Copyright (C) 2003-2005 M. Bakker, Nero 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** 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"
24**
25** Commercial non-GPL licensing of this software is possible.
26** For more info contact Nero AG through Mpeg4AAClicense@nero.com.
27**
28** $Id: mdct.c,v 1.47 2007/11/01 12:33:31 menno Exp $
29**/
30
31/*
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,
35 *
36 * As described in:
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’," IEEE Proc. on ICASSP‘91, 1991, pp. 2209-2212.
40 *
41 *
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.
44 *
45 */
46#include <stdlib.h>
47#include "common.h"
48#include "structs.h"
49
50
51#ifdef _WIN32_WCE
52#define assert(x)
53#else
54#include <assert.h>
55#endif
56
57#include "cfft.h"
58#include "mdct.h"
59#include "mdct_tab.h"
60
61
62mdct_info *faad_mdct_init(uint16_t N)
63{
64 mdct_info *mdct = (mdct_info*)faad_malloc(sizeof(mdct_info));
65
66 assert(N % 8 == 0);
67
68 mdct->N = N;
69
70 /* NOTE: For "small framelengths" in FIXED_POINT the coefficients need to be
71 * scaled by sqrt("(nearest power of 2) > N" / N) */
72
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 */
76 switch (N) {
77 case 2048:
78 mdct->sincos = (complex_t*)mdct_tab_2048;
79 break;
80 case 256:
81 mdct->sincos = (complex_t*)mdct_tab_256;
82 break;
83#ifdef LD_DEC
84 case 1024:
85 mdct->sincos = (complex_t*)mdct_tab_1024;
86 break;
87#endif
88#ifdef ALLOW_SMALL_FRAMELENGTH
89 case 1920:
90 mdct->sincos = (complex_t*)mdct_tab_1920;
91 break;
92 case 240:
93 mdct->sincos = (complex_t*)mdct_tab_240;
94 break;
95#ifdef LD_DEC
96 case 960:
97 mdct->sincos = (complex_t*)mdct_tab_960;
98 break;
99#endif
100#endif
101#ifdef SSR_DEC
102 case 512:
103 mdct->sincos = (complex_t*)mdct_tab_512;
104 break;
105 case 64:
106 mdct->sincos = (complex_t*)mdct_tab_64;
107 break;
108#endif
109 }
110
111 /* initialise fft */
112 mdct->cfft = cffti(N / 4);
113
114#ifdef PROFILE
115 mdct->cycles = 0;
116 mdct->fft_cycles = 0;
117#endif
118
119 return mdct;
120}
121
122void faad_mdct_end(mdct_info *mdct)
123{
124 if (mdct != NULL) {
125#ifdef PROFILE
126 printf("MDCT[%.4d]: %I64d cycles\n", mdct->N, mdct->cycles);
127 printf("CFFT[%.4d]: %I64d cycles\n", mdct->N / 4, mdct->fft_cycles);
128#endif
129
130 cfftu(mdct->cfft);
131
132 faad_free(mdct);
133 }
134}
135
136void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
137{
138 uint16_t k;
139
140 complex_t x;
141#ifdef ALLOW_SMALL_FRAMELENGTH
142#ifdef FIXED_POINT
143 real_t scale, b_scale = 0;
144#endif
145#endif
146 ALIGN complex_t Z1[512];
147 complex_t *sincos = mdct->sincos;
148
149 uint16_t N = mdct->N;
150 uint16_t N2 = N >> 1;
151 uint16_t N4 = N >> 2;
152 uint16_t N8 = N >> 3;
153
154#ifdef PROFILE
155 int64_t count1, count2 = faad_get_ts();
156#endif
157
158#ifdef ALLOW_SMALL_FRAMELENGTH
159#ifdef FIXED_POINT
160 /* detect non-power of 2 */
161 if (N & (N - 1)) {
162 /* adjust scale for non-power of 2 MDCT */
163 /* 2048/1920 */
164 b_scale = 1;
165 scale = COEF_CONST(1.0666666666666667);
166 }
167#endif
168#endif
169
170 /* pre-IFFT complex multiplication */
171 for (k = 0; k < N4; k++) {
172 ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
173 X_in[2 * k], X_in[N2 - 1 - 2 * k], RE(sincos[k]), IM(sincos[k]));
174 }
175
176#ifdef PROFILE
177 count1 = faad_get_ts();
178#endif
179
180 /* complex IFFT, any non-scaling FFT can be used here */
181 cfftb(mdct->cfft, Z1);
182
183#ifdef PROFILE
184 count1 = faad_get_ts() - count1;
185#endif
186
187 /* post-IFFT complex multiplication */
188 for (k = 0; k < N4; k++) {
189 RE(x) = RE(Z1[k]);
190 IM(x) = IM(Z1[k]);
191 ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
192 IM(x), RE(x), RE(sincos[k]), IM(sincos[k]));
193
194#ifdef ALLOW_SMALL_FRAMELENGTH
195#ifdef FIXED_POINT
196 /* non-power of 2 MDCT scaling */
197 if (b_scale) {
198 RE(Z1[k]) = MUL_C(RE(Z1[k]), scale);
199 IM(Z1[k]) = MUL_C(IM(Z1[k]), scale);
200 }
201#endif
202#endif
203 }
204
205 /* reordering */
206 for (k = 0; k < N8; k += 2) {
207 X_out[ 2 * k] = IM(Z1[N8 + k]);
208 X_out[ 2 + 2 * k] = IM(Z1[N8 + 1 + k]);
209
210 X_out[ 1 + 2 * k] = -RE(Z1[N8 - 1 - k]);
211 X_out[ 3 + 2 * k] = -RE(Z1[N8 - 2 - k]);
212
213 X_out[N4 + 2 * k] = RE(Z1[ k]);
214 X_out[N4 + + 2 + 2 * k] = RE(Z1[ 1 + k]);
215
216 X_out[N4 + 1 + 2 * k] = -IM(Z1[N4 - 1 - k]);
217 X_out[N4 + 3 + 2 * k] = -IM(Z1[N4 - 2 - k]);
218
219 X_out[N2 + 2 * k] = RE(Z1[N8 + k]);
220 X_out[N2 + + 2 + 2 * k] = RE(Z1[N8 + 1 + k]);
221
222 X_out[N2 + 1 + 2 * k] = -IM(Z1[N8 - 1 - k]);
223 X_out[N2 + 3 + 2 * k] = -IM(Z1[N8 - 2 - k]);
224
225 X_out[N2 + N4 + 2 * k] = -IM(Z1[ k]);
226 X_out[N2 + N4 + 2 + 2 * k] = -IM(Z1[ 1 + k]);
227
228 X_out[N2 + N4 + 1 + 2 * k] = RE(Z1[N4 - 1 - k]);
229 X_out[N2 + N4 + 3 + 2 * k] = RE(Z1[N4 - 2 - k]);
230 }
231
232#ifdef PROFILE
233 count2 = faad_get_ts() - count2;
234 mdct->fft_cycles += count1;
235 mdct->cycles += (count2 - count1);
236#endif
237}
238
239#ifdef LTP_DEC
240void faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
241{
242 uint16_t k;
243
244 complex_t x;
245 ALIGN complex_t Z1[512];
246 complex_t *sincos = mdct->sincos;
247
248 uint16_t N = mdct->N;
249 uint16_t N2 = N >> 1;
250 uint16_t N4 = N >> 2;
251 uint16_t N8 = N >> 3;
252
253#ifndef FIXED_POINT
254 real_t scale = REAL_CONST(N);
255#else
256 real_t scale = REAL_CONST(4.0 / N);
257#endif
258
259#ifdef ALLOW_SMALL_FRAMELENGTH
260#ifdef FIXED_POINT
261 /* detect non-power of 2 */
262 if (N & (N - 1)) {
263 /* adjust scale for non-power of 2 MDCT */
264 /* *= sqrt(2048/1920) */
265 scale = MUL_C(scale, COEF_CONST(1.0327955589886444));
266 }
267#endif
268#endif
269
270 /* pre-FFT complex multiplication */
271 for (k = 0; k < N8; k++) {
272 uint16_t n = k << 1;
273 RE(x) = X_in[N - N4 - 1 - n] + X_in[N - N4 + n];
274 IM(x) = X_in[ N4 + n] - X_in[ N4 - 1 - n];
275
276 ComplexMult(&RE(Z1[k]), &IM(Z1[k]),
277 RE(x), IM(x), RE(sincos[k]), IM(sincos[k]));
278
279 RE(Z1[k]) = MUL_R(RE(Z1[k]), scale);
280 IM(Z1[k]) = MUL_R(IM(Z1[k]), scale);
281
282 RE(x) = X_in[N2 - 1 - n] - X_in[ n];
283 IM(x) = X_in[N2 + n] + X_in[N - 1 - n];
284
285 ComplexMult(&RE(Z1[k + N8]), &IM(Z1[k + N8]),
286 RE(x), IM(x), RE(sincos[k + N8]), IM(sincos[k + N8]));
287
288 RE(Z1[k + N8]) = MUL_R(RE(Z1[k + N8]), scale);
289 IM(Z1[k + N8]) = MUL_R(IM(Z1[k + N8]), scale);
290 }
291
292 /* complex FFT, any non-scaling FFT can be used here */
293 cfftf(mdct->cfft, Z1);
294
295 /* post-FFT complex multiplication */
296 for (k = 0; k < N4; k++) {
297 uint16_t n = k << 1;
298 ComplexMult(&RE(x), &IM(x),
299 RE(Z1[k]), IM(Z1[k]), RE(sincos[k]), IM(sincos[k]));
300
301 X_out[ n] = -RE(x);
302 X_out[N2 - 1 - n] = IM(x);
303 X_out[N2 + n] = -IM(x);
304 X_out[N - 1 - n] = RE(x);
305 }
306}
307#endif
308
generated by cgit at 2024-05-18 19:53:02 (GMT)