mdct.c 15.1 KB
Newer Older
Mike Whitson's avatar
Mike Whitson committed
1 2
/********************************************************************
 *                                                                  *
Monty's avatar
 
Monty committed
3
 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
Monty's avatar
 
Monty committed
4 5 6
 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
Mike Whitson's avatar
Mike Whitson committed
7
 *                                                                  *
Monty's avatar
Monty committed
8
 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
9
 * by the Xiph.Org Foundation http://www.xiph.org/                  *
10
 *                                                                  *
Mike Whitson's avatar
Mike Whitson committed
11 12
 ********************************************************************

Monty's avatar
 
Monty committed
13
 function: normalized modified discrete cosine transform
Monty's avatar
 
Monty committed
14
           power of two length transform only [64 <= n ]
Mike Whitson's avatar
Mike Whitson committed
15

Monty's avatar
 
Monty committed
16 17 18 19 20
 Original algorithm adapted long ago from _The use of multirate filter
 banks for coding of high quality digital audio_, by T. Sporer,
 K. Brandenburg and B. Edler, collection of the European Signal
 Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp
 211-214
Mike Whitson's avatar
Mike Whitson committed
21

Monty's avatar
 
Monty committed
22 23 24
 The below code implements an algorithm that no longer looks much like
 that presented in the paper, but the basic structure remains if you
 dig deep enough to see it.
Mike Whitson's avatar
Mike Whitson committed
25

Monty's avatar
 
Monty committed
26 27
 This module DOES NOT INCLUDE code to generate/apply the window
 function.  Everybody has their own weird favorite including me... I
Monty's avatar
 
Monty committed
28
 happen to like the properties of y=sin(.5PI*sin^2(x)), but others may
Monty's avatar
 
Monty committed
29
 vehemently disagree.
Mike Whitson's avatar
Mike Whitson committed
30 31 32

 ********************************************************************/

Monty's avatar
 
Monty committed
33 34 35 36 37 38
/* this can also be run as an integer transform by uncommenting a
   define in mdct.h; the integerization is a first pass and although
   it's likely stable for Vorbis, the dynamic range is constrained and
   roundoff isn't done (so it's noisy).  Consider it functional, but
   only a starting point.  There's no point on a machine with an FPU */

Monty's avatar
 
Monty committed
39
#include <stdio.h>
Mike Whitson's avatar
Mike Whitson committed
40
#include <stdlib.h>
Monty's avatar
 
Monty committed
41
#include <string.h>
Mike Whitson's avatar
Mike Whitson committed
42
#include <math.h>
Monty's avatar
 
Monty committed
43
#include "vorbis/codec.h"
Mike Whitson's avatar
Mike Whitson committed
44
#include "mdct.h"
Monty's avatar
 
Monty committed
45
#include "os.h"
Monty's avatar
 
Monty committed
46
#include "misc.h"
Mike Whitson's avatar
Mike Whitson committed
47 48 49 50

/* build lookups for trig functions; also pre-figure scaling and
   some window function algebra. */

Monty's avatar
 
Monty committed
51
void mdct_init(mdct_lookup *lookup,int n){
52 53
  int   *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4));
  DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4));
Monty's avatar
Monty committed
54

Mike Whitson's avatar
Mike Whitson committed
55
  int i;
Monty's avatar
 
Monty committed
56
  int n2=n>>1;
57
  int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
Mike Whitson's avatar
Mike Whitson committed
58
  lookup->n=n;
Monty's avatar
 
Monty committed
59
  lookup->trig=T;
Mike Whitson's avatar
Mike Whitson committed
60 61
  lookup->bitrev=bitrev;

Monty's avatar
 
Monty committed
62
/* trig lookups... */
Mike Whitson's avatar
Mike Whitson committed
63 64

  for(i=0;i<n/4;i++){
Monty's avatar
 
Monty committed
65 66 67 68
    T[i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i)));
    T[i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i)));
    T[n2+i*2]=FLOAT_CONV(cos((M_PI/(2*n))*(2*i+1)));
    T[n2+i*2+1]=FLOAT_CONV(sin((M_PI/(2*n))*(2*i+1)));
Mike Whitson's avatar
Mike Whitson committed
69 70
  }
  for(i=0;i<n/8;i++){
Monty's avatar
 
Monty committed
71 72
    T[n+i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i+2))*.5);
    T[n+i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i+2))*.5);
