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-1999 *
9 * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: modified discrete cosine transform
15 power of two length transform only [16 <= n ]
17 author: Monty <xiphmont@mit.edu>
18 modifications by: Monty
19 last modification date: Jun 04 1999
21 Algorithm adapted from _The use of multirate filter banks for coding
22 of high quality digital audio_, by T. Sporer, K. Brandenburg and
23 B. Edler, collection of the European Signal Processing Conference
24 (EUSIPCO), Amsterdam, June 1992, Vol.1, pp 211-214
26 Note that the below code won't make much sense without the paper;
27 The presented algorithm was already fairly polished, and the code
28 once followed it closely. The current code both corrects several
29 typos in the paper and goes beyond the presented optimizations
30 (steps 4 through 6 are, for example, entirely eliminated).
32 This module DOES NOT INCLUDE code to generate the window function.
33 Everybody has their own weird favorite including me... I happen to
34 like the properties of y=sin(2PI*sin^2(x)), but others may vehemently
37 ********************************************************************/
39 static const char rcsid[] = "$Id";
46 static double oPI = 3.14159265358979323846;
48 /* build lookups for trig functions; also pre-figure scaling and
49 some window function algebra. */
51 MDCT_lookup *MDCT_init(int n){
52 MDCT_lookup *lookup=malloc(sizeof(MDCT_lookup));
53 int *bitrev=malloc(sizeof(int)*(n/4));
54 double *trig=malloc(sizeof(double)*(n+n/4));
66 int log2n=lookup->log2n=rint(log(n)/log(2));
69 lookup->bitrev=bitrev;
74 AE[i]=cos((oPI/n)*(4*i));
75 AO[i]=-sin((oPI/n)*(4*i));
76 BE[i]=cos((oPI/(2*n))*(2*i+1));
77 BO[i]=sin((oPI/(2*n))*(2*i+1));
80 CE[i]=cos((oPI/n)*(4*i+2));
81 CO[i]=-sin((oPI/n)*(4*i+2));
84 /* bitreverse lookup... */
87 int mask=(1<<(log2n-1))-1,i,j;
92 if((msb>>j)&i)acc|=1<<j;
93 bitA[i]=((~acc)&mask)*2;
101 void MDCT_free(MDCT_lookup *l){
103 if(l->trig)free(l->trig);
104 if(l->bitrev)free(l->bitrev);
109 static inline void _MDCT_kernel(double *x,
110 int n, int n2, int n4, int n8,
112 double *w=x+1; /* interleaved access improves cache locality */
120 double *AE=init->trig+n4;
124 double x0=xA[i]-xB[i];
125 double x1=xA[i+2]-xB[i+2];
128 w[i] =x0 * *AE + x1 * *AO;
131 w[i] =x1 * *AE - x0 * *AO;
141 for(i=0;i<init->log2n-3;i++){
145 double *AE=init->trig;
149 for(r=0;r<(n4>>i);r+=4){
151 int w2=wbase-(k0>>1);
154 for(s=0;s<(2<<i);s++){
155 x[w1+2]=w[w1+2]+w[w2+2];
157 x[w2+2]=(w[w1+2]-w[w2+2])* *AE-(w[w1]-w[w2])* *AO;
158 x[w2] =(w[w1]-w[w2])* *AE+(w[w1+2]-w[w2+2])* *AO;
173 /* step 4, 5, 6, 7 */
175 double *CE=init->trig+n;
177 int *bitA=init->bitrev;
187 double wA=w[t1]-w[t2];
188 double wB=w[t3]+w[t4];
189 double wC=w[t1]+w[t2];
190 double wD=w[t3]-w[t4];
193 double wBCO=wB* *(CO++);
195 double wBCE=wB* *(CE++);
197 *x1 =( wC+wACO+wBCE)*.5;
198 *(x2-2)=( wC-wACO-wBCE)*.5;
199 *(x1+2)=( wD+wBCO-wACE)*.5;
200 *x2 =(-wD+wBCO-wACE)*.5;
207 void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
209 double *x=alloca(n*sizeof(double));
215 /* window + rotate + step 1 */
220 double *AE=init->trig+n4;
226 tempA= in[in1+2]*window[in1+2] + in[in2]*window[in2];
227 tempB= in[in1]*window[in1] + in[in2+2]*window[in2+2];
230 x[i]= tempB* *AO + tempA* *AE;
231 x[i+2]= tempB* *AE - tempA* *AO;
237 tempA= in[in1+2]*window[in1+2] - in[in2]*window[in2];
238 tempB= in[in1]*window[in1] - in[in2+2]*window[in2+2];
241 x[i]= tempB* *AO + tempA* *AE;
242 x[i+2]= tempB* *AE - tempA* *AO;
248 tempA= -in[in1+2]*window[in1+2] - in[in2]*window[in2];
249 tempB= -in[in1]*window[in1] - in[in2+2]*window[in2+2];
252 x[i]= tempB* *AO + tempA* *AE;
253 x[i+2]= tempB* *AE - tempA* *AO;
257 _MDCT_kernel(x,n,n2,n4,n8,init);
262 double *BE=init->trig+n2;
266 out[i] =x[0]* *BE+x[2]* *BO;
267 *(--out2)=x[0]* *BO-x[2]* *BE;
275 void iMDCT(double *in, double *out, MDCT_lookup *init, double *window){
277 double *x=alloca(n*sizeof(double));
283 /* window + rotate + step 1 */
287 double *AE=init->trig+n4;
292 *xO=-*(inO+2)* *AO - *inO * *AE;
294 *xO= *inO * *AO - *(inO+2)* *AE;
303 *xO=*inO * *AO + *(inO+2) * *AE;
305 *xO=*inO * *AE - *(inO+2) * *AO;
312 _MDCT_kernel(x,n,n2,n4,n8,init);
317 double *BE=init->trig+n2;
320 int o3=n4+n2,o4=o3-1;
324 double temp1= (*x * *BO - *(x+2) * *BE)* scale;
325 double temp2= (*x * *BE + *(x+2) * *BO)* -scale;
327 out[o1]=-temp1*window[o1];
328 out[o2]=temp1*window[o2];
329 out[o3]=temp2*window[o3];
330 out[o4]=temp2*window[o4];