summaryrefslogtreecommitdiff
path: root/core/multimedia/opieplayer/libmpeg3/audio/dct.c
Unidiff
Diffstat (limited to 'core/multimedia/opieplayer/libmpeg3/audio/dct.c') (more/less context) (ignore whitespace changes)
-rw-r--r--core/multimedia/opieplayer/libmpeg3/audio/dct.c1135
1 files changed, 1135 insertions, 0 deletions
diff --git a/core/multimedia/opieplayer/libmpeg3/audio/dct.c b/core/multimedia/opieplayer/libmpeg3/audio/dct.c
new file mode 100644
index 0000000..1fd52ce
--- a/dev/null
+++ b/core/multimedia/opieplayer/libmpeg3/audio/dct.c
@@ -0,0 +1,1135 @@
1/*
2 *
3 * This file is part of libmpeg3
4 *
5 * libmpeg3 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, or (at your option)
8 * any later version.
9 *
10 * libmpeg3 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 GNU Make; see the file COPYING. If not, write to
17 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
18 *
19 */
20
21/*
22 * Discrete Cosine Tansform (DCT) for subband synthesis
23 * optimized for machines with no auto-increment.
24 * The performance is highly compiler dependend. Maybe
25 * the dct64.c version for 'normal' processor may be faster
26 * even for Intel processors.
27 */
28
29#include "mpeg3audio.h"
30#include "../libmpeg3.h"
31#include "../mpeg3protos.h"
32#include "tables.h"
33
34#include <math.h>
35
36int mpeg3audio_dct64_1(mpeg3_real_t *out0, mpeg3_real_t *out1, mpeg3_real_t *b1, mpeg3_real_t *b2, mpeg3_real_t *samples)
37{
38 register mpeg3_real_t *costab = mpeg3_pnts[0];
39
40 b1[0x00] = samples[0x00] + samples[0x1F];
41 b1[0x01] = samples[0x01] + samples[0x1E];
42 b1[0x1F] = (samples[0x00] - samples[0x1F]) * costab[0x0];
43 b1[0x1E] = (samples[0x01] - samples[0x1E]) * costab[0x1];
44
45 b1[0x02] = samples[0x02] + samples[0x1D];
46 b1[0x03] = samples[0x03] + samples[0x1C];
47 b1[0x1D] = (samples[0x02] - samples[0x1D]) * costab[0x2];
48 b1[0x1C] = (samples[0x03] - samples[0x1C]) * costab[0x3];
49
50 b1[0x04] = samples[0x04] + samples[0x1B];
51 b1[0x05] = samples[0x05] + samples[0x1A];
52 b1[0x1B] = (samples[0x04] - samples[0x1B]) * costab[0x4];
53 b1[0x1A] = (samples[0x05] - samples[0x1A]) * costab[0x5];
54
55 b1[0x06] = samples[0x06] + samples[0x19];
56 b1[0x07] = samples[0x07] + samples[0x18];
57 b1[0x19] = (samples[0x06] - samples[0x19]) * costab[0x6];
58 b1[0x18] = (samples[0x07] - samples[0x18]) * costab[0x7];
59
60 b1[0x08] = samples[0x08] + samples[0x17];
61 b1[0x09] = samples[0x09] + samples[0x16];
62 b1[0x17] = (samples[0x08] - samples[0x17]) * costab[0x8];
63 b1[0x16] = (samples[0x09] - samples[0x16]) * costab[0x9];
64
65 b1[0x0A] = samples[0x0A] + samples[0x15];
66 b1[0x0B] = samples[0x0B] + samples[0x14];
67 b1[0x15] = (samples[0x0A] - samples[0x15]) * costab[0xA];
68 b1[0x14] = (samples[0x0B] - samples[0x14]) * costab[0xB];
69
70 b1[0x0C] = samples[0x0C] + samples[0x13];
71 b1[0x0D] = samples[0x0D] + samples[0x12];
72 b1[0x13] = (samples[0x0C] - samples[0x13]) * costab[0xC];
73 b1[0x12] = (samples[0x0D] - samples[0x12]) * costab[0xD];
74
75 b1[0x0E] = samples[0x0E] + samples[0x11];
76 b1[0x0F] = samples[0x0F] + samples[0x10];
77 b1[0x11] = (samples[0x0E] - samples[0x11]) * costab[0xE];
78 b1[0x10] = (samples[0x0F] - samples[0x10]) * costab[0xF];
79
80 costab = mpeg3_pnts[1];
81
82 b2[0x00] = b1[0x00] + b1[0x0F];
83 b2[0x01] = b1[0x01] + b1[0x0E];
84 b2[0x0F] = (b1[0x00] - b1[0x0F]) * costab[0];
85 b2[0x0E] = (b1[0x01] - b1[0x0E]) * costab[1];
86
87 b2[0x02] = b1[0x02] + b1[0x0D];
88 b2[0x03] = b1[0x03] + b1[0x0C];
89 b2[0x0D] = (b1[0x02] - b1[0x0D]) * costab[2];
90 b2[0x0C] = (b1[0x03] - b1[0x0C]) * costab[3];
91
92 b2[0x04] = b1[0x04] + b1[0x0B];
93 b2[0x05] = b1[0x05] + b1[0x0A];
94 b2[0x0B] = (b1[0x04] - b1[0x0B]) * costab[4];
95 b2[0x0A] = (b1[0x05] - b1[0x0A]) * costab[5];
96
97 b2[0x06] = b1[0x06] + b1[0x09];
98 b2[0x07] = b1[0x07] + b1[0x08];
99 b2[0x09] = (b1[0x06] - b1[0x09]) * costab[6];
100 b2[0x08] = (b1[0x07] - b1[0x08]) * costab[7];
101
102 /* */
103
104 b2[0x10] = b1[0x10] + b1[0x1F];
105 b2[0x11] = b1[0x11] + b1[0x1E];
106 b2[0x1F] = (b1[0x1F] - b1[0x10]) * costab[0];
107 b2[0x1E] = (b1[0x1E] - b1[0x11]) * costab[1];
108
109 b2[0x12] = b1[0x12] + b1[0x1D];
110 b2[0x13] = b1[0x13] + b1[0x1C];
111 b2[0x1D] = (b1[0x1D] - b1[0x12]) * costab[2];
112 b2[0x1C] = (b1[0x1C] - b1[0x13]) * costab[3];
113
114 b2[0x14] = b1[0x14] + b1[0x1B];
115 b2[0x15] = b1[0x15] + b1[0x1A];
116 b2[0x1B] = (b1[0x1B] - b1[0x14]) * costab[4];
117 b2[0x1A] = (b1[0x1A] - b1[0x15]) * costab[5];
118
119 b2[0x16] = b1[0x16] + b1[0x19];
120 b2[0x17] = b1[0x17] + b1[0x18];
121 b2[0x19] = (b1[0x19] - b1[0x16]) * costab[6];
122 b2[0x18] = (b1[0x18] - b1[0x17]) * costab[7];
123
124 costab = mpeg3_pnts[2];
125
126 b1[0x00] = b2[0x00] + b2[0x07];
127 b1[0x07] = (b2[0x00] - b2[0x07]) * costab[0];
128 b1[0x01] = b2[0x01] + b2[0x06];
129 b1[0x06] = (b2[0x01] - b2[0x06]) * costab[1];
130 b1[0x02] = b2[0x02] + b2[0x05];
131 b1[0x05] = (b2[0x02] - b2[0x05]) * costab[2];
132 b1[0x03] = b2[0x03] + b2[0x04];
133 b1[0x04] = (b2[0x03] - b2[0x04]) * costab[3];
134
135 b1[0x08] = b2[0x08] + b2[0x0F];
136 b1[0x0F] = (b2[0x0F] - b2[0x08]) * costab[0];
137 b1[0x09] = b2[0x09] + b2[0x0E];
138 b1[0x0E] = (b2[0x0E] - b2[0x09]) * costab[1];
139 b1[0x0A] = b2[0x0A] + b2[0x0D];
140 b1[0x0D] = (b2[0x0D] - b2[0x0A]) * costab[2];
141 b1[0x0B] = b2[0x0B] + b2[0x0C];
142 b1[0x0C] = (b2[0x0C] - b2[0x0B]) * costab[3];
143
144 b1[0x10] = b2[0x10] + b2[0x17];
145 b1[0x17] = (b2[0x10] - b2[0x17]) * costab[0];
146 b1[0x11] = b2[0x11] + b2[0x16];
147 b1[0x16] = (b2[0x11] - b2[0x16]) * costab[1];
148 b1[0x12] = b2[0x12] + b2[0x15];
149 b1[0x15] = (b2[0x12] - b2[0x15]) * costab[2];
150 b1[0x13] = b2[0x13] + b2[0x14];
151 b1[0x14] = (b2[0x13] - b2[0x14]) * costab[3];
152
153 b1[0x18] = b2[0x18] + b2[0x1F];
154 b1[0x1F] = (b2[0x1F] - b2[0x18]) * costab[0];
155 b1[0x19] = b2[0x19] + b2[0x1E];
156 b1[0x1E] = (b2[0x1E] - b2[0x19]) * costab[1];
157 b1[0x1A] = b2[0x1A] + b2[0x1D];
158 b1[0x1D] = (b2[0x1D] - b2[0x1A]) * costab[2];
159 b1[0x1B] = b2[0x1B] + b2[0x1C];
160 b1[0x1C] = (b2[0x1C] - b2[0x1B]) * costab[3];
161
162 {
163 register mpeg3_real_t const cos0 = mpeg3_pnts[3][0];
164 register mpeg3_real_t const cos1 = mpeg3_pnts[3][1];
165
166 b2[0x00] = b1[0x00] + b1[0x03];
167 b2[0x03] = (b1[0x00] - b1[0x03]) * cos0;
168 b2[0x01] = b1[0x01] + b1[0x02];
169 b2[0x02] = (b1[0x01] - b1[0x02]) * cos1;
170
171 b2[0x04] = b1[0x04] + b1[0x07];
172 b2[0x07] = (b1[0x07] - b1[0x04]) * cos0;
173 b2[0x05] = b1[0x05] + b1[0x06];
174 b2[0x06] = (b1[0x06] - b1[0x05]) * cos1;
175
176 b2[0x08] = b1[0x08] + b1[0x0B];
177 b2[0x0B] = (b1[0x08] - b1[0x0B]) * cos0;
178 b2[0x09] = b1[0x09] + b1[0x0A];
179 b2[0x0A] = (b1[0x09] - b1[0x0A]) * cos1;
180
181 b2[0x0C] = b1[0x0C] + b1[0x0F];
182 b2[0x0F] = (b1[0x0F] - b1[0x0C]) * cos0;
183 b2[0x0D] = b1[0x0D] + b1[0x0E];
184 b2[0x0E] = (b1[0x0E] - b1[0x0D]) * cos1;
185
186 b2[0x10] = b1[0x10] + b1[0x13];
187 b2[0x13] = (b1[0x10] - b1[0x13]) * cos0;
188 b2[0x11] = b1[0x11] + b1[0x12];
189 b2[0x12] = (b1[0x11] - b1[0x12]) * cos1;
190
191 b2[0x14] = b1[0x14] + b1[0x17];
192 b2[0x17] = (b1[0x17] - b1[0x14]) * cos0;
193 b2[0x15] = b1[0x15] + b1[0x16];
194 b2[0x16] = (b1[0x16] - b1[0x15]) * cos1;
195
196 b2[0x18] = b1[0x18] + b1[0x1B];
197 b2[0x1B] = (b1[0x18] - b1[0x1B]) * cos0;
198 b2[0x19] = b1[0x19] + b1[0x1A];
199 b2[0x1A] = (b1[0x19] - b1[0x1A]) * cos1;
200
201 b2[0x1C] = b1[0x1C] + b1[0x1F];
202 b2[0x1F] = (b1[0x1F] - b1[0x1C]) * cos0;
203 b2[0x1D] = b1[0x1D] + b1[0x1E];
204 b2[0x1E] = (b1[0x1E] - b1[0x1D]) * cos1;
205 }
206
207 {
208 register mpeg3_real_t const cos0 = mpeg3_pnts[4][0];
209
210 b1[0x00] = b2[0x00] + b2[0x01];
211 b1[0x01] = (b2[0x00] - b2[0x01]) * cos0;
212 b1[0x02] = b2[0x02] + b2[0x03];
213 b1[0x03] = (b2[0x03] - b2[0x02]) * cos0;
214 b1[0x02] += b1[0x03];
215
216 b1[0x04] = b2[0x04] + b2[0x05];
217 b1[0x05] = (b2[0x04] - b2[0x05]) * cos0;
218 b1[0x06] = b2[0x06] + b2[0x07];
219 b1[0x07] = (b2[0x07] - b2[0x06]) * cos0;
220 b1[0x06] += b1[0x07];
221 b1[0x04] += b1[0x06];
222 b1[0x06] += b1[0x05];
223 b1[0x05] += b1[0x07];
224
225 b1[0x08] = b2[0x08] + b2[0x09];
226 b1[0x09] = (b2[0x08] - b2[0x09]) * cos0;
227 b1[0x0A] = b2[0x0A] + b2[0x0B];
228 b1[0x0B] = (b2[0x0B] - b2[0x0A]) * cos0;
229 b1[0x0A] += b1[0x0B];
230
231 b1[0x0C] = b2[0x0C] + b2[0x0D];
232 b1[0x0D] = (b2[0x0C] - b2[0x0D]) * cos0;
233 b1[0x0E] = b2[0x0E] + b2[0x0F];
234 b1[0x0F] = (b2[0x0F] - b2[0x0E]) * cos0;
235 b1[0x0E] += b1[0x0F];
236 b1[0x0C] += b1[0x0E];
237 b1[0x0E] += b1[0x0D];
238 b1[0x0D] += b1[0x0F];
239
240 b1[0x10] = b2[0x10] + b2[0x11];
241 b1[0x11] = (b2[0x10] - b2[0x11]) * cos0;
242 b1[0x12] = b2[0x12] + b2[0x13];
243 b1[0x13] = (b2[0x13] - b2[0x12]) * cos0;
244 b1[0x12] += b1[0x13];
245
246 b1[0x14] = b2[0x14] + b2[0x15];
247 b1[0x15] = (b2[0x14] - b2[0x15]) * cos0;
248 b1[0x16] = b2[0x16] + b2[0x17];
249 b1[0x17] = (b2[0x17] - b2[0x16]) * cos0;
250 b1[0x16] += b1[0x17];
251 b1[0x14] += b1[0x16];
252 b1[0x16] += b1[0x15];
253 b1[0x15] += b1[0x17];
254
255 b1[0x18] = b2[0x18] + b2[0x19];
256 b1[0x19] = (b2[0x18] - b2[0x19]) * cos0;
257 b1[0x1A] = b2[0x1A] + b2[0x1B];
258 b1[0x1B] = (b2[0x1B] - b2[0x1A]) * cos0;
259 b1[0x1A] += b1[0x1B];
260
261 b1[0x1C] = b2[0x1C] + b2[0x1D];
262 b1[0x1D] = (b2[0x1C] - b2[0x1D]) * cos0;
263 b1[0x1E] = b2[0x1E] + b2[0x1F];
264 b1[0x1F] = (b2[0x1F] - b2[0x1E]) * cos0;
265 b1[0x1E] += b1[0x1F];
266 b1[0x1C] += b1[0x1E];
267 b1[0x1E] += b1[0x1D];
268 b1[0x1D] += b1[0x1F];
269 }
270
271 out0[0x10*16] = b1[0x00];
272 out0[0x10*12] = b1[0x04];
273 out0[0x10* 8] = b1[0x02];
274 out0[0x10* 4] = b1[0x06];
275 out0[0x10* 0] = b1[0x01];
276 out1[0x10* 0] = b1[0x01];
277 out1[0x10* 4] = b1[0x05];
278 out1[0x10* 8] = b1[0x03];
279 out1[0x10*12] = b1[0x07];
280
281 out0[0x10*14] = b1[0x08] + b1[0x0C];
282 out0[0x10*10] = b1[0x0C] + b1[0x0a];
283 out0[0x10* 6] = b1[0x0A] + b1[0x0E];
284 out0[0x10* 2] = b1[0x0E] + b1[0x09];
285 out1[0x10* 2] = b1[0x09] + b1[0x0D];
286 out1[0x10* 6] = b1[0x0D] + b1[0x0B];
287 out1[0x10*10] = b1[0x0B] + b1[0x0F];
288 out1[0x10*14] = b1[0x0F];
289
290 {
291 register mpeg3_real_t tmp;
292 tmp = b1[0x18] + b1[0x1C];
293 out0[0x10*15] = tmp + b1[0x10];
294 out0[0x10*13] = tmp + b1[0x14];
295 tmp = b1[0x1C] + b1[0x1A];
296 out0[0x10*11] = tmp + b1[0x14];
297 out0[0x10* 9] = tmp + b1[0x12];
298 tmp = b1[0x1A] + b1[0x1E];
299 out0[0x10* 7] = tmp + b1[0x12];
300 out0[0x10* 5] = tmp + b1[0x16];
301 tmp = b1[0x1E] + b1[0x19];
302 out0[0x10* 3] = tmp + b1[0x16];
303 out0[0x10* 1] = tmp + b1[0x11];
304 tmp = b1[0x19] + b1[0x1D];
305 out1[0x10* 1] = tmp + b1[0x11];
306 out1[0x10* 3] = tmp + b1[0x15];
307 tmp = b1[0x1D] + b1[0x1B];
308 out1[0x10* 5] = tmp + b1[0x15];
309 out1[0x10* 7] = tmp + b1[0x13];
310 tmp = b1[0x1B] + b1[0x1F];
311 out1[0x10* 9] = tmp + b1[0x13];
312 out1[0x10*11] = tmp + b1[0x17];
313 out1[0x10*13] = b1[0x17] + b1[0x1F];
314 out1[0x10*15] = b1[0x1F];
315 }
316 return 0;
317}
318
319/*
320 * the call via dct64 is a trick to force GCC to use
321 * (new) registers for the b1,b2 pointer to the bufs[xx] field
322 */
323int mpeg3audio_dct64(mpeg3_real_t *a, mpeg3_real_t *b, mpeg3_real_t *c)
324{
325 mpeg3_real_t bufs[0x40];
326 return mpeg3audio_dct64_1(a, b, bufs, bufs + 0x20, c);
327}
328
329/*//////////////////////////////////////////////////////////////// */
330/* */
331/* 9 Point Inverse Discrete Cosine Transform */
332/* */
333/* This piece of code is Copyright 1997 Mikko Tommila and is freely usable */
334/* by anybody. The algorithm itself is of course in the public domain. */
335/* */
336/* Again derived heuristically from the 9-point WFTA. */
337/* */
338/* The algorithm is optimized (?) for speed, not for small rounding errors or */
339/* good readability. */
340/* */
341/* 36 additions, 11 multiplications */
342/* */
343/* Again this is very likely sub-optimal. */
344/* */
345/* The code is optimized to use a minimum number of temporary variables, */
346/* so it should compile quite well even on 8-register Intel x86 processors. */
347/* This makes the code quite obfuscated and very difficult to understand. */
348/* */
349/* References: */
350/* [1] S. Winograd: "On Computing the Discrete Fourier Transform", */
351/* Mathematics of Computation, Volume 32, Number 141, January 1978, */
352/* Pages 175-199 */
353
354
355/*------------------------------------------------------------------*/
356/* */
357/* Function: Calculation of the inverse MDCT */
358/* */
359/*------------------------------------------------------------------*/
360
361int mpeg3audio_dct36(mpeg3_real_t *inbuf, mpeg3_real_t *o1, mpeg3_real_t *o2, mpeg3_real_t *wintab, mpeg3_real_t *tsbuf)
362{
363 mpeg3_real_t tmp[18];
364
365 {
366 register mpeg3_real_t *in = inbuf;
367
368 in[17]+=in[16]; in[16]+=in[15]; in[15]+=in[14];
369 in[14]+=in[13]; in[13]+=in[12]; in[12]+=in[11];
370 in[11]+=in[10]; in[10]+=in[9]; in[9] +=in[8];
371 in[8] +=in[7]; in[7] +=in[6]; in[6] +=in[5];
372 in[5] +=in[4]; in[4] +=in[3]; in[3] +=in[2];
373 in[2] +=in[1]; in[1] +=in[0];
374
375 in[17]+=in[15]; in[15]+=in[13]; in[13]+=in[11]; in[11]+=in[9];
376 in[9] +=in[7]; in[7] +=in[5]; in[5] +=in[3]; in[3] +=in[1];
377
378
379 {
380 mpeg3_real_t t3;
381 {
382 mpeg3_real_t t0, t1, t2;
383
384 t0 = mpeg3_COS6_2 * (in[8] + in[16] - in[4]);
385 t1 = mpeg3_COS6_2 * in[12];
386
387 t3 = in[0];
388 t2 = t3 - t1 - t1;
389 tmp[1] = tmp[7] = t2 - t0;
390 tmp[4] = t2 + t0 + t0;
391 t3 += t1;
392
393 t2 = mpeg3_COS6_1 * (in[10] + in[14] - in[2]);
394 tmp[1] -= t2;
395 tmp[7] += t2;
396 }
397 {
398 mpeg3_real_t t0, t1, t2;
399
400 t0 = mpeg3_cos9[0] * (in[4] + in[8] );
401 t1 = mpeg3_cos9[1] * (in[8] - in[16]);
402 t2 = mpeg3_cos9[2] * (in[4] + in[16]);
403
404 tmp[2] = tmp[6] = t3 - t0 - t2;
405 tmp[0] = tmp[8] = t3 + t0 + t1;
406 tmp[3] = tmp[5] = t3 - t1 + t2;
407 }
408 }
409 {
410 mpeg3_real_t t1, t2, t3;
411
412 t1 = mpeg3_cos18[0] * (in[2] + in[10]);
413 t2 = mpeg3_cos18[1] * (in[10] - in[14]);
414 t3 = mpeg3_COS6_1 * in[6];
415
416 {
417 mpeg3_real_t t0 = t1 + t2 + t3;
418 tmp[0] += t0;
419 tmp[8] -= t0;
420 }
421
422 t2 -= t3;
423 t1 -= t3;
424
425 t3 = mpeg3_cos18[2] * (in[2] + in[14]);
426
427 t1 += t3;
428 tmp[3] += t1;
429 tmp[5] -= t1;
430
431 t2 -= t3;
432 tmp[2] += t2;
433 tmp[6] -= t2;
434 }
435
436
437 {
438 mpeg3_real_t t0, t1, t2, t3, t4, t5, t6, t7;
439
440 t1 = mpeg3_COS6_2 * in[13];
441 t2 = mpeg3_COS6_2 * (in[9] + in[17] - in[5]);
442
443 t3 = in[1] + t1;
444 t4 = in[1] - t1 - t1;
445 t5 = t4 - t2;
446
447 t0 = mpeg3_cos9[0] * (in[5] + in[9]);
448 t1 = mpeg3_cos9[1] * (in[9] - in[17]);
449
450 tmp[13] = (t4 + t2 + t2) * mpeg3_tfcos36[17-13];
451 t2 = mpeg3_cos9[2] * (in[5] + in[17]);
452
453 t6 = t3 - t0 - t2;
454 t0 += t3 + t1;
455 t3 += t2 - t1;
456
457 t2 = mpeg3_cos18[0] * (in[3] + in[11]);
458 t4 = mpeg3_cos18[1] * (in[11] - in[15]);
459 t7 = mpeg3_COS6_1 * in[7];
460
461 t1 = t2 + t4 + t7;
462 tmp[17] = (t0 + t1) * mpeg3_tfcos36[17-17];
463 tmp[9] = (t0 - t1) * mpeg3_tfcos36[17-9];
464 t1 = mpeg3_cos18[2] * (in[3] + in[15]);
465 t2 += t1 - t7;
466
467 tmp[14] = (t3 + t2) * mpeg3_tfcos36[17-14];
468 t0 = mpeg3_COS6_1 * (in[11] + in[15] - in[3]);
469 tmp[12] = (t3 - t2) * mpeg3_tfcos36[17-12];
470
471 t4 -= t1 + t7;
472
473 tmp[16] = (t5 - t0) * mpeg3_tfcos36[17-16];
474 tmp[10] = (t5 + t0) * mpeg3_tfcos36[17-10];
475 tmp[15] = (t6 + t4) * mpeg3_tfcos36[17-15];
476 tmp[11] = (t6 - t4) * mpeg3_tfcos36[17-11];
477 }
478
479#define MACRO(v) \
480 { \
481 mpeg3_real_t tmpval; \
482 tmpval = tmp[(v)] + tmp[17-(v)]; \
483 out2[9+(v)] = tmpval * w[27+(v)]; \
484 out2[8-(v)] = tmpval * w[26-(v)]; \
485 tmpval = tmp[(v)] - tmp[17-(v)]; \
486 ts[SBLIMIT*(8-(v))] = out1[8-(v)] + tmpval * w[8-(v)]; \
487 ts[SBLIMIT*(9+(v))] = out1[9+(v)] + tmpval * w[9+(v)]; \
488 }
489
490 {
491 register mpeg3_real_t *out2 = o2;
492 register mpeg3_real_t *w = wintab;
493 register mpeg3_real_t *out1 = o1;
494 register mpeg3_real_t *ts = tsbuf;
495
496 MACRO(0);
497 MACRO(1);
498 MACRO(2);
499 MACRO(3);
500 MACRO(4);
501 MACRO(5);
502 MACRO(6);
503 MACRO(7);
504 MACRO(8);
505 }
506 }
507 return 0;
508}
509
510/*
511 * new DCT12
512 */
513int mpeg3audio_dct12(mpeg3_real_t *in,mpeg3_real_t *rawout1,mpeg3_real_t *rawout2,register mpeg3_real_t *wi,register mpeg3_real_t *ts)
514{
515#define DCT12_PART1 \
516 in5 = in[5*3]; \
517 in5 += (in4 = in[4*3]); \
518 in4 += (in3 = in[3*3]); \
519 in3 += (in2 = in[2*3]); \
520 in2 += (in1 = in[1*3]); \
521 in1 += (in0 = in[0*3]); \
522 \
523 in5 += in3; in3 += in1; \
524 \
525 in2 *= mpeg3_COS6_1; \
526 in3 *= mpeg3_COS6_1; \
527
528#define DCT12_PART2 \
529 in0 += in4 * mpeg3_COS6_2; \
530 \
531 in4 = in0 + in2; \
532 in0 -= in2; \
533 \
534 in1 += in5 * mpeg3_COS6_2; \
535 \
536 in5 = (in1 + in3) * mpeg3_tfcos12[0]; \
537 in1 = (in1 - in3) * mpeg3_tfcos12[2]; \
538 \
539 in3 = in4 + in5; \
540 in4 -= in5; \
541 \
542 in2 = in0 + in1; \
543 in0 -= in1;
544
545
546 {
547 mpeg3_real_t in0,in1,in2,in3,in4,in5;
548 register mpeg3_real_t *out1 = rawout1;
549 ts[SBLIMIT*0] = out1[0]; ts[SBLIMIT*1] = out1[1]; ts[SBLIMIT*2] = out1[2];
550 ts[SBLIMIT*3] = out1[3]; ts[SBLIMIT*4] = out1[4]; ts[SBLIMIT*5] = out1[5];
551
552 DCT12_PART1
553
554 {
555 mpeg3_real_t tmp0,tmp1 = (in0 - in4);
556 {
557 mpeg3_real_t tmp2 = (in1 - in5) * mpeg3_tfcos12[1];
558 tmp0 = tmp1 + tmp2;
559 tmp1 -= tmp2;
560 }
561 ts[(17-1)*SBLIMIT] = out1[17-1] + tmp0 * wi[11-1];
562 ts[(12+1)*SBLIMIT] = out1[12+1] + tmp0 * wi[6+1];
563 ts[(6 +1)*SBLIMIT] = out1[6 +1] + tmp1 * wi[1];
564 ts[(11-1)*SBLIMIT] = out1[11-1] + tmp1 * wi[5-1];
565 }
566
567 DCT12_PART2
568
569 ts[(17-0)*SBLIMIT] = out1[17-0] + in2 * wi[11-0];
570 ts[(12+0)*SBLIMIT] = out1[12+0] + in2 * wi[6+0];
571 ts[(12+2)*SBLIMIT] = out1[12+2] + in3 * wi[6+2];
572 ts[(17-2)*SBLIMIT] = out1[17-2] + in3 * wi[11-2];
573
574 ts[(6+0)*SBLIMIT] = out1[6+0] + in0 * wi[0];
575 ts[(11-0)*SBLIMIT] = out1[11-0] + in0 * wi[5-0];
576 ts[(6+2)*SBLIMIT] = out1[6+2] + in4 * wi[2];
577 ts[(11-2)*SBLIMIT] = out1[11-2] + in4 * wi[5-2];
578 }
579
580 in++;
581
582 {
583 mpeg3_real_t in0,in1,in2,in3,in4,in5;
584 register mpeg3_real_t *out2 = rawout2;
585
586 DCT12_PART1
587
588 {
589 mpeg3_real_t tmp0,tmp1 = (in0 - in4);
590 {
591 mpeg3_real_t tmp2 = (in1 - in5) * mpeg3_tfcos12[1];
592 tmp0 = tmp1 + tmp2;
593 tmp1 -= tmp2;
594 }
595 out2[5-1] = tmp0 * wi[11-1];
596 out2[0+1] = tmp0 * wi[6+1];
597 ts[(12+1)*SBLIMIT] += tmp1 * wi[1];
598 ts[(17-1)*SBLIMIT] += tmp1 * wi[5-1];
599 }
600
601 DCT12_PART2
602
603 out2[5-0] = in2 * wi[11-0];
604 out2[0+0] = in2 * wi[6+0];
605 out2[0+2] = in3 * wi[6+2];
606 out2[5-2] = in3 * wi[11-2];
607
608 ts[(12+0)*SBLIMIT] += in0 * wi[0];
609 ts[(17-0)*SBLIMIT] += in0 * wi[5-0];
610 ts[(12+2)*SBLIMIT] += in4 * wi[2];
611 ts[(17-2)*SBLIMIT] += in4 * wi[5-2];
612 }
613
614 in++;
615
616 {
617 mpeg3_real_t in0,in1,in2,in3,in4,in5;
618 register mpeg3_real_t *out2 = rawout2;
619 out2[12]=out2[13]=out2[14]=out2[15]=out2[16]=out2[17]=0.0;
620
621 DCT12_PART1
622
623 {
624 mpeg3_real_t tmp0,tmp1 = (in0 - in4);
625 {
626 mpeg3_real_t tmp2 = (in1 - in5) * mpeg3_tfcos12[1];
627 tmp0 = tmp1 + tmp2;
628 tmp1 -= tmp2;
629 }
630 out2[11-1] = tmp0 * wi[11-1];
631 out2[6 +1] = tmp0 * wi[6+1];
632 out2[0+1] += tmp1 * wi[1];
633 out2[5-1] += tmp1 * wi[5-1];
634 }
635
636 DCT12_PART2
637
638 out2[11-0] = in2 * wi[11-0];
639 out2[6 +0] = in2 * wi[6+0];
640 out2[6 +2] = in3 * wi[6+2];
641 out2[11-2] = in3 * wi[11-2];
642
643 out2[0+0] += in0 * wi[0];
644 out2[5-0] += in0 * wi[5-0];
645 out2[0+2] += in4 * wi[2];
646 out2[5-2] += in4 * wi[5-2];
647 }
648 return 0;
649}
650
651/* AC3 IMDCT tables */
652
653/* Twiddle factors for IMDCT */
654#if !defined(USE_FIXED_POINT) || defined(PRINT_FIXED_POINT_TABLES)
655static mpeg3_real_t mpeg3_xcos1[AC3_N / 4];
656static mpeg3_real_t mpeg3_xsin1[AC3_N / 4];
657static mpeg3_real_t mpeg3_xcos2[AC3_N / 8];
658static mpeg3_real_t mpeg3_xsin2[AC3_N / 8];
659#else
660#define USE_FP_TABLES
661#include "fptables.h"
662#endif
663
664/* 128 point bit-reverse LUT */
665static unsigned char mpeg3_bit_reverse_512[] =
666{
667 0x00, 0x40, 0x20, 0x60, 0x10, 0x50, 0x30, 0x70,
668 0x08, 0x48, 0x28, 0x68, 0x18, 0x58, 0x38, 0x78,
669 0x04, 0x44, 0x24, 0x64, 0x14, 0x54, 0x34, 0x74,
670 0x0c, 0x4c, 0x2c, 0x6c, 0x1c, 0x5c, 0x3c, 0x7c,
671 0x02, 0x42, 0x22, 0x62, 0x12, 0x52, 0x32, 0x72,
672 0x0a, 0x4a, 0x2a, 0x6a, 0x1a, 0x5a, 0x3a, 0x7a,
673 0x06, 0x46, 0x26, 0x66, 0x16, 0x56, 0x36, 0x76,
674 0x0e, 0x4e, 0x2e, 0x6e, 0x1e, 0x5e, 0x3e, 0x7e,
675 0x01, 0x41, 0x21, 0x61, 0x11, 0x51, 0x31, 0x71,
676 0x09, 0x49, 0x29, 0x69, 0x19, 0x59, 0x39, 0x79,
677 0x05, 0x45, 0x25, 0x65, 0x15, 0x55, 0x35, 0x75,
678 0x0d, 0x4d, 0x2d, 0x6d, 0x1d, 0x5d, 0x3d, 0x7d,
679 0x03, 0x43, 0x23, 0x63, 0x13, 0x53, 0x33, 0x73,
680 0x0b, 0x4b, 0x2b, 0x6b, 0x1b, 0x5b, 0x3b, 0x7b,
681 0x07, 0x47, 0x27, 0x67, 0x17, 0x57, 0x37, 0x77,
682 0x0f, 0x4f, 0x2f, 0x6f, 0x1f, 0x5f, 0x3f, 0x7f
683};
684
685static unsigned char mpeg3_bit_reverse_256[] =
686{
687 0x00, 0x20, 0x10, 0x30, 0x08, 0x28, 0x18, 0x38,
688 0x04, 0x24, 0x14, 0x34, 0x0c, 0x2c, 0x1c, 0x3c,
689 0x02, 0x22, 0x12, 0x32, 0x0a, 0x2a, 0x1a, 0x3a,
690 0x06, 0x26, 0x16, 0x36, 0x0e, 0x2e, 0x1e, 0x3e,
691 0x01, 0x21, 0x11, 0x31, 0x09, 0x29, 0x19, 0x39,
692 0x05, 0x25, 0x15, 0x35, 0x0d, 0x2d, 0x1d, 0x3d,
693 0x03, 0x23, 0x13, 0x33, 0x0b, 0x2b, 0x1b, 0x3b,
694 0x07, 0x27, 0x17, 0x37, 0x0f, 0x2f, 0x1f, 0x3f
695};
696
697/* Windowing function for Modified DCT - Thank you acroread */
698static mpeg3_real_t mpeg3_window[] =
699{
700 0.00014, 0.00024, 0.00037, 0.00051, 0.00067, 0.00086, 0.00107, 0.00130,
701 0.00157, 0.00187, 0.00220, 0.00256, 0.00297, 0.00341, 0.00390, 0.00443,
702 0.00501, 0.00564, 0.00632, 0.00706, 0.00785, 0.00871, 0.00962, 0.01061,
703 0.01166, 0.01279, 0.01399, 0.01526, 0.01662, 0.01806, 0.01959, 0.02121,
704 0.02292, 0.02472, 0.02662, 0.02863, 0.03073, 0.03294, 0.03527, 0.03770,
705 0.04025, 0.04292, 0.04571, 0.04862, 0.05165, 0.05481, 0.05810, 0.06153,
706 0.06508, 0.06878, 0.07261, 0.07658, 0.08069, 0.08495, 0.08935, 0.09389,
707 0.09859, 0.10343, 0.10842, 0.11356, 0.11885, 0.12429, 0.12988, 0.13563,
708 0.14152, 0.14757, 0.15376, 0.16011, 0.16661, 0.17325, 0.18005, 0.18699,
709 0.19407, 0.20130, 0.20867, 0.21618, 0.22382, 0.23161, 0.23952, 0.24757,
710 0.25574, 0.26404, 0.27246, 0.28100, 0.28965, 0.29841, 0.30729, 0.31626,
711 0.32533, 0.33450, 0.34376, 0.35311, 0.36253, 0.37204, 0.38161, 0.39126,
712 0.40096, 0.41072, 0.42054, 0.43040, 0.44030, 0.45023, 0.46020, 0.47019,
713 0.48020, 0.49022, 0.50025, 0.51028, 0.52031, 0.53033, 0.54033, 0.55031,
714 0.56026, 0.57019, 0.58007, 0.58991, 0.59970, 0.60944, 0.61912, 0.62873,
715 0.63827, 0.64774, 0.65713, 0.66643, 0.67564, 0.68476, 0.69377, 0.70269,
716 0.71150, 0.72019, 0.72877, 0.73723, 0.74557, 0.75378, 0.76186, 0.76981,
717 0.77762, 0.78530, 0.79283, 0.80022, 0.80747, 0.81457, 0.82151, 0.82831,
718 0.83496, 0.84145, 0.84779, 0.85398, 0.86001, 0.86588, 0.87160, 0.87716,
719 0.88257, 0.88782, 0.89291, 0.89785, 0.90264, 0.90728, 0.91176, 0.91610,
720 0.92028, 0.92432, 0.92822, 0.93197, 0.93558, 0.93906, 0.94240, 0.94560,
721 0.94867, 0.95162, 0.95444, 0.95713, 0.95971, 0.96217, 0.96451, 0.96674,
722 0.96887, 0.97089, 0.97281, 0.97463, 0.97635, 0.97799, 0.97953, 0.98099,
723 0.98236, 0.98366, 0.98488, 0.98602, 0.98710, 0.98811, 0.98905, 0.98994,
724 0.99076, 0.99153, 0.99225, 0.99291, 0.99353, 0.99411, 0.99464, 0.99513,
725 0.99558, 0.99600, 0.99639, 0.99674, 0.99706, 0.99736, 0.99763, 0.99788,
726 0.99811, 0.99831, 0.99850, 0.99867, 0.99882, 0.99895, 0.99908, 0.99919,
727 0.99929, 0.99938, 0.99946, 0.99953, 0.99959, 0.99965, 0.99969, 0.99974,
728 0.99978, 0.99981, 0.99984, 0.99986, 0.99988, 0.99990, 0.99992, 0.99993,
729 0.99994, 0.99995, 0.99996, 0.99997, 0.99998, 0.99998, 0.99998, 0.99999,
730 0.99999, 0.99999, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
731 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000
732};
733
734mpeg3_complex_t cmplx_mult(mpeg3_complex_t a, mpeg3_complex_t b)
735{
736 mpeg3_complex_t ret;
737
738 ret.real = a.real * b.real - a.imag * b.imag;
739 ret.imag = a.real * b.imag + a.imag * b.real;
740
741 return ret;
742}
743
744int mpeg3audio_imdct_init(mpeg3audio_t *audio)
745{
746 int i, k;
747 mpeg3_complex_t angle_step;
748 mpeg3_complex_t current_angle;
749
750/* Twiddle factors to turn IFFT into IMDCT */
751 for(i = 0; i < AC3_N / 4; i++)
752 {
753 mpeg3_xcos1[i] = -cos(2.0f * M_PI * (8 * i + 1 ) / ( 8 * AC3_N));
754 mpeg3_xsin1[i] = -sin(2.0f * M_PI * (8 * i + 1 ) / ( 8 * AC3_N));
755 }
756
757/* More twiddle factors to turn IFFT into IMDCT */
758 for(i = 0; i < AC3_N / 8; i++)
759 {
760 mpeg3_xcos2[i] = -cos(2.0f * M_PI * (8 * i + 1 ) / ( 4 * AC3_N));
761 mpeg3_xsin2[i] = -sin(2.0f * M_PI * (8 * i + 1 ) / ( 4 * AC3_N));
762 }
763
764/* Canonical twiddle factors for FFT */
765#if defined(USE_FIXED_POINT) && !defined(PRINT_FIXED_POINT_TABLES)
766 for(i = 0; i < 7; i++)
767 {
768 audio->ac3_w[i] = (mpeg3_complex_t*)ac3_w_fixedpoints[i];
769 }
770#else
771 audio->ac3_w[0] = audio->ac3_w_1;
772 audio->ac3_w[1] = audio->ac3_w_2;
773 audio->ac3_w[2] = audio->ac3_w_4;
774 audio->ac3_w[3] = audio->ac3_w_8;
775 audio->ac3_w[4] = audio->ac3_w_16;
776 audio->ac3_w[5] = audio->ac3_w_32;
777 audio->ac3_w[6] = audio->ac3_w_64;
778
779 for(i = 0; i < 7; i++)
780 {
781 angle_step.real = cos(-2.0f * M_PI / (1 << (i + 1)));
782 angle_step.imag = sin(-2.0f * M_PI / (1 << (i + 1)));
783
784 current_angle.real = 1.0f;
785 current_angle.imag = 0.0f;
786
787 for (k = 0; k < 1 << i; k++)
788 {
789 audio->ac3_w[i][k] = current_angle;
790 current_angle = cmplx_mult(current_angle, angle_step);
791 }
792 }
793
794#ifdef PRINT_FIXED_POINT_TABLES
795 printf("#ifdef USE_FP_TABLES\n");
796 printf("static long mpeg3_xcos1_fixedpoints[] = {");
797 for(i = 0; i < AC3_N / 4; i++) {
798 printf("%c0x%08x,", i%8?' ':'\n', mpeg3_xcos1[i].fixedPoint());
799 }
800 printf("\n};\nstatic mpeg3_real_t *mpeg3_xcos1 = \n"
801 "(mpeg3_real_t*)mpeg3_xcos1_fixedpoints;\n");
802
803 printf("static long mpeg3_xsin1_fixedpoints[] = {");
804 for(i = 0; i < AC3_N / 4; i++) {
805 printf("%c0x%08x,", i%8?' ':'\n', mpeg3_xsin1[i].fixedPoint());
806 }
807 printf("\n};\nstatic mpeg3_real_t *mpeg3_xsin1 = \n"
808 "(mpeg3_real_t*)mpeg3_xsin1_fixedpoints;\n");
809
810
811 printf("static long mpeg3_xcos2_fixedpoints[] = {");
812 for(i = 0; i < AC3_N / 4; i++) {
813 printf("%c0x%08x,", i%8?' ':'\n', mpeg3_xcos2[i].fixedPoint());
814 }
815 printf("\n};\nstatic mpeg3_real_t *mpeg3_xcos2 = \n"
816 "(mpeg3_real_t*)mpeg3_xcos2_fixedpoints;\n");
817
818 printf("static long mpeg3_xsin2_fixedpoints[] = {");
819 for(i = 0; i < AC3_N / 4; i++) {
820 printf("%c0x%08x,", i%8?' ':'\n', mpeg3_xsin2[i].fixedPoint());
821 }
822 printf("\n};\nstatic mpeg3_real_t *mpeg3_xsin2 = \n"
823 "(mpeg3_real_t*)mpeg3_xsin2_fixedpoints;\n");
824
825
826 printf("typedef struct { long r, i; } fixed_cmplx;\n");
827 for(i = 0; i < 7; i++)
828 {
829 printf("fixed_cmplx ac3_w_d%d[] = { ", 1<<i);
830 for (k = 0; k < 1 << i; k++)
831 {
832 printf("%s{ 0x%08x, 0x%08x },", k%4?" ":"\n ",
833 audio->ac3_w[i][k].real.fixedPoint(),
834 audio->ac3_w[i][k].imag.fixedPoint());
835 }
836 printf("};\n");
837 }
838
839 printf("fixed_cmplx *ac3_w_fixedpoints[] = {\n");
840 for(i = 0; i < 7; i++)
841 {
842 printf("ac3_w_d%d, ", 1<<i);
843 }
844 printf("};\n");
845 printf("#endif\n");
846#endif
847#endif
848
849 return 0;
850}
851
852
853inline void swap_cmplx(mpeg3_complex_t *a, mpeg3_complex_t *b)
854{
855 mpeg3_complex_t tmp;
856
857 tmp = *a;
858 *a = *b;
859 *b = tmp;
860}
861
862void mpeg3audio_ac3_imdct_do_512(mpeg3audio_t *audio,
863 mpeg3_real_t data[],
864 mpeg3_real_t *y,
865 int step,
866 mpeg3_real_t *delay)
867{
868 int i, k;
869 int p, q;
870 int m;
871 int two_m;
872 int two_m_plus_one;
873
874 mpeg3_real_t tmp_a_i;
875 mpeg3_real_t tmp_a_r;
876 mpeg3_real_t tmp_b_i;
877 mpeg3_real_t tmp_b_r;
878
879 mpeg3_real_t *y_ptr;
880 mpeg3_real_t *delay_ptr;
881 mpeg3_real_t *window_ptr;
882 mpeg3_complex_t *buf = audio->ac3_imdct_buf;
883
884/* Pre IFFT complex multiply plus IFFT cmplx conjugate */
885 for(i = 0; i < AC3_N / 4; i++)
886 {
887 buf[i].real = (data[AC3_N / 2 - 2 * i - 1] * mpeg3_xcos1[i]) - (data[2 * i] * mpeg3_xsin1[i]);
888 buf[i].imag = -((data[2 * i] * mpeg3_xcos1[i]) + (data[AC3_N / 2 - 2 * i - 1] * mpeg3_xsin1[i]));
889 }
890
891/* Bit reversed shuffling */
892 for(i = 0; i < AC3_N / 4; i++)
893 {
894 k = mpeg3_bit_reverse_512[i];
895 if(k < i)
896 swap_cmplx(&buf[i], &buf[k]);
897 }
898
899/* FFT Merge */
900 for(m = 0; m < 7; m++)
901 {
902 if(m)
903 two_m = (1 << m);
904 else
905 two_m = 1;
906
907 two_m_plus_one = (1 << (m + 1));
908
909 for(k = 0; k < two_m; k++)
910 {
911 for(i = 0; i < AC3_N / 4; i += two_m_plus_one)
912 {
913 p = k + i;
914 q = p + two_m;
915 tmp_a_r = buf[p].real;
916 tmp_a_i = buf[p].imag;
917 tmp_b_r = buf[q].real * audio->ac3_w[m][k].real - buf[q].imag * audio->ac3_w[m][k].imag;
918 tmp_b_i = buf[q].imag * audio->ac3_w[m][k].real + buf[q].real * audio->ac3_w[m][k].imag;
919 buf[p].real = tmp_a_r + tmp_b_r;
920 buf[p].imag = tmp_a_i + tmp_b_i;
921 buf[q].real = tmp_a_r - tmp_b_r;
922 buf[q].imag = tmp_a_i - tmp_b_i;
923 }
924 }
925 }
926
927/* Post IFFT complex multiply plus IFFT complex conjugate*/
928 for(i = 0; i < AC3_N / 4; i++)
929 {
930 tmp_a_r = buf[i].real;
931 tmp_a_i = -buf[i].imag;
932 buf[i].real = (tmp_a_r * mpeg3_xcos1[i]) - (tmp_a_i * mpeg3_xsin1[i]);
933 buf[i].imag = (tmp_a_r * mpeg3_xsin1[i]) + (tmp_a_i * mpeg3_xcos1[i]);
934 }
935
936 y_ptr = y;
937 delay_ptr = delay;
938 window_ptr = mpeg3_window;
939
940/* Window and convert to real valued signal */
941 for(i = 0; i < AC3_N / 8; i++)
942 {
943 *y_ptr = -buf[AC3_N / 8 + i].imag * *window_ptr++ + *delay_ptr++;
944 y_ptr += step;
945 *y_ptr = buf[AC3_N / 8 - i - 1].real * *window_ptr++ + *delay_ptr++;
946 y_ptr += step;
947 }
948
949 for(i = 0; i < AC3_N / 8; i++)
950 {
951 *y_ptr = -buf[i].real * *window_ptr++ + *delay_ptr++;
952 y_ptr += step;
953 *y_ptr = buf[AC3_N / 4 - i - 1].imag * *window_ptr++ + *delay_ptr++;
954 y_ptr += step;
955 }
956
957/* The trailing edge of the window goes into the delay line */
958 delay_ptr = delay;
959
960 for(i = 0; i < AC3_N / 8; i++)
961 {
962 *delay_ptr++ = -buf[AC3_N / 8 + i].real * *--window_ptr;
963 *delay_ptr++ = buf[AC3_N / 8 - i - 1].imag * *--window_ptr;
964 }
965
966 for(i = 0; i < AC3_N / 8; i++)
967 {
968 *delay_ptr++ = buf[i].imag * *--window_ptr;
969 *delay_ptr++ = -buf[AC3_N / 4 - i - 1].real * *--window_ptr;
970 }
971}
972
973void mpeg3audio_ac3_imdct_do_256(mpeg3audio_t *audio,
974 mpeg3_real_t data[],
975 mpeg3_real_t *y,
976 int step,
977 mpeg3_real_t *delay)
978{
979 int i, k;
980 int p, q;
981 int m;
982 int two_m;
983 int two_m_plus_one;
984 mpeg3_complex_t *buf = audio->ac3_imdct_buf;
985 mpeg3_real_t *y_ptr;
986 mpeg3_real_t *delay_ptr;
987 mpeg3_real_t *window_ptr;
988
989 mpeg3_real_t tmp_a_i;
990 mpeg3_real_t tmp_a_r;
991 mpeg3_real_t tmp_b_i;
992 mpeg3_real_t tmp_b_r;
993
994 mpeg3_complex_t *buf_1, *buf_2;
995
996 buf_1 = &buf[0];
997 buf_2 = &buf[64];
998
999/* Pre IFFT complex multiply plus IFFT cmplx conjugate */
1000 for(k = 0; k < AC3_N / 8; k++)
1001 {
1002 p = 2 * (AC3_N / 4 - 2 * k - 1);
1003 q = 2 * (2 * k);
1004
1005 buf_1[k].real = data[p] * mpeg3_xcos2[k] - data[q] * mpeg3_xsin2[k];
1006 buf_1[k].imag = - (data[q] * mpeg3_xcos2[k] + data[p] * mpeg3_xsin2[k]);
1007 buf_2[k].real = data[p + 1] * mpeg3_xcos2[k] - data[q + 1] * mpeg3_xsin2[k];
1008 buf_2[k].imag = - (data[q + 1] * mpeg3_xcos2[k] + data[p + 1] * mpeg3_xsin2[k]);
1009 }
1010
1011/* IFFT Bit reversed shuffling */
1012 for(i = 0; i < AC3_N / 8; i++)
1013 {
1014 k = mpeg3_bit_reverse_256[i];
1015 if(k < i)
1016 {
1017 swap_cmplx(&buf_1[i], &buf_1[k]);
1018 swap_cmplx(&buf_2[i], &buf_2[k]);
1019 }
1020 }
1021
1022/* FFT Merge */
1023 for(m = 0; m < 6; m++)
1024 {
1025 if(m)
1026 two_m = (1 << m);
1027 else
1028 two_m = 1;
1029
1030 two_m_plus_one = (1 << (m + 1));
1031
1032 for(k = 0; k < two_m; k++)
1033 {
1034 for(i = 0; i < AC3_N / 8; i += two_m_plus_one)
1035 {
1036 p = k + i;
1037 q = p + two_m;
1038/* Do block 1 */
1039 tmp_a_r = buf_1[p].real;
1040 tmp_a_i = buf_1[p].imag;
1041 tmp_b_r = buf_1[q].real * audio->ac3_w[m][k].real - buf_1[q].imag * audio->ac3_w[m][k].imag;
1042 tmp_b_i = buf_1[q].imag * audio->ac3_w[m][k].real + buf_1[q].real * audio->ac3_w[m][k].imag;
1043 buf_1[p].real = tmp_a_r + tmp_b_r;
1044 buf_1[p].imag = tmp_a_i + tmp_b_i;
1045 buf_1[q].real = tmp_a_r - tmp_b_r;
1046 buf_1[q].imag = tmp_a_i - tmp_b_i;
1047
1048/* Do block 2 */
1049 tmp_a_r = buf_2[p].real;
1050 tmp_a_i = buf_2[p].imag;
1051 tmp_b_r = buf_2[q].real * audio->ac3_w[m][k].real - buf_2[q].imag * audio->ac3_w[m][k].imag;
1052 tmp_b_i = buf_2[q].imag * audio->ac3_w[m][k].real + buf_2[q].real * audio->ac3_w[m][k].imag;
1053 buf_2[p].real = tmp_a_r + tmp_b_r;
1054 buf_2[p].imag = tmp_a_i + tmp_b_i;
1055 buf_2[q].real = tmp_a_r - tmp_b_r;
1056 buf_2[q].imag = tmp_a_i - tmp_b_i;
1057 }
1058 }
1059 }
1060
1061/* Post IFFT complex multiply */
1062 for(i = 0; i < AC3_N / 8; i++)
1063 {
1064 tmp_a_r = buf_1[i].real;
1065 tmp_a_i = -buf_1[i].imag;
1066 buf_1[i].real = (tmp_a_r * mpeg3_xcos2[i]) - (tmp_a_i * mpeg3_xsin2[i]);
1067 buf_1[i].imag = (tmp_a_r * mpeg3_xsin2[i]) + (tmp_a_i * mpeg3_xcos2[i]);
1068 tmp_a_r = buf_2[i].real;
1069 tmp_a_i = -buf_2[i].imag;
1070 buf_2[i].real = (tmp_a_r * mpeg3_xcos2[i]) - (tmp_a_i * mpeg3_xsin2[i]);
1071 buf_2[i].imag = (tmp_a_r * mpeg3_xsin2[i]) + (tmp_a_i * mpeg3_xcos2[i]);
1072 }
1073
1074/* Window and convert to real valued signal */
1075 y_ptr = y;
1076 delay_ptr = delay;
1077 window_ptr = mpeg3_window;
1078
1079 for(i = 0; i < AC3_N / 8; i++)
1080 {
1081 *y_ptr = -buf[AC3_N / 8 + i].imag * *window_ptr++ + *delay_ptr++;
1082 y_ptr += step;
1083 *y_ptr = buf[AC3_N / 8 - i - 1].real * *window_ptr++ + *delay_ptr++;
1084 y_ptr += step;
1085 }
1086
1087 for(i = 0; i < AC3_N / 8; i++)
1088 {
1089 *y_ptr = -buf[i].real * *window_ptr++ + *delay_ptr++;
1090 y_ptr += step;
1091 *y_ptr = buf[AC3_N / 4 - i - 1].imag * *window_ptr++ + *delay_ptr++;
1092 y_ptr += step;
1093 }
1094
1095/* The trailing edge of the window goes into the delay line */
1096 delay_ptr = delay;
1097
1098 for(i = 0; i < AC3_N / 8; i++)
1099 {
1100 *delay_ptr++ = -buf[AC3_N / 8 + i].real * *--window_ptr;
1101 *delay_ptr++ = buf[AC3_N / 8 - i - 1].imag * *--window_ptr;
1102 }
1103
1104 for(i = 0; i < AC3_N / 8; i++)
1105 {
1106 *delay_ptr++ = buf[i].imag * *--window_ptr;
1107 *delay_ptr++ = -buf[AC3_N / 4 - i - 1].real * *--window_ptr;
1108 }
1109}
1110
1111int mpeg3audio_ac3_imdct(mpeg3audio_t *audio,
1112 mpeg3_ac3bsi_t *bsi,
1113 mpeg3_ac3audblk_t *audblk,
1114 mpeg3ac3_stream_samples_t samples)
1115{
1116 int i;
1117
1118 for(i = 0; i < bsi->nfchans; i++)
1119 {
1120 if(audblk->blksw[i])
1121 mpeg3audio_ac3_imdct_do_256(audio,
1122 samples[i],
1123 audio->pcm_sample + audio->pcm_point + i,
1124 bsi->nfchans,
1125 audio->ac3_delay[i]);
1126 else
1127 mpeg3audio_ac3_imdct_do_512(audio,
1128 samples[i],
1129 audio->pcm_sample + audio->pcm_point + i,
1130 bsi->nfchans,
1131 audio->ac3_delay[i]);
1132 }
1133 audio->pcm_point += AC3_N / 2 * bsi->nfchans;
1134 return 0;
1135}