Mike Whitson's avatar
Mike Whitson committed
73 74 75 76 77 78 79 80 81 82
  }

  /* bitreverse lookup... */

  {
    int mask=(1<<(log2n-1))-1,i,j;
    int msb=1<<(log2n-2);
    for(i=0;i<n/8;i++){
      int acc=0;
      for(j=0;msb>>j;j++)
83
        if((msb>>j)&i)acc|=1<<j;
Monty's avatar
 
Monty committed
84
      bitrev[i*2]=((~acc)&mask)-1;
Monty's avatar
 
Monty committed
85
      bitrev[i*2+1]=acc;
Monty's avatar
 
Monty committed
86

Mike Whitson's avatar
Mike Whitson committed
87 88
    }
  }
Monty's avatar
 
Monty committed
89
  lookup->scale=FLOAT_CONV(4.f/n);
Monty's avatar
 
Monty committed
90 91 92
}

/* 8 point butterfly (in place, 4 register) */
Monty's avatar
 
Monty committed
93
STIN void mdct_butterfly_8(DATA_TYPE *x){
Monty's avatar
 
Monty committed
94 95 96 97 98
  REG_TYPE r0   = x[6] + x[2];
  REG_TYPE r1   = x[6] - x[2];
  REG_TYPE r2   = x[4] + x[0];
  REG_TYPE r3   = x[4] - x[0];

99 100
           x[6] = r0   + r2;
           x[4] = r0   - r2;
Monty's avatar
Monty committed
101

102 103 104 105
           r0   = x[5] - x[1];
           r2   = x[7] - x[3];
           x[0] = r1   + r0;
           x[2] = r1   - r0;
Monty's avatar
Monty committed
106

107 108 109 110 111 112
           r0   = x[5] + x[1];
           r1   = x[7] + x[3];
           x[3] = r2   + r3;
           x[1] = r2   - r3;
           x[7] = r1   + r0;
           x[5] = r1   - r0;
Monty's avatar
Monty committed
113

Monty's avatar
 
Monty committed
114 115 116
}

/* 16 point butterfly (in place, 4 register) */
Monty's avatar
 
Monty committed
117
STIN void mdct_butterfly_16(DATA_TYPE *x){
Monty's avatar
 
Monty committed
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
  REG_TYPE r0     = x[1]  - x[9];
  REG_TYPE r1     = x[0]  - x[8];

           x[8]  += x[0];
           x[9]  += x[1];
           x[0]   = MULT_NORM((r0   + r1) * cPI2_8);
           x[1]   = MULT_NORM((r0   - r1) * cPI2_8);

           r0     = x[3]  - x[11];
           r1     = x[10] - x[2];
           x[10] += x[2];
           x[11] += x[3];
           x[2]   = r0;
           x[3]   = r1;

           r0     = x[12] - x[4];
           r1     = x[13] - x[5];
           x[12] += x[4];
           x[13] += x[5];
           x[4]   = MULT_NORM((r0   - r1) * cPI2_8);
           x[5]   = MULT_NORM((r0   + r1) * cPI2_8);

           r0     = x[14] - x[6];
           r1     = x[15] - x[7];
           x[14] += x[6];
           x[15] += x[7];
           x[6]  = r0;
           x[7]  = r1;

147 148
           mdct_butterfly_8(x);
           mdct_butterfly_8(x+8);
Monty's avatar
 
Monty committed
149 150 151
}

/* 32 point butterfly (in place, 4 register) */
Monty's avatar
 
