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-2009 *
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 16227 2009-07-08 06:58:46Z 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
mdct_init(mdct_lookup * lookup,int n)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) */
mdct_butterfly_8(DATA_TYPE * x)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) */
mdct_butterfly_16(DATA_TYPE * x)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) */
mdct_butterfly_32(DATA_TYPE * x)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) */
mdct_butterfly_first(DATA_TYPE * T,DATA_TYPE * x,int points)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) */
mdct_butterfly_generic(DATA_TYPE * T,DATA_TYPE * x,int points,int trigint)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
mdct_butterflies(mdct_lookup * init,DATA_TYPE * x,int points)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
mdct_clear(mdct_lookup * l)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
mdct_bitreverse(mdct_lookup * init,DATA_TYPE * x)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
mdct_backward(mdct_lookup * init,DATA_TYPE * in,DATA_TYPE * out)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
mdct_forward(mdct_lookup * init,DATA_TYPE * in,DATA_TYPE * out)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