There's nothing, but it now compiles
[opus.git] / libcelt / mdct.c
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-2002             *
9  * by the XIPHOPHORUS Company 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 7187 2004-07-20 07:24:27Z 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 #define STIN static inline
50
51 /* build lookups for trig functions; also pre-figure scaling and
52    some window function algebra. */
53
54 void mdct_init(mdct_lookup *lookup,int n){
55   int   *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4));
56   DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4));
57   
58   int i;
59   int n2=n>>1;
60   int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
61   lookup->n=n;
62   lookup->trig=T;
63   lookup->bitrev=bitrev;
64
65 /* trig lookups... */
66
67   for(i=0;i<n/4;i++){
68     T[i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i)));
69     T[i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i)));
70     T[n2+i*2]=FLOAT_CONV(cos((M_PI/(2*n))*(2*i+1)));
71     T[n2+i*2+1]=FLOAT_CONV(sin((M_PI/(2*n))*(2*i+1)));
72   }
73   for(i=0;i<n/8;i++){
74     T[n+i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i+2))*.5);
75     T[n+i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i+2))*.5);
76   }
77
78   /* bitreverse lookup... */
79
80   {
81     int mask=(1<<(log2n-1))-1,i,j;
82     int msb=1<<(log2n-2);
83     for(i=0;i<n/8;i++){
84       int acc=0;
85       for(j=0;msb>>j;j++)
86         if((msb>>j)&i)acc|=1<<j;
87       bitrev[i*2]=((~acc)&mask)-1;
88       bitrev[i*2+1]=acc;
89
90     }
91   }
92   lookup->scale=FLOAT_CONV(4.f/n);
93 }
94
95 /* 8 point butterfly (in place, 4 register) */
96 STIN void mdct_butterfly_8(DATA_TYPE *x){
97   REG_TYPE r0   = x[6] + x[2];
98   REG_TYPE r1   = x[6] - x[2];
99   REG_TYPE r2   = x[4] + x[0];
100   REG_TYPE r3   = x[4] - x[0];
101
102            x[6] = r0   + r2;
103            x[4] = r0   - r2;
104            
105            r0   = x[5] - x[1];
106            r2   = x[7] - x[3];
107            x[0] = r1   + r0;
108            x[2] = r1   - r0;
109            
110            r0   = x[5] + x[1];
111            r1   = x[7] + x[3];
112            x[3] = r2   + r3;
113            x[1] = r2   - r3;
114            x[7] = r1   + r0;
115            x[5] = r1   - r0;
116            
117 }
118
119 /* 16 point butterfly (in place, 4 register) */
120 STIN void mdct_butterfly_16(DATA_TYPE *x){
121   REG_TYPE r0     = x[1]  - x[9];
122   REG_TYPE r1     = x[0]  - x[8];
123
124            x[8]  += x[0];
125            x[9]  += x[1];
126            x[0]   = MULT_NORM((r0   + r1) * cPI2_8);
127            x[1]   = MULT_NORM((r0   - r1) * cPI2_8);
128
129            r0     = x[3]  - x[11];
130            r1     = x[10] - x[2];
131            x[10] += x[2];
132            x[11] += x[3];
133            x[2]   = r0;
134            x[3]   = r1;
135
136            r0     = x[12] - x[4];
137            r1     = x[13] - x[5];
138            x[12] += x[4];
139            x[13] += x[5];
140            x[4]   = MULT_NORM((r0   - r1) * cPI2_8);
141            x[5]   = MULT_NORM((r0   + r1) * cPI2_8);
142
143            r0     = x[14] - x[6];
144            r1     = x[15] - x[7];
145            x[14] += x[6];
146            x[15] += x[7];
147            x[6]  = r0;
148            x[7]  = r1;
149
150            mdct_butterfly_8(x);
151            mdct_butterfly_8(x+8);
152 }
153
154 /* 32 point butterfly (in place, 4 register) */
155 STIN void mdct_butterfly_32(DATA_TYPE *x){
156   REG_TYPE r0     = x[30] - x[14];
157   REG_TYPE r1     = x[31] - x[15];
158
159            x[30] +=         x[14];           
160            x[31] +=         x[15];
161            x[14]  =         r0;              
162            x[15]  =         r1;
163
164            r0     = x[28] - x[12];   
165            r1     = x[29] - x[13];
166            x[28] +=         x[12];           
167            x[29] +=         x[13];
168            x[12]  = MULT_NORM( r0 * cPI1_8  -  r1 * cPI3_8 );
169            x[13]  = MULT_NORM( r0 * cPI3_8  +  r1 * cPI1_8 );
170
171            r0     = x[26] - x[10];
172            r1     = x[27] - x[11];
173            x[26] +=         x[10];
174            x[27] +=         x[11];
175            x[10]  = MULT_NORM(( r0  - r1 ) * cPI2_8);
176            x[11]  = MULT_NORM(( r0  + r1 ) * cPI2_8);
177
178            r0     = x[24] - x[8];
179            r1     = x[25] - x[9];
180            x[24] += x[8];
181            x[25] += x[9];
182            x[8]   = MULT_NORM( r0 * cPI3_8  -  r1 * cPI1_8 );
183            x[9]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );
184
185            r0     = x[22] - x[6];
186            r1     = x[7]  - x[23];
187            x[22] += x[6];
188            x[23] += x[7];
189            x[6]   = r1;
190            x[7]   = r0;
191
192            r0     = x[4]  - x[20];
193            r1     = x[5]  - x[21];
194            x[20] += x[4];
195            x[21] += x[5];
196            x[4]   = MULT_NORM( r1 * cPI1_8  +  r0 * cPI3_8 );
197            x[5]   = MULT_NORM( r1 * cPI3_8  -  r0 * cPI1_8 );
198
199            r0     = x[2]  - x[18];
200            r1     = x[3]  - x[19];
201            x[18] += x[2];
202            x[19] += x[3];
203            x[2]   = MULT_NORM(( r1  + r0 ) * cPI2_8);
204            x[3]   = MULT_NORM(( r1  - r0 ) * cPI2_8);
205
206            r0     = x[0]  - x[16];
207            r1     = x[1]  - x[17];
208            x[16] += x[0];
209            x[17] += x[1];
210            x[0]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );
211            x[1]   = MULT_NORM( r1 * cPI1_8  -  r0 * cPI3_8 );
212
213            mdct_butterfly_16(x);
214            mdct_butterfly_16(x+16);
215
216 }
217
218 /* N point first stage butterfly (in place, 2 register) */
219 STIN void mdct_butterfly_first(DATA_TYPE *T,
220                                         DATA_TYPE *x,
221                                         int points){
222   
223   DATA_TYPE *x1        = x          + points      - 8;
224   DATA_TYPE *x2        = x          + (points>>1) - 8;
225   REG_TYPE   r0;
226   REG_TYPE   r1;
227
228   do{
229     
230                r0      = x1[6]      -  x2[6];
231                r1      = x1[7]      -  x2[7];
232                x1[6]  += x2[6];
233                x1[7]  += x2[7];
234                x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
235                x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
236                
237                r0      = x1[4]      -  x2[4];
238                r1      = x1[5]      -  x2[5];
239                x1[4]  += x2[4];
240                x1[5]  += x2[5];
241                x2[4]   = MULT_NORM(r1 * T[5]  +  r0 * T[4]);
242                x2[5]   = MULT_NORM(r1 * T[4]  -  r0 * T[5]);
243                
244                r0      = x1[2]      -  x2[2];
245                r1      = x1[3]      -  x2[3];
246                x1[2]  += x2[2];
247                x1[3]  += x2[3];
248                x2[2]   = MULT_NORM(r1 * T[9]  +  r0 * T[8]);
249                x2[3]   = MULT_NORM(r1 * T[8]  -  r0 * T[9]);
250                
251                r0      = x1[0]      -  x2[0];
252                r1      = x1[1]      -  x2[1];
253                x1[0]  += x2[0];
254                x1[1]  += x2[1];
255                x2[0]   = MULT_NORM(r1 * T[13] +  r0 * T[12]);
256                x2[1]   = MULT_NORM(r1 * T[12] -  r0 * T[13]);
257                
258     x1-=8;
259     x2-=8;
260     T+=16;
261
262   }while(x2>=x);
263 }
264
265 /* N/stage point generic N stage butterfly (in place, 2 register) */
266 STIN void mdct_butterfly_generic(DATA_TYPE *T,
267                                           DATA_TYPE *x,
268                                           int points,
269                                           int trigint){
270   
271   DATA_TYPE *x1        = x          + points      - 8;
272   DATA_TYPE *x2        = x          + (points>>1) - 8;
273   REG_TYPE   r0;
274   REG_TYPE   r1;
275
276   do{
277     
278                r0      = x1[6]      -  x2[6];
279                r1      = x1[7]      -  x2[7];
280                x1[6]  += x2[6];
281                x1[7]  += x2[7];
282                x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
283                x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
284                
285                T+=trigint;
286                
287                r0      = x1[4]      -  x2[4];
288                r1      = x1[5]      -  x2[5];
289                x1[4]  += x2[4];
290                x1[5]  += x2[5];
291                x2[4]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
292                x2[5]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
293                
294                T+=trigint;
295                
296                r0      = x1[2]      -  x2[2];
297                r1      = x1[3]      -  x2[3];
298                x1[2]  += x2[2];
299                x1[3]  += x2[3];
300                x2[2]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
301                x2[3]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
302                
303                T+=trigint;
304                
305                r0      = x1[0]      -  x2[0];
306                r1      = x1[1]      -  x2[1];
307                x1[0]  += x2[0];
308                x1[1]  += x2[1];
309                x2[0]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
310                x2[1]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
311
312                T+=trigint;
313     x1-=8;
314     x2-=8;
315
316   }while(x2>=x);
317 }
318
319 STIN void mdct_butterflies(mdct_lookup *init,
320                              DATA_TYPE *x,
321                              int points){
322   
323   DATA_TYPE *T=init->trig;
324   int stages=init->log2n-5;
325   int i,j;
326   
327   if(--stages>0){
328     mdct_butterfly_first(T,x,points);
329   }
330
331   for(i=1;--stages>0;i++){
332     for(j=0;j<(1<<i);j++)
333       mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
334   }
335
336   for(j=0;j<points;j+=32)
337     mdct_butterfly_32(x+j);
338
339 }
340
341 void mdct_clear(mdct_lookup *l){
342   if(l){
343     if(l->trig)_ogg_free(l->trig);
344     if(l->bitrev)_ogg_free(l->bitrev);
345     memset(l,0,sizeof(*l));
346   }
347 }
348
349 STIN void mdct_bitreverse(mdct_lookup *init, 
350                             DATA_TYPE *x){
351   int        n       = init->n;
352   int       *bit     = init->bitrev;
353   DATA_TYPE *w0      = x;
354   DATA_TYPE *w1      = x = w0+(n>>1);
355   DATA_TYPE *T       = init->trig+n;
356
357   do{
358     DATA_TYPE *x0    = x+bit[0];
359     DATA_TYPE *x1    = x+bit[1];
360
361     REG_TYPE  r0     = x0[1]  - x1[1];
362     REG_TYPE  r1     = x0[0]  + x1[0];
363     REG_TYPE  r2     = MULT_NORM(r1     * T[0]   + r0 * T[1]);
364     REG_TYPE  r3     = MULT_NORM(r1     * T[1]   - r0 * T[0]);
365
366               w1    -= 4;
367
368               r0     = HALVE(x0[1] + x1[1]);
369               r1     = HALVE(x0[0] - x1[0]);
370       
371               w0[0]  = r0     + r2;
372               w1[2]  = r0     - r2;
373               w0[1]  = r1     + r3;
374               w1[3]  = r3     - r1;
375
376               x0     = x+bit[2];
377               x1     = x+bit[3];
378
379               r0     = x0[1]  - x1[1];
380               r1     = x0[0]  + x1[0];
381               r2     = MULT_NORM(r1     * T[2]   + r0 * T[3]);
382               r3     = MULT_NORM(r1     * T[3]   - r0 * T[2]);
383
384               r0     = HALVE(x0[1] + x1[1]);
385               r1     = HALVE(x0[0] - x1[0]);
386       
387               w0[2]  = r0     + r2;
388               w1[0]  = r0     - r2;
389               w0[3]  = r1     + r3;
390               w1[1]  = r3     - r1;
391
392               T     += 4;
393               bit   += 4;
394               w0    += 4;
395
396   }while(w0<w1);
397 }
398
399 void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
400   int n=init->n;
401   int n2=n>>1;
402   int n4=n>>2;
403
404   /* rotate */
405
406   DATA_TYPE *iX = in+n2-7;
407   DATA_TYPE *oX = out+n2+n4;
408   DATA_TYPE *T  = init->trig+n4;
409
410   do{
411     oX         -= 4;
412     oX[0]       = MULT_NORM(-iX[2] * T[3] - iX[0]  * T[2]);
413     oX[1]       = MULT_NORM (iX[0] * T[3] - iX[2]  * T[2]);
414     oX[2]       = MULT_NORM(-iX[6] * T[1] - iX[4]  * T[0]);
415     oX[3]       = MULT_NORM (iX[4] * T[1] - iX[6]  * T[0]);
416     iX         -= 8;
417     T          += 4;
418   }while(iX>=in);
419
420   iX            = in+n2-8;
421   oX            = out+n2+n4;
422   T             = init->trig+n4;
423
424   do{
425     T          -= 4;
426     oX[0]       =  MULT_NORM (iX[4] * T[3] + iX[6] * T[2]);
427     oX[1]       =  MULT_NORM (iX[4] * T[2] - iX[6] * T[3]);
428     oX[2]       =  MULT_NORM (iX[0] * T[1] + iX[2] * T[0]);
429     oX[3]       =  MULT_NORM (iX[0] * T[0] - iX[2] * T[1]);
430     iX         -= 8;
431     oX         += 4;
432   }while(iX>=in);
433
434   mdct_butterflies(init,out+n2,n2);
435   mdct_bitreverse(init,out);
436
437   /* roatate + window */
438
439   {
440     DATA_TYPE *oX1=out+n2+n4;
441     DATA_TYPE *oX2=out+n2+n4;
442     DATA_TYPE *iX =out;
443     T             =init->trig+n2;
444     
445     do{
446       oX1-=4;
447
448       oX1[3]  =  MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
449       oX2[0]  = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
450
451       oX1[2]  =  MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
452       oX2[1]  = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
453
454       oX1[1]  =  MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
455       oX2[2]  = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
456
457       oX1[0]  =  MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
458       oX2[3]  = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
459
460       oX2+=4;
461       iX    +=   8;
462       T     +=   8;
463     }while(iX<oX1);
464
465     iX=out+n2+n4;
466     oX1=out+n4;
467     oX2=oX1;
468
469     do{
470       oX1-=4;
471       iX-=4;
472
473       oX2[0] = -(oX1[3] = iX[3]);
474       oX2[1] = -(oX1[2] = iX[2]);
475       oX2[2] = -(oX1[1] = iX[1]);
476       oX2[3] = -(oX1[0] = iX[0]);
477
478       oX2+=4;
479     }while(oX2<iX);
480
481     iX=out+n2+n4;
482     oX1=out+n2+n4;
483     oX2=out+n2;
484     do{
485       oX1-=4;
486       oX1[0]= iX[3];
487       oX1[1]= iX[2];
488       oX1[2]= iX[1];
489       oX1[3]= iX[0];
490       iX+=4;
491     }while(oX1>oX2);
492   }
493 }
494
495 void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
496   int n=init->n;
497   int n2=n>>1;
498   int n4=n>>2;
499   int n8=n>>3;
500   DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */
501   DATA_TYPE *w2=w+n2;
502
503   /* rotate */
504
505   /* window + rotate + step 1 */
506   
507   REG_TYPE r0;
508   REG_TYPE r1;
509   DATA_TYPE *x0=in+n2+n4;
510   DATA_TYPE *x1=x0+1;
511   DATA_TYPE *T=init->trig+n2;
512   
513   int i=0;
514   
515   for(i=0;i<n8;i+=2){
516     x0 -=4;
517     T-=2;
518     r0= x0[2] + x1[0];
519     r1= x0[0] + x1[2];       
520     w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
521     w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
522     x1 +=4;
523   }
524
525   x1=in+1;
526   
527   for(;i<n2-n8;i+=2){
528     T-=2;
529     x0 -=4;
530     r0= x0[2] - x1[0];
531     r1= x0[0] - x1[2];       
532     w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
533     w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
534     x1 +=4;
535   }
536     
537   x0=in+n;
538
539   for(;i<n2;i+=2){
540     T-=2;
541     x0 -=4;
542     r0= -x0[2] - x1[0];
543     r1= -x0[0] - x1[2];       
544     w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
545     w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
546     x1 +=4;
547   }
548
549
550   mdct_butterflies(init,w+n2,n2);
551   mdct_bitreverse(init,w);
552
553   /* roatate + window */
554
555   T=init->trig+n2;
556   x0=out+n2;
557
558   for(i=0;i<n4;i++){
559     x0--;
560     out[i] =MULT_NORM((w[0]*T[0]+w[1]*T[1])*init->scale);
561     x0[0]  =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale);
562     w+=2;
563     T+=2;
564   }
565 }
566