Monty committed
152
STIN void mdct_butterfly_32(DATA_TYPE *x){
Monty's avatar
 
Monty committed
153 154 155
  REG_TYPE r0     = x[30] - x[14];
  REG_TYPE r1     = x[31] - x[15];

Monty's avatar
Monty committed
156
           x[30] +=         x[14];
157
           x[31] +=         x[15];
Monty's avatar
Monty committed
158
           x[14]  =         r0;
159
           x[15]  =         r1;
Monty's avatar
 
Monty committed
160

Monty's avatar
Monty committed
161
           r0     = x[28] - x[12];
162
           r1     = x[29] - x[13];
Monty's avatar
Monty committed
163
           x[28] +=         x[12];
164
           x[29] +=         x[13];
Monty's avatar
 
Monty committed
165
           x[12]  = MULT_NORM( r0 * cPI1_8  -  r1 * cPI3_8 );
166
           x[13]  = MULT_NORM( r0 * cPI3_8  +  r1 * cPI1_8 );
Monty's avatar
 
Monty committed
167 168

           r0     = x[26] - x[10];
169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
           r1     = x[27] - x[11];
           x[26] +=         x[10];
           x[27] +=         x[11];
           x[10]  = MULT_NORM(( r0  - r1 ) * cPI2_8);
           x[11]  = MULT_NORM(( r0  + r1 ) * cPI2_8);

           r0     = x[24] - x[8];
           r1     = x[25] - x[9];
           x[24] += x[8];
           x[25] += x[9];
           x[8]   = MULT_NORM( r0 * cPI3_8  -  r1 * cPI1_8 );
           x[9]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );

           r0     = x[22] - x[6];
           r1     = x[7]  - x[23];
           x[22] += x[6];
           x[23] += x[7];
           x[6]   = r1;
           x[7]   = r0;

           r0     = x[4]  - x[20];
           r1     = x[5]  - x[21];
           x[20] += x[4];
           x[21] += x[5];
           x[4]   = MULT_NORM( r1 * cPI1_8  +  r0 * cPI3_8 );
           x[5]   = MULT_NORM( r1 * cPI3_8  -  r0 * cPI1_8 );

           r0     = x[2]  - x[18];
           r1     = x[3]  - x[19];
           x[18] += x[2];
           x[19] += x[3];
           x[2]   = MULT_NORM(( r1  + r0 ) * cPI2_8);
           x[3]   = MULT_NORM(( r1  - r0 ) * cPI2_8);

           r0     = x[0]  - x[16];
           r1     = x[1]  - x[17];
           x[16] += x[0];
           x[17] += x[1];
           x[0]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );
           x[1]   = MULT_NORM( r1 * cPI1_8  -  r0 * cPI3_8 );

           mdct_butterfly_16(x);
           mdct_butterfly_16(x+16);
Monty's avatar
 
Monty committed
212 213 214 215

}

/* N point first stage butterfly (in place, 2 register) */
Monty's avatar
 
Monty committed
216
STIN void mdct_butterfly_first(DATA_TYPE *T,
217 218
                                        DATA_TYPE *x,
                                        int points){
Monty's avatar
Monty committed
219

Monty's avatar
 
Monty committed
220 221 222 223 224 225
  DATA_TYPE *x1        = x          + points      - 8;
  DATA_TYPE *x2        = x          + (points>>1) - 8;
  REG_TYPE   r0;
  REG_TYPE   r1;

  do{
Monty's avatar
Monty committed
226

Monty's avatar
 
Monty committed
227
               r0      = x1[6]      -  x2[6];
228 229 230 231 232
               r1      = x1[7]      -  x2[7];
               x1[6]  += x2[6];
               x1[7]  += x2[7];
               x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
               x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
Monty's avatar
Monty committed
233

234 235 236 237 238 239
               r0      = x1[4]      -  x2[4];
               r1      = x1[5]      -  x2[5];
               x1[4]  += x2[4];
               x1[5]  += x2[5];
               x2[4]   = MULT_NORM(r1 * T[5]  +  r0 * T[4]);
               x2[5]   = MULT_NORM(r1 * T[4]  -  r0 * T[5]);
Monty's avatar
Monty committed
240

241 242 243 244 245 246
               r0      = x1[2]      -  x2[2];
               r1      = x1[3]      -  x2[3];
               x1[2]  += x2[2];
               x1[3]  += x2[3];
               x2[2]   = MULT_NORM(r1 * T[9]  +  r0 * T[8]);
               x2[3]   = MULT_NORM(r1 * T[8]  -  r0 * T[9]);
Monty's avatar
Monty committed
247

248 249 250 251 252 253
               r0      = x1[0]      -  x2[0];
               r1      = x1[1]      -  x2[1];
               x1[0]  += x2[0];
               x1[1]  += x2[1];
               x2[0]   = MULT_NORM(r1 * T[13] +  r0 * T[12]);
               x2[1]   = MULT_NORM(r1 * T[12] -  r0 * T[13]);
Monty's avatar
Monty committed
254

Monty's avatar
 
Monty committed
255 256 257 258 259 260 261 262
    x1-=8;
    x2-=8;
    T+=16;

  }while(x2>=x);
}

/* N/stage point generic N stage butterfly (in place, 2 register) */
Monty's avatar
 
