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 ********************************************************************/
44 static double oPI = 3.14159265358979323846;
46 /* build lookups for trig functions; also pre-figure scaling and
47 some window function algebra. */
49 MDCT_lookup *MDCT_init(int n){
50 MDCT_lookup *lookup=malloc(sizeof(MDCT_lookup));
51 int *bitrev=malloc(sizeof(int)*(n/4));
52 double *trig=malloc(sizeof(double)*(n+n/4));
64 int log2n=lookup->log2n=rint(log(n)/log(2));
67 lookup->bitrev=bitrev;
72 AE[i]=cos((oPI/n)*(4*i));
73 AO[i]=-sin((oPI/n)*(4*i));
74 BE[i]=cos((oPI/(2*n))*(2*i+1));
75 BO[i]=sin((oPI/(2*n))*(2*i+1));
78 CE[i]=cos((oPI/n)*(4*i+2));
79 CO[i]=-sin((oPI/n)*(4*i+2));
82 /* bitreverse lookup... */
85 int mask=(1<<(log2n-1))-1,i,j;
90 if((msb>>j)&i)acc|=1<<j;
91 bitA[i]=((~acc)&mask)*2;
99 void MDCT_free(MDCT_lookup *l){
101 if(l->trig)free(l->trig);
102 if(l->bitrev)free(l->bitrev);
107 static inline void _MDCT_kernel(double *x,
108 int n, int n2, int n4, int n8,
110 double *w=x+1; /* interleaved access improves cache locality */
118 double *AE=init->trig+n4;
122 double x0=xA[i]-xB[i];
123 double x1=xA[i+2]-xB[i+2];
126 w[i] =x0 * *AE + x1 * *AO;
129 w[i] =x1 * *AE - x0 * *AO;
139 for(i=0;i<init->log2n-3;i++){
143 double *AE=init->trig;
147 for(r=0;r<(n4>>i);r+=4){
149 int w2=wbase-(k0>>1);
152 for(s=0;s<(2<<i);s++){
153 x[w1+2]=w[w1+2]+w[w2+2];
155 x[w2+2]=(w[w1+2]-w[w2+2])* *AE-(w[w1]-w[w2])* *AO;
156 x[w2] =(w[w1]-w[w2])* *AE+(w[w1+2]-w[w2+2])* *AO;
171 /* step 4, 5, 6, 7 */
173 double *CE=init->trig+n;
175 int *bitA=init->bitrev;
185 double wA=w[t1]-w[t2];
186 double wB=w[t3]+w[t4];
187 double wC=w[t1]+w[t2];
188 double wD=w[t3]-w[t4];
191 double wBCO=wB* *(CO++);
193 double wBCE=wB* *(CE++);
195 *x1 =( wC+wACO+wBCE)*.5;
196 *(x2-2)=( wC-wACO-wBCE)*.5;
197 *(x1+2)=( wD+wBCO-wACE)*.5;
198 *x2 =(-wD+wBCO-wACE)*.5;
205 void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
207 double *x=alloca(n*sizeof(double));
213 /* window + rotate + step 1 */
218 double *AE=init->trig+n4;
224 tempA= in[in1+2]*window[in1+2] + in[in2]*window[in2];
225 tempB= in[in1]*window[in1] + in[in2+2]*window[in2+2];
228 x[i]= tempB* *AO + tempA* *AE;
229 x[i+2]= tempB* *AE - tempA* *AO;
235 tempA= in[in1+2]*window[in1+2] - in[in2]*window[in2];
236 tempB= in[in1]*window[in1] - in[in2+2]*window[in2+2];
239 x[i]= tempB* *AO + tempA* *AE;
240 x[i+2]= tempB* *AE - tempA* *AO;
246 tempA= -in[in1+2]*window[in1+2] - in[in2]*window[in2];
247 tempB= -in[in1]*window[in1] - in[in2+2]*window[in2+2];
250 x[i]= tempB* *AO + tempA* *AE;
251 x[i+2]= tempB* *AE - tempA* *AO;
255 _MDCT_kernel(x,n,n2,n4,n8,init);
260 double *BE=init->trig+n2;
264 out[i] =x[0]* *BE+x[2]* *BO;
265 *(--out2)=x[0]* *BO-x[2]* *BE;
273 void iMDCT(double *in, double *out, MDCT_lookup *init, double *window){
275 double *x=alloca(n*sizeof(double));
281 /* window + rotate + step 1 */
285 double *AE=init->trig+n4;
290 *xO=-*(inO+2)* *AO - *inO * *AE;
292 *xO= *inO * *AO - *(inO+2)* *AE;
301 *xO=*inO * *AO + *(inO+2) * *AE;
303 *xO=*inO * *AE - *(inO+2) * *AO;
310 _MDCT_kernel(x,n,n2,n4,n8,init);
315 double *BE=init->trig+n2;
318 int o3=n4+n2,o4=o3-1;
322 double temp1= (*x * *BO - *(x+2) * *BE)* scale;
323 double temp2= (*x * *BE + *(x+2) * *BO)* -scale;
325 out[o1]=-temp1*window[o1];
326 out[o2]=temp1*window[o2];
327 out[o3]=temp2*window[o3];
328 out[o4]=temp2*window[o4];