1 | /******************************************************************** |
---|
2 | * * |
---|
3 | * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * |
---|
4 | * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * |
---|
5 | * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * |
---|
6 | * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * |
---|
7 | * * |
---|
8 | * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2007 * |
---|
9 | * by the Xiph.Org Foundation http://www.xiph.org/ * |
---|
10 | * * |
---|
11 | ******************************************************************** |
---|
12 | |
---|
13 | function: normalized modified discrete cosine transform |
---|
14 | power of two length transform only [64 <= n ] |
---|
15 | last mod: $Id: mdct.c 13293 2007-07-24 00:09:47Z xiphmont $ |
---|
16 | |
---|
17 | Original algorithm adapted long ago from _The use of multirate filter |
---|
18 | banks for coding of high quality digital audio_, by T. Sporer, |
---|
19 | K. Brandenburg and B. Edler, collection of the European Signal |
---|
20 | Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp |
---|
21 | 211-214 |
---|
22 | |
---|
23 | The below code implements an algorithm that no longer looks much like |
---|
24 | that presented in the paper, but the basic structure remains if you |
---|
25 | dig deep enough to see it. |
---|
26 | |
---|
27 | This module DOES NOT INCLUDE code to generate/apply the window |
---|
28 | function. Everybody has their own weird favorite including me... I |
---|
29 | happen to like the properties of y=sin(.5PI*sin^2(x)), but others may |
---|
30 | vehemently disagree. |
---|
31 | |
---|
32 | ********************************************************************/ |
---|
33 | |
---|
34 | /* this can also be run as an integer transform by uncommenting a |
---|
35 | define in mdct.h; the integerization is a first pass and although |
---|
36 | it's likely stable for Vorbis, the dynamic range is constrained and |
---|
37 | roundoff isn't done (so it's noisy). Consider it functional, but |
---|
38 | only a starting point. There's no point on a machine with an FPU */ |
---|
39 | |
---|
40 | #include <stdio.h> |
---|
41 | #include <stdlib.h> |
---|
42 | #include <string.h> |
---|
43 | #include <math.h> |
---|
44 | #include "vorbis/codec.h" |
---|
45 | #include "mdct.h" |
---|
46 | #include "os.h" |
---|
47 | #include "misc.h" |
---|
48 | |
---|
49 | /* build lookups for trig functions; also pre-figure scaling and |
---|
50 | some window function algebra. */ |
---|
51 | |
---|
52 | void mdct_init(mdct_lookup *lookup,int n){ |
---|
53 | int *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4)); |
---|
54 | DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4)); |
---|
55 | |
---|
56 | int i; |
---|
57 | int n2=n>>1; |
---|
58 | int log2n=lookup->log2n=rint(log((float)n)/log(2.f)); |
---|
59 | lookup->n=n; |
---|
60 | lookup->trig=T; |
---|
61 | lookup->bitrev=bitrev; |
---|
62 | |
---|
63 | /* trig lookups... */ |
---|
64 | |
---|
65 | for(i=0;i<n/4;i++){ |
---|
66 | T[i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i))); |
---|
67 | T[i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i))); |
---|
68 | T[n2+i*2]=FLOAT_CONV(cos((M_PI/(2*n))*(2*i+1))); |
---|
69 | T[n2+i*2+1]=FLOAT_CONV(sin((M_PI/(2*n))*(2*i+1))); |
---|
70 | } |
---|
71 | for(i=0;i<n/8;i++){ |
---|
72 | T[n+i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i+2))*.5); |
---|
73 | T[n+i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i+2))*.5); |
---|
74 | } |
---|
75 | |
---|
76 | /* bitreverse lookup... */ |
---|
77 | |
---|
78 | { |
---|
79 | int mask=(1<<(log2n-1))-1,i,j; |
---|
80 | int msb=1<<(log2n-2); |
---|
81 | for(i=0;i<n/8;i++){ |
---|
82 | int acc=0; |
---|
83 | for(j=0;msb>>j;j++) |
---|
84 | if((msb>>j)&i)acc|=1<<j; |
---|
85 | bitrev[i*2]=((~acc)&mask)-1; |
---|
86 | bitrev[i*2+1]=acc; |
---|
87 | |
---|
88 | } |
---|
89 | } |
---|
90 | lookup->scale=FLOAT_CONV(4.f/n); |
---|
91 | } |
---|
92 | |
---|
93 | /* 8 point butterfly (in place, 4 register) */ |
---|
94 | STIN void mdct_butterfly_8(DATA_TYPE *x){ |
---|
95 | REG_TYPE r0 = x[6] + x[2]; |
---|
96 | REG_TYPE r1 = x[6] - x[2]; |
---|
97 | REG_TYPE r2 = x[4] + x[0]; |
---|
98 | REG_TYPE r3 = x[4] - x[0]; |
---|
99 | |
---|
100 | x[6] = r0 + r2; |
---|
101 | x[4] = r0 - r2; |
---|
102 | |
---|
103 | r0 = x[5] - x[1]; |
---|
104 | r2 = x[7] - x[3]; |
---|
105 | x[0] = r1 + r0; |
---|
106 | x[2] = r1 - r0; |
---|
107 | |
---|
108 | r0 = x[5] + x[1]; |
---|
109 | r1 = x[7] + x[3]; |
---|
110 | x[3] = r2 + r3; |
---|
111 | x[1] = r2 - r3; |
---|
112 | x[7] = r1 + r0; |
---|
113 | x[5] = r1 - r0; |
---|
114 | |
---|
115 | } |
---|
116 | |
---|
117 | /* 16 point butterfly (in place, 4 register) */ |
---|
118 | STIN void mdct_butterfly_16(DATA_TYPE *x){ |
---|
119 | REG_TYPE r0 = x[1] - x[9]; |
---|
120 | REG_TYPE r1 = x[0] - x[8]; |
---|
121 | |
---|
122 | x[8] += x[0]; |
---|
123 | x[9] += x[1]; |
---|
124 | x[0] = MULT_NORM((r0 + r1) * cPI2_8); |
---|
125 | x[1] = MULT_NORM((r0 - r1) * cPI2_8); |
---|
126 | |
---|
127 | r0 = x[3] - x[11]; |
---|
128 | r1 = x[10] - x[2]; |
---|
129 | x[10] += x[2]; |
---|
130 | x[11] += x[3]; |
---|
131 | x[2] = r0; |
---|
132 | x[3] = r1; |
---|
133 | |
---|
134 | r0 = x[12] - x[4]; |
---|
135 | r1 = x[13] - x[5]; |
---|
136 | x[12] += x[4]; |
---|
137 | x[13] += x[5]; |
---|
138 | x[4] = MULT_NORM((r0 - r1) * cPI2_8); |
---|
139 | x[5] = MULT_NORM((r0 + r1) * cPI2_8); |
---|
140 | |
---|
141 | r0 = x[14] - x[6]; |
---|
142 | r1 = x[15] - x[7]; |
---|
143 | x[14] += x[6]; |
---|
144 | x[15] += x[7]; |
---|
145 | x[6] = r0; |
---|
146 | x[7] = r1; |
---|
147 | |
---|
148 | mdct_butterfly_8(x); |
---|
149 | mdct_butterfly_8(x+8); |
---|
150 | } |
---|
151 | |
---|
152 | /* 32 point butterfly (in place, 4 register) */ |
---|
153 | STIN void mdct_butterfly_32(DATA_TYPE *x){ |
---|
154 | REG_TYPE r0 = x[30] - x[14]; |
---|
155 | REG_TYPE r1 = x[31] - x[15]; |
---|
156 | |
---|
157 | x[30] += x[14]; |
---|
158 | x[31] += x[15]; |
---|
159 | x[14] = r0; |
---|
160 | x[15] = r1; |
---|
161 | |
---|
162 | r0 = x[28] - x[12]; |
---|
163 | r1 = x[29] - x[13]; |
---|
164 | x[28] += x[12]; |
---|
165 | x[29] += x[13]; |
---|
166 | x[12] = MULT_NORM( r0 * cPI1_8 - r1 * cPI3_8 ); |
---|
167 | x[13] = MULT_NORM( r0 * cPI3_8 + r1 * cPI1_8 ); |
---|
168 | |
---|
169 | r0 = x[26] - x[10]; |
---|
170 | r1 = x[27] - x[11]; |
---|
171 | x[26] += x[10]; |
---|
172 | x[27] += x[11]; |
---|
173 | x[10] = MULT_NORM(( r0 - r1 ) * cPI2_8); |
---|
174 | x[11] = MULT_NORM(( r0 + r1 ) * cPI2_8); |
---|
175 | |
---|
176 | r0 = x[24] - x[8]; |
---|
177 | r1 = x[25] - x[9]; |
---|
178 | x[24] += x[8]; |
---|
179 | x[25] += x[9]; |
---|
180 | x[8] = MULT_NORM( r0 * cPI3_8 - r1 * cPI1_8 ); |
---|
181 | x[9] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 ); |
---|
182 | |
---|
183 | r0 = x[22] - x[6]; |
---|
184 | r1 = x[7] - x[23]; |
---|
185 | x[22] += x[6]; |
---|
186 | x[23] += x[7]; |
---|
187 | x[6] = r1; |
---|
188 | x[7] = r0; |
---|
189 | |
---|
190 | r0 = x[4] - x[20]; |
---|
191 | r1 = x[5] - x[21]; |
---|
192 | x[20] += x[4]; |
---|
193 | x[21] += x[5]; |
---|
194 | x[4] = MULT_NORM( r1 * cPI1_8 + r0 * cPI3_8 ); |
---|
195 | x[5] = MULT_NORM( r1 * cPI3_8 - r0 * cPI1_8 ); |
---|
196 | |
---|
197 | r0 = x[2] - x[18]; |
---|
198 | r1 = x[3] - x[19]; |
---|
199 | x[18] += x[2]; |
---|
200 | x[19] += x[3]; |
---|
201 | x[2] = MULT_NORM(( r1 + r0 ) * cPI2_8); |
---|
202 | x[3] = MULT_NORM(( r1 - r0 ) * cPI2_8); |
---|
203 | |
---|
204 | r0 = x[0] - x[16]; |
---|
205 | r1 = x[1] - x[17]; |
---|
206 | x[16] += x[0]; |
---|
207 | x[17] += x[1]; |
---|
208 | x[0] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 ); |
---|
209 | x[1] = MULT_NORM( r1 * cPI1_8 - r0 * cPI3_8 ); |
---|
210 | |
---|
211 | mdct_butterfly_16(x); |
---|
212 | mdct_butterfly_16(x+16); |
---|
213 | |
---|
214 | } |
---|
215 | |
---|
216 | /* N point first stage butterfly (in place, 2 register) */ |
---|
217 | STIN void mdct_butterfly_first(DATA_TYPE *T, |
---|
218 | DATA_TYPE *x, |
---|
219 | int points){ |
---|
220 | |
---|
221 | DATA_TYPE *x1 = x + points - 8; |
---|
222 | DATA_TYPE *x2 = x + (points>>1) - 8; |
---|
223 | REG_TYPE r0; |
---|
224 | REG_TYPE r1; |
---|
225 | |
---|
226 | do{ |
---|
227 | |
---|
228 | r0 = x1[6] - x2[6]; |
---|
229 | r1 = x1[7] - x2[7]; |
---|
230 | x1[6] += x2[6]; |
---|
231 | x1[7] += x2[7]; |
---|
232 | x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]); |
---|
233 | x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]); |
---|
234 | |
---|
235 | r0 = x1[4] - x2[4]; |
---|
236 | r1 = x1[5] - x2[5]; |
---|
237 | x1[4] += x2[4]; |
---|
238 | x1[5] += x2[5]; |
---|
239 | x2[4] = MULT_NORM(r1 * T[5] + r0 * T[4]); |
---|
240 | x2[5] = MULT_NORM(r1 * T[4] - r0 * T[5]); |
---|
241 | |
---|
242 | r0 = x1[2] - x2[2]; |
---|
243 | r1 = x1[3] - x2[3]; |
---|
244 | x1[2] += x2[2]; |
---|
245 | x1[3] += x2[3]; |
---|
246 | x2[2] = MULT_NORM(r1 * T[9] + r0 * T[8]); |
---|
247 | x2[3] = MULT_NORM(r1 * T[8] - r0 * T[9]); |
---|
248 | |
---|
249 | r0 = x1[0] - x2[0]; |
---|
250 | r1 = x1[1] - x2[1]; |
---|
251 | x1[0] += x2[0]; |
---|
252 | x1[1] += x2[1]; |
---|
253 | x2[0] = MULT_NORM(r1 * T[13] + r0 * T[12]); |
---|
254 | x2[1] = MULT_NORM(r1 * T[12] - r0 * T[13]); |
---|
255 | |
---|
256 | x1-=8; |
---|
257 | x2-=8; |
---|
258 | T+=16; |
---|
259 | |
---|
260 | }while(x2>=x); |
---|
261 | } |
---|
262 | |
---|
263 | /* N/stage point generic N stage butterfly (in place, 2 register) */ |
---|
264 | STIN void mdct_butterfly_generic(DATA_TYPE *T, |
---|
265 | DATA_TYPE *x, |
---|
266 | int points, |
---|
267 | int trigint){ |
---|
268 | |
---|
269 | DATA_TYPE *x1 = x + points - 8; |
---|
270 | DATA_TYPE *x2 = x + (points>>1) - 8; |
---|
271 | REG_TYPE r0; |
---|
272 | REG_TYPE r1; |
---|
273 | |
---|
274 | do{ |
---|
275 | |
---|
276 | r0 = x1[6] - x2[6]; |
---|
277 | r1 = x1[7] - x2[7]; |
---|
278 | x1[6] += x2[6]; |
---|
279 | x1[7] += x2[7]; |
---|
280 | x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]); |
---|
281 | x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]); |
---|
282 | |
---|
283 | T+=trigint; |
---|
284 | |
---|
285 | r0 = x1[4] - x2[4]; |
---|
286 | r1 = x1[5] - x2[5]; |
---|
287 | x1[4] += x2[4]; |
---|
288 | x1[5] += x2[5]; |
---|
289 | x2[4] = MULT_NORM(r1 * T[1] + r0 * T[0]); |
---|
290 | x2[5] = MULT_NORM(r1 * T[0] - r0 * T[1]); |
---|
291 | |
---|
292 | T+=trigint; |
---|
293 | |
---|
294 | r0 = x1[2] - x2[2]; |
---|
295 | r1 = x1[3] - x2[3]; |
---|
296 | x1[2] += x2[2]; |
---|
297 | x1[3] += x2[3]; |
---|
298 | x2[2] = MULT_NORM(r1 * T[1] + r0 * T[0]); |
---|
299 | x2[3] = MULT_NORM(r1 * T[0] - r0 * T[1]); |
---|
300 | |
---|
301 | T+=trigint; |
---|
302 | |
---|
303 | r0 = x1[0] - x2[0]; |
---|
304 | r1 = x1[1] - x2[1]; |
---|
305 | x1[0] += x2[0]; |
---|
306 | x1[1] += x2[1]; |
---|
307 | x2[0] = MULT_NORM(r1 * T[1] + r0 * T[0]); |
---|
308 | x2[1] = MULT_NORM(r1 * T[0] - r0 * T[1]); |
---|
309 | |
---|
310 | T+=trigint; |
---|
311 | x1-=8; |
---|
312 | x2-=8; |
---|
313 | |
---|
314 | }while(x2>=x); |
---|
315 | } |
---|
316 | |
---|
317 | STIN void mdct_butterflies(mdct_lookup *init, |
---|
318 | DATA_TYPE *x, |
---|
319 | int points){ |
---|
320 | |
---|
321 | DATA_TYPE *T=init->trig; |
---|
322 | int stages=init->log2n-5; |
---|
323 | int i,j; |
---|
324 | |
---|
325 | if(--stages>0){ |
---|
326 | mdct_butterfly_first(T,x,points); |
---|
327 | } |
---|
328 | |
---|
329 | for(i=1;--stages>0;i++){ |
---|
330 | for(j=0;j<(1<<i);j++) |
---|
331 | mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i); |
---|
332 | } |
---|
333 | |
---|
334 | for(j=0;j<points;j+=32) |
---|
335 | mdct_butterfly_32(x+j); |
---|
336 | |
---|
337 | } |
---|
338 | |
---|
339 | void mdct_clear(mdct_lookup *l){ |
---|
340 | if(l){ |
---|
341 | if(l->trig)_ogg_free(l->trig); |
---|
342 | if(l->bitrev)_ogg_free(l->bitrev); |
---|
343 | memset(l,0,sizeof(*l)); |
---|
344 | } |
---|
345 | } |
---|
346 | |
---|
347 | STIN void mdct_bitreverse(mdct_lookup *init, |
---|
348 | DATA_TYPE *x){ |
---|
349 | int n = init->n; |
---|
350 | int *bit = init->bitrev; |
---|
351 | DATA_TYPE *w0 = x; |
---|
352 | DATA_TYPE *w1 = x = w0+(n>>1); |
---|
353 | DATA_TYPE *T = init->trig+n; |
---|
354 | |
---|
355 | do{ |
---|
356 | DATA_TYPE *x0 = x+bit[0]; |
---|
357 | DATA_TYPE *x1 = x+bit[1]; |
---|
358 | |
---|
359 | REG_TYPE r0 = x0[1] - x1[1]; |
---|
360 | REG_TYPE r1 = x0[0] + x1[0]; |
---|
361 | REG_TYPE r2 = MULT_NORM(r1 * T[0] + r0 * T[1]); |
---|
362 | REG_TYPE r3 = MULT_NORM(r1 * T[1] - r0 * T[0]); |
---|
363 | |
---|
364 | w1 -= 4; |
---|
365 | |
---|
366 | r0 = HALVE(x0[1] + x1[1]); |
---|
367 | r1 = HALVE(x0[0] - x1[0]); |
---|
368 | |
---|
369 | w0[0] = r0 + r2; |
---|
370 | w1[2] = r0 - r2; |
---|
371 | w0[1] = r1 + r3; |
---|
372 | w1[3] = r3 - r1; |
---|
373 | |
---|
374 | x0 = x+bit[2]; |
---|
375 | x1 = x+bit[3]; |
---|
376 | |
---|
377 | r0 = x0[1] - x1[1]; |
---|
378 | r1 = x0[0] + x1[0]; |
---|
379 | r2 = MULT_NORM(r1 * T[2] + r0 * T[3]); |
---|
380 | r3 = MULT_NORM(r1 * T[3] - r0 * T[2]); |
---|
381 | |
---|
382 | r0 = HALVE(x0[1] + x1[1]); |
---|
383 | r1 = HALVE(x0[0] - x1[0]); |
---|
384 | |
---|
385 | w0[2] = r0 + r2; |
---|
386 | w1[0] = r0 - r2; |
---|
387 | w0[3] = r1 + r3; |
---|
388 | w1[1] = r3 - r1; |
---|
389 | |
---|
390 | T += 4; |
---|
391 | bit += 4; |
---|
392 | w0 += 4; |
---|
393 | |
---|
394 | }while(w0<w1); |
---|
395 | } |
---|
396 | |
---|
397 | void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){ |
---|
398 | int n=init->n; |
---|
399 | int n2=n>>1; |
---|
400 | int n4=n>>2; |
---|
401 | |
---|
402 | /* rotate */ |
---|
403 | |
---|
404 | DATA_TYPE *iX = in+n2-7; |
---|
405 | DATA_TYPE *oX = out+n2+n4; |
---|
406 | DATA_TYPE *T = init->trig+n4; |
---|
407 | |
---|
408 | do{ |
---|
409 | oX -= 4; |
---|
410 | oX[0] = MULT_NORM(-iX[2] * T[3] - iX[0] * T[2]); |
---|
411 | oX[1] = MULT_NORM (iX[0] * T[3] - iX[2] * T[2]); |
---|
412 | oX[2] = MULT_NORM(-iX[6] * T[1] - iX[4] * T[0]); |
---|
413 | oX[3] = MULT_NORM (iX[4] * T[1] - iX[6] * T[0]); |
---|
414 | iX -= 8; |
---|
415 | T += 4; |
---|
416 | }while(iX>=in); |
---|
417 | |
---|
418 | iX = in+n2-8; |
---|
419 | oX = out+n2+n4; |
---|
420 | T = init->trig+n4; |
---|
421 | |
---|
422 | do{ |
---|
423 | T -= 4; |
---|
424 | oX[0] = MULT_NORM (iX[4] * T[3] + iX[6] * T[2]); |
---|
425 | oX[1] = MULT_NORM (iX[4] * T[2] - iX[6] * T[3]); |
---|
426 | oX[2] = MULT_NORM (iX[0] * T[1] + iX[2] * T[0]); |
---|
427 | oX[3] = MULT_NORM (iX[0] * T[0] - iX[2] * T[1]); |
---|
428 | iX -= 8; |
---|
429 | oX += 4; |
---|
430 | }while(iX>=in); |
---|
431 | |
---|
432 | mdct_butterflies(init,out+n2,n2); |
---|
433 | mdct_bitreverse(init,out); |
---|
434 | |
---|
435 | /* roatate + window */ |
---|
436 | |
---|
437 | { |
---|
438 | DATA_TYPE *oX1=out+n2+n4; |
---|
439 | DATA_TYPE *oX2=out+n2+n4; |
---|
440 | DATA_TYPE *iX =out; |
---|
441 | T =init->trig+n2; |
---|
442 | |
---|
443 | do{ |
---|
444 | oX1-=4; |
---|
445 | |
---|
446 | oX1[3] = MULT_NORM (iX[0] * T[1] - iX[1] * T[0]); |
---|
447 | oX2[0] = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]); |
---|
448 | |
---|
449 | oX1[2] = MULT_NORM (iX[2] * T[3] - iX[3] * T[2]); |
---|
450 | oX2[1] = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]); |
---|
451 | |
---|
452 | oX1[1] = MULT_NORM (iX[4] * T[5] - iX[5] * T[4]); |
---|
453 | oX2[2] = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]); |
---|
454 | |
---|
455 | oX1[0] = MULT_NORM (iX[6] * T[7] - iX[7] * T[6]); |
---|
456 | oX2[3] = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]); |
---|
457 | |
---|
458 | oX2+=4; |
---|
459 | iX += 8; |
---|
460 | T += 8; |
---|
461 | }while(iX<oX1); |
---|
462 | |
---|
463 | iX=out+n2+n4; |
---|
464 | oX1=out+n4; |
---|
465 | oX2=oX1; |
---|
466 | |
---|
467 | do{ |
---|
468 | oX1-=4; |
---|
469 | iX-=4; |
---|
470 | |
---|
471 | oX2[0] = -(oX1[3] = iX[3]); |
---|
472 | oX2[1] = -(oX1[2] = iX[2]); |
---|
473 | oX2[2] = -(oX1[1] = iX[1]); |
---|
474 | oX2[3] = -(oX1[0] = iX[0]); |
---|
475 | |
---|
476 | oX2+=4; |
---|
477 | }while(oX2<iX); |
---|
478 | |
---|
479 | iX=out+n2+n4; |
---|
480 | oX1=out+n2+n4; |
---|
481 | oX2=out+n2; |
---|
482 | do{ |
---|
483 | oX1-=4; |
---|
484 | oX1[0]= iX[3]; |
---|
485 | oX1[1]= iX[2]; |
---|
486 | oX1[2]= iX[1]; |
---|
487 | oX1[3]= iX[0]; |
---|
488 | iX+=4; |
---|
489 | }while(oX1>oX2); |
---|
490 | } |
---|
491 | } |
---|
492 | |
---|
493 | void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){ |
---|
494 | int n=init->n; |
---|
495 | int n2=n>>1; |
---|
496 | int n4=n>>2; |
---|
497 | int n8=n>>3; |
---|
498 | DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */ |
---|
499 | DATA_TYPE *w2=w+n2; |
---|
500 | |
---|
501 | /* rotate */ |
---|
502 | |
---|
503 | /* window + rotate + step 1 */ |
---|
504 | |
---|
505 | REG_TYPE r0; |
---|
506 | REG_TYPE r1; |
---|
507 | DATA_TYPE *x0=in+n2+n4; |
---|
508 | DATA_TYPE *x1=x0+1; |
---|
509 | DATA_TYPE *T=init->trig+n2; |
---|
510 | |
---|
511 | int i=0; |
---|
512 | |
---|
513 | for(i=0;i<n8;i+=2){ |
---|
514 | x0 -=4; |
---|
515 | T-=2; |
---|
516 | r0= x0[2] + x1[0]; |
---|
517 | r1= x0[0] + x1[2]; |
---|
518 | w2[i]= MULT_NORM(r1*T[1] + r0*T[0]); |
---|
519 | w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]); |
---|
520 | x1 +=4; |
---|
521 | } |
---|
522 | |
---|
523 | x1=in+1; |
---|
524 | |
---|
525 | for(;i<n2-n8;i+=2){ |
---|
526 | T-=2; |
---|
527 | x0 -=4; |
---|
528 | r0= x0[2] - x1[0]; |
---|
529 | r1= x0[0] - x1[2]; |
---|
530 | w2[i]= MULT_NORM(r1*T[1] + r0*T[0]); |
---|
531 | w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]); |
---|
532 | x1 +=4; |
---|
533 | } |
---|
534 | |
---|
535 | x0=in+n; |
---|
536 | |
---|
537 | for(;i<n2;i+=2){ |
---|
538 | T-=2; |
---|
539 | x0 -=4; |
---|
540 | r0= -x0[2] - x1[0]; |
---|
541 | r1= -x0[0] - x1[2]; |
---|
542 | w2[i]= MULT_NORM(r1*T[1] + r0*T[0]); |
---|
543 | w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]); |
---|
544 | x1 +=4; |
---|
545 | } |
---|
546 | |
---|
547 | |
---|
548 | mdct_butterflies(init,w+n2,n2); |
---|
549 | mdct_bitreverse(init,w); |
---|
550 | |
---|
551 | /* roatate + window */ |
---|
552 | |
---|
553 | T=init->trig+n2; |
---|
554 | x0=out+n2; |
---|
555 | |
---|
556 | for(i=0;i<n4;i++){ |
---|
557 | x0--; |
---|
558 | out[i] =MULT_NORM((w[0]*T[0]+w[1]*T[1])*init->scale); |
---|
559 | x0[0] =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale); |
---|
560 | w+=2; |
---|
561 | T+=2; |
---|
562 | } |
---|
563 | } |
---|
564 | |
---|