Monty committed
263
STIN void mdct_butterfly_generic(DATA_TYPE *T,
264 265 266
                                          DATA_TYPE *x,
                                          int points,
                                          int trigint){
Monty's avatar
Monty committed
267

Monty's avatar
 
Monty committed
268 269 270 271 272 273
  DATA_TYPE *x1        = x          + points      - 8;
  DATA_TYPE *x2        = x          + (points>>1) - 8;
  REG_TYPE   r0;
  REG_TYPE   r1;

  do{
Monty's avatar
Monty committed
274

Monty's avatar
 
Monty committed
275
               r0      = x1[6]      -  x2[6];
276 277 278 279 280
               r1      = x1[7]      -  x2[7];
               x1[6]  += x2[6];
               x1[7]  += x2[7];
               x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
               x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
Monty's avatar
Monty committed
281

282
               T+=trigint;
Monty's avatar
Monty committed
283

284 285 286 287 288 289
               r0      = x1[4]      -  x2[4];
               r1      = x1[5]      -  x2[5];
               x1[4]  += x2[4];
               x1[5]  += x2[5];
               x2[4]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
               x2[5]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
Monty's avatar
Monty committed
290

291
               T+=trigint;
Monty's avatar
Monty committed
292

293 294 295 296 297 298
               r0      = x1[2]      -  x2[2];
               r1      = x1[3]      -  x2[3];
               x1[2]  += x2[2];
               x1[3]  += x2[3];
               x2[2]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
               x2[3]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
Monty's avatar
Monty committed
299

300
               T+=trigint;
Monty's avatar
Monty committed
301

302 303 304 305 306 307 308 309
               r0      = x1[0]      -  x2[0];
               r1      = x1[1]      -  x2[1];
               x1[0]  += x2[0];
               x1[1]  += x2[1];
               x2[0]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
               x2[1]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);

               T+=trigint;
Monty's avatar
 
Monty committed
310 311 312 313 314 315
    x1-=8;
    x2-=8;

  }while(x2>=x);
}

Monty's avatar
 
Monty committed
316
STIN void mdct_butterflies(mdct_lookup *init,
317 318
                             DATA_TYPE *x,
                             int points){
Monty's avatar
Monty committed
319

Monty's avatar
 
Monty committed
320 321 322
  DATA_TYPE *T=init->trig;
  int stages=init->log2n-5;
  int i,j;
Monty's avatar
Monty committed
323

Monty's avatar
 
Monty committed
324 325 326 327 328 329 330 331 332 333 334 335
  if(--stages>0){
    mdct_butterfly_first(T,x,points);
  }

  for(i=1;--stages>0;i++){
    for(j=0;j<(1<<i);j++)
      mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
  }

  for(j=0;j<points;j+=32)
    mdct_butterfly_32(x+j);

Mike Whitson's avatar
Mike Whitson committed
336 337
}

Monty's avatar
 
Monty committed
338
void mdct_clear(mdct_lookup *l){
Mike Whitson's avatar
Mike Whitson committed
339
  if(l){
Monty's avatar
 
Monty committed
340 341
    if(l->trig)_ogg_free(l->trig);
    if(l->bitrev)_ogg_free(l->bitrev);
342
    memset(l,0,sizeof(*l));
Mike Whitson's avatar
Mike Whitson committed
343 344 345
  }
}

