1 /********************************************************************
3 * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
5 * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
6 * PLEASE READ THESE TERMS DISTRIBUTING. *
8 * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000 *
9 * by Monty <monty@xiph.org> and The XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: normalized modified discrete cosine transform
15 power of two length transform only [16 <= n ]
16 last mod: $Id: mdct.c,v 1.17 2000/10/12 03:12:53 xiphmont Exp $
18 Algorithm adapted from _The use of multirate filter banks for coding
19 of high quality digital audio_, by T. Sporer, K. Brandenburg and
20 B. Edler, collection of the European Signal Processing Conference
21 (EUSIPCO), Amsterdam, June 1992, Vol.1, pp 211-214
23 Note that the below code won't make much sense without the paper;
24 The presented algorithm was already fairly polished, and the code
25 once followed it closely. The current code both corrects several
26 typos in the paper and goes beyond the presented optimizations
27 (steps 4 through 6 are, for example, entirely eliminated).
29 This module DOES NOT INCLUDE code to generate the window function.
30 Everybody has their own weird favorite including me... I happen to
31 like the properties of y=sin(2PI*sin^2(x)), but others may vehemently
34 ********************************************************************/
44 /* build lookups for trig functions; also pre-figure scaling and
45 some window function algebra. */
47 void mdct_init(mdct_lookup *lookup,int n){
48 int *bitrev=malloc(sizeof(int)*(n/4));
49 float *trig=malloc(sizeof(float)*(n+n/4));
58 int log2n=lookup->log2n=rint(log(n)/log(2));
61 lookup->bitrev=bitrev;
66 AE[i*2]=cos((M_PI/n)*(4*i));
67 AO[i*2]=-sin((M_PI/n)*(4*i));
68 BE[i*2]=cos((M_PI/(2*n))*(2*i+1));
69 BO[i*2]=sin((M_PI/(2*n))*(2*i+1));
72 CE[i*2]=cos((M_PI/n)*(4*i+2));
73 CO[i*2]=-sin((M_PI/n)*(4*i+2));
76 /* bitreverse lookup... */
79 int mask=(1<<(log2n-1))-1,i,j;
84 if((msb>>j)&i)acc|=1<<j;
85 bitrev[i*2]=((~acc)&mask);
91 void mdct_clear(mdct_lookup *l){
93 if(l->trig)free(l->trig);
94 if(l->bitrev)free(l->bitrev);
95 memset(l,0,sizeof(mdct_lookup));
99 static float *_mdct_kernel(float *x, float *w,
100 int n, int n2, int n4, int n8,
109 float *A=init->trig+n2;
114 w2[i]= *xA++ + *xB++;
120 w[i++]= x0 * A[0] + x1 * A[1];
121 w[i]= x1 * A[0] - x0 * A[1];
123 w2[i++]= *xA++ + *xB++;
132 for(i=0;i<init->log2n-3;i++){
139 for(r=0;r<(k0>>2);r++){
147 for(s=0;s<(2<<i);s++){
151 wA =w[++w1] -w[++w2];
154 x[w2] =wA*AEv - wB*AOv;
155 x[w2-1]=wB*AEv + wA*AOv;
171 /* step 4, 5, 6, 7 */
173 float *C=init->trig+n;
174 int *bit=init->bitrev;
181 float wA=w[t1]-w[t2+1];
182 float wB=w[t1-1]+w[t2];
183 float wC=w[t1]+w[t2+1];
184 float wD=w[t1-1]-w[t2];
191 *x1++=( wC+wACO+wBCE)*.5;
192 *x2--=(-wD+wBCO-wACE)*.5;
193 *x1++=( wD+wBCO-wACE)*.5;
194 *x2--=( wC-wACO-wBCE)*.5;
200 void mdct_forward(mdct_lookup *init, float *in, float *out){
202 float *x=alloca(sizeof(float)*(n/2));
203 float *w=alloca(sizeof(float)*(n/2));
210 /* window + rotate + step 1 */
215 float *A=init->trig+n2;
221 tempA= in[in1+2] + in[in2];
222 tempB= in[in1] + in[in2+2];
224 x[i]= tempB*A[1] + tempA*A[0];
225 x[i+1]= tempB*A[0] - tempA*A[1];
232 tempA= in[in1+2] - in[in2];
233 tempB= in[in1] - in[in2+2];
235 x[i]= tempB*A[1] + tempA*A[0];
236 x[i+1]= tempB*A[0] - tempA*A[1];
243 tempA= -in[in1+2] - in[in2];
244 tempB= -in[in1] - in[in2+2];
246 x[i]= tempB*A[1] + tempA*A[0];
247 x[i+1]= tempB*A[0] - tempA*A[1];
251 xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
256 float *B=init->trig+n2;
260 out[i] =(xx[0]*B[0]+xx[1]*B[1])*scale;
261 *(--out2)=(xx[0]*B[1]-xx[1]*B[0])*scale;
269 void mdct_backward(mdct_lookup *init, float *in, float *out){
271 float *x=alloca(sizeof(float)*(n/2));
272 float *w=alloca(sizeof(float)*(n/2));
279 /* rotate + step 1 */
283 float *A=init->trig+n2;
287 *xO++=-*(inO+2)*A[1] - *inO*A[0];
288 *xO++= *inO*A[1] - *(inO+2)*A[0];
296 *xO++=*inO*A[1] + *(inO+2)*A[0];
297 *xO++=*inO*A[0] - *(inO+2)*A[1];
303 xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
308 float *B=init->trig+n2;
310 int o3=n4+n2,o4=o3-1;
313 float temp1= (*xx * B[1] - *(xx+1) * B[0]);
314 float temp2=-(*xx * B[0] + *(xx+1) * B[1]);