1 /********************************************************************
3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
5 * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH *
6 * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
8 * THE OggVorbis 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.18 2000/11/06 00:07:01 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=_ogg_malloc(sizeof(int)*(n/4));
49 float *trig=_ogg_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;
115 w2[i]= *xA++ + *xB++;
118 w[i++]= x0 * A[0] + x1 * A[1];
119 w[i]= x1 * A[0] - x0 * A[1];
120 w2[i++]= *xA++ + *xB++;
128 for(i=0;i<init->log2n-3;i++){
135 for(r=0;r<(k0>>2);r++){
151 wA =w[++w1] -w[++w2];
153 x[w2] =wA*AEv - wB*AOv;
154 x[w2-1]=wB*AEv + wA*AOv;
159 wA =w[++w1] -w[++w2];
161 x[w2] =wA*AEv - wB*AOv;
162 x[w2-1]=wB*AEv + wA*AOv;
167 wA =w[++w1] -w[++w2];
169 x[w2] =wA*AEv - wB*AOv;
170 x[w2-1]=wB*AEv + wA*AOv;
175 wA =w[++w1] -w[++w2];
177 x[w2] =wA*AEv - wB*AOv;
178 x[w2-1]=wB*AEv + wA*AOv;
187 wA =w[++w1] -w[++w2];
189 x[w2] =wA*AEv - wB*AOv;
190 x[w2-1]=wB*AEv + wA*AOv;
206 /* step 4, 5, 6, 7 */
208 float *C=init->trig+n;
209 int *bit=init->bitrev;
217 float wA=w[t1]-w[t2+1];
218 float wB=w[t1-1]+w[t2];
219 float wC=w[t1]+w[t2+1];
220 float wD=w[t1-1]-w[t2];
227 *x1++=( wC+wACO+wBCE)*.5;
228 *x2--=(-wD+wBCO-wACE)*.5;
229 *x1++=( wD+wBCO-wACE)*.5;
230 *x2--=( wC-wACO-wBCE)*.5;
236 void mdct_forward(mdct_lookup *init, float *in, float *out){
238 float *x=alloca(sizeof(float)*(n/2));
239 float *w=alloca(sizeof(float)*(n/2));
246 /* window + rotate + step 1 */
251 float *A=init->trig+n2;
257 tempA= in[in1+2] + in[in2];
258 tempB= in[in1] + in[in2+2];
260 x[i]= tempB*A[1] + tempA*A[0];
261 x[i+1]= tempB*A[0] - tempA*A[1];
268 tempA= in[in1+2] - in[in2];
269 tempB= in[in1] - in[in2+2];
271 x[i]= tempB*A[1] + tempA*A[0];
272 x[i+1]= tempB*A[0] - tempA*A[1];
279 tempA= -in[in1+2] - in[in2];
280 tempB= -in[in1] - in[in2+2];
282 x[i]= tempB*A[1] + tempA*A[0];
283 x[i+1]= tempB*A[0] - tempA*A[1];
287 xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
292 float *B=init->trig+n2;
296 out[i] =(xx[0]*B[0]+xx[1]*B[1])*scale;
297 *(--out2)=(xx[0]*B[1]-xx[1]*B[0])*scale;
305 void mdct_backward(mdct_lookup *init, float *in, float *out){
307 float *x=alloca(sizeof(float)*(n/2));
308 float *w=alloca(sizeof(float)*(n/2));
315 /* rotate + step 1 */
319 float *A=init->trig+n2;
323 *xO++=-*(inO+2)*A[1] - *inO*A[0];
324 *xO++= *inO*A[1] - *(inO+2)*A[0];
332 *xO++=*inO*A[1] + *(inO+2)*A[0];
333 *xO++=*inO*A[0] - *(inO+2)*A[1];
339 xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
344 float *B=init->trig+n2;
346 int o3=n4+n2,o4=o3-1;
349 float temp1= (*xx * B[1] - *(xx+1) * B[0]);
350 float temp2=-(*xx * B[0] + *(xx+1) * B[1]);