Monty's avatar
Monty committed
346
STIN void mdct_bitreverse(mdct_lookup *init,
347
                            DATA_TYPE *x){
Monty's avatar
 
Monty committed
348 349 350 351 352
  int        n       = init->n;
  int       *bit     = init->bitrev;
  DATA_TYPE *w0      = x;
  DATA_TYPE *w1      = x = w0+(n>>1);
  DATA_TYPE *T       = init->trig+n;
Mike Whitson's avatar
Mike Whitson committed
353

Monty's avatar
 
Monty committed
354 355 356
  do{
    DATA_TYPE *x0    = x+bit[0];
    DATA_TYPE *x1    = x+bit[1];
Monty's avatar
 
Monty committed
357

Monty's avatar
 
Monty committed
358 359 360 361
    REG_TYPE  r0     = x0[1]  - x1[1];
    REG_TYPE  r1     = x0[0]  + x1[0];
    REG_TYPE  r2     = MULT_NORM(r1     * T[0]   + r0 * T[1]);
    REG_TYPE  r3     = MULT_NORM(r1     * T[1]   - r0 * T[0]);
Mike Whitson's avatar
Mike Whitson committed
362

363
              w1    -= 4;
Mike Whitson's avatar
Mike Whitson committed
364

Monty's avatar
 
Monty committed
365 366
              r0     = HALVE(x0[1] + x1[1]);
              r1     = HALVE(x0[0] - x1[0]);
Monty's avatar
Monty committed
367

368 369 370 371
              w0[0]  = r0     + r2;
              w1[2]  = r0     - r2;
              w0[1]  = r1     + r3;
              w1[3]  = r3     - r1;
Mike Whitson's avatar
Mike Whitson committed
372

Monty's avatar
 
Monty committed
373 374 375 376 377 378 379 380 381 382
              x0     = x+bit[2];
              x1     = x+bit[3];

              r0     = x0[1]  - x1[1];
              r1     = x0[0]  + x1[0];
              r2     = MULT_NORM(r1     * T[2]   + r0 * T[3]);
              r3     = MULT_NORM(r1     * T[3]   - r0 * T[2]);

              r0     = HALVE(x0[1] + x1[1]);
              r1     = HALVE(x0[0] - x1[0]);
Monty's avatar
Monty committed
383

384 385 386 387 388 389 390 391
              w0[2]  = r0     + r2;
              w1[0]  = r0     - r2;
              w0[3]  = r1     + r3;
              w1[1]  = r3     - r1;

              T     += 4;
              bit   += 4;
              w0    += 4;
Monty's avatar
 
Monty committed
392 393

  }while(w0<w1);
Mike Whitson's avatar
Mike Whitson committed
394 395
}

Monty's avatar
 
Monty committed
396
void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
Mike Whitson's avatar
Mike Whitson committed
397 398 399 400
  int n=init->n;
  int n2=n>>1;
  int n4=n>>2;

Monty's avatar
 
Monty committed
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434
  /* rotate */

  DATA_TYPE *iX = in+n2-7;
  DATA_TYPE *oX = out+n2+n4;
  DATA_TYPE *T  = init->trig+n4;

  do{
    oX         -= 4;
    oX[0]       = MULT_NORM(-iX[2] * T[3] - iX[0]  * T[2]);
    oX[1]       = MULT_NORM (iX[0] * T[3] - iX[2]  * T[2]);
    oX[2]       = MULT_NORM(-iX[6] * T[1] - iX[4]  * T[0]);
    oX[3]       = MULT_NORM (iX[4] * T[1] - iX[6]  * T[0]);
    iX         -= 8;
    T          += 4;
  }while(iX>=in);

  iX            = in+n2-8;
  oX            = out+n2+n4;
  T             = init->trig+n4;

  do{
    T          -= 4;
    oX[0]       =  MULT_NORM (iX[4] * T[3] + iX[6] * T[2]);
    oX[1]       =  MULT_NORM (iX[4] * T[2] - iX[6] * T[3]);
    oX[2]       =  MULT_NORM (iX[0] * T[1] + iX[2] * T[0]);
    oX[3]       =  MULT_NORM (iX[0] * T[0] - iX[2] * T[1]);
    iX         -= 8;
    oX         += 4;
  }while(iX>=in);

  mdct_butterflies(init,out+n2,n2);
  mdct_bitreverse(init,out);

  /* roatate + window */
Mike Whitson's avatar
Mike Whitson committed
435

Monty's avatar
 
Monty committed
436 437 438 439 440
  {
    DATA_TYPE *oX1=out+n2+n4;
    DATA_TYPE *oX2=out+n2+n4;
    DATA_TYPE *iX =out;
    T             =init->trig+n2;
Monty's avatar
Monty committed
441

Monty's avatar
 
Monty committed
442 443
    do{
      oX1-=4;
Mike Whitson's avatar
Mike Whitson committed
444

Monty's avatar
 
Monty committed
445 446
      oX1[3]  =  MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
      oX2[0]  = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
Mike Whitson's avatar
Mike Whitson committed
447

Monty's avatar
 
Monty committed
448 449
      oX1[2]  =  MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
      oX2[1]  = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
Mike Whitson's avatar
Mike Whitson committed
450

Monty's avatar
 
Monty committed
451 452
      oX1[1]  =  MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
      oX2[2]  = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
Mike Whitson's avatar
Mike Whitson committed
453

Monty's avatar
 
Monty committed
454 455
      oX1[0]  =  MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
      oX2[3]  = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
Mike Whitson's avatar
Mike Whitson committed
456

Monty's avatar
 
Monty committed
457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488
      oX2+=4;
      iX    +=   8;
      T     +=   8;
    }while(iX<oX1);

    iX=out+n2+n4;
    oX1=out+n4;
    oX2=oX1;

    do{
      oX1-=4;
      iX-=4;

      oX2[0] = -(oX1[3] = iX[3]);
      oX2[1] = -(oX1[2] = iX[2]);
      oX2[2] = -(oX1[1] = iX[1]);
      oX2[3] = -(oX1[0] = iX[0]);

      oX2+=4;
    }while(oX2<iX);

    iX=out+n2+n4;
    oX1=out+n2+n4;
    oX2=out+n2;
    do{
      oX1-=4;
      oX1[0]= iX[3];
      oX1[1]= iX[2];
      oX1[2]= iX[1];
      oX1[3]= iX[0];
      iX+=4;
    }while(oX1>oX2);
Mike Whitson's avatar
Mike Whitson committed
489 490 491
  }
}

Monty's avatar
 
Monty committed
492
void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
Mike Whitson's avatar
Mike Whitson committed
493 494 495 496
  int n=init->n;
  int n2=n>>1;
  int n4=n>>2;
  int n8=n>>3;
497
  DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */
Monty's avatar
 
Monty committed
498
  DATA_TYPE *w2=w+n2;
Mike Whitson's avatar
Mike Whitson committed
499

Monty's avatar
 
Monty committed
500
  /* rotate */
Mike Whitson's avatar
Mike Whitson committed
501

Monty's avatar
 
Monty committed
502
  /* window + rotate + step 1 */
Monty's avatar
Monty committed
503

Monty's avatar
 
Monty committed
504 505 506 507 508
  REG_TYPE r0;
  REG_TYPE r1;
  DATA_TYPE *x0=in+n2+n4;
  DATA_TYPE *x1=x0+1;
  DATA_TYPE *T=init->trig+n2;
Monty's avatar
Monty committed
509

Monty's avatar
 
Monty committed
510
  int i=0;
Monty's avatar
Monty committed
511

Monty's avatar
 
Monty committed
512 513 514 515
  for(i=0;i<n8;i+=2){
    x0 -=4;
    T-=2;
    r0= x0[2] + x1[0];
Monty's avatar
Monty committed
516
    r1= x0[0] + x1[2];
Monty's avatar
 
Monty committed
517 518 519 520
    w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
    w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
    x1 +=4;
  }
Mike Whitson's avatar
Mike Whitson committed
521

Monty's avatar
 
Monty committed
522
  x1=in+1;
Monty's avatar
Monty committed
523

Monty's avatar
 
Monty committed
524 525 526 527
  for(;i<n2-n8;i+=2){
    T-=2;
    x0 -=4;
    r0= x0[2] - x1[0];
Monty's avatar
Monty committed
528
    r1= x0[0] - x1[2];
Monty's avatar
 
Monty committed
529 530 531 532
    w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
    w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
    x1 +=4;
  }
Monty's avatar
Monty committed
533

Monty's avatar
 
Monty committed
534 535 536 537 538 539
  x0=in+n;

  for(;i<n2;i+=2){
    T-=2;
    x0 -=4;
    r0= -x0[2] - x1[0];
Monty's avatar
Monty committed
540
    r1= -x0[0] - x1[2];
Monty's avatar
 
Monty committed
541 542 543
    w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
    w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
    x1 +=4;
Mike Whitson's avatar
Mike Whitson committed
544 545 546
  }


Monty's avatar
 
Monty committed
547 548
  mdct_butterflies(init,w+n2,n2);
  mdct_bitreverse(init,w);
Mike Whitson's avatar
Mike Whitson committed
549

Monty's avatar
 
Monty committed
550 551 552 553 554 555 556 557 558 559 560
  /* roatate + window */

  T=init->trig+n2;
  x0=out+n2;

  for(i=0;i<n4;i++){
    x0--;
    out[i] =MULT_NORM((w[0]*T[0]+w[1]*T[1])*init->scale);
    x0[0]  =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale);
    w+=2;
    T+=2;
Mike Whitson's avatar
Mike Whitson committed
561 562
  }
}