Coverage Report

Created: 2025-08-26 06:54

/src/vorbis/lib/mdct.c
Line
Count
Source (jump to first uncovered line)
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-2009             *
9
 * by the Xiph.Org Foundation https://xiph.org/                     *
10
 *                                                                  *
11
 ********************************************************************
12
13
 function: normalized modified discrete cosine transform
14
           power of two length transform only [64 <= n ]
15
16
 Original algorithm adapted long ago from _The use of multirate filter
17
 banks for coding of high quality digital audio_, by T. Sporer,
18
 K. Brandenburg and B. Edler, collection of the European Signal
19
 Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp
20
 211-214
21
22
 The below code implements an algorithm that no longer looks much like
23
 that presented in the paper, but the basic structure remains if you
24
 dig deep enough to see it.
25
26
 This module DOES NOT INCLUDE code to generate/apply the window
27
 function.  Everybody has their own weird favorite including me... I
28
 happen to like the properties of y=sin(.5PI*sin^2(x)), but others may
29
 vehemently disagree.
30
31
 ********************************************************************/
32
33
/* this can also be run as an integer transform by uncommenting a
34
   define in mdct.h; the integerization is a first pass and although
35
   it's likely stable for Vorbis, the dynamic range is constrained and
36
   roundoff isn't done (so it's noisy).  Consider it functional, but
37
   only a starting point.  There's no point on a machine with an FPU */
38
39
#include <stdio.h>
40
#include <stdlib.h>
41
#include <string.h>
42
#include <math.h>
43
#include "vorbis/codec.h"
44
#include "mdct.h"
45
#include "os.h"
46
#include "misc.h"
47
48
/* build lookups for trig functions; also pre-figure scaling and
49
   some window function algebra. */
50
51
18.7k
void mdct_init(mdct_lookup *lookup,int n){
52
18.7k
  int   *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4));
53
18.7k
  DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4));
54
55
18.7k
  int i;
56
18.7k
  int n2=n>>1;
57
18.7k
  int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
58
18.7k
  lookup->n=n;
59
18.7k
  lookup->trig=T;
60
18.7k
  lookup->bitrev=bitrev;
61
62
/* trig lookups... */
63
64
6.13M
  for(i=0;i<n/4;i++){
65
6.12M
    T[i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i)));
66
6.12M
    T[i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i)));
67
6.12M
    T[n2+i*2]=FLOAT_CONV(cos((M_PI/(2*n))*(2*i+1)));
68
6.12M
    T[n2+i*2+1]=FLOAT_CONV(sin((M_PI/(2*n))*(2*i+1)));
69
6.12M
  }
70
3.07M
  for(i=0;i<n/8;i++){
71
3.06M
    T[n+i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i+2))*.5);
72
3.06M
    T[n+i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i+2))*.5);
73
3.06M
  }
74
75
  /* bitreverse lookup... */
76
77
18.7k
  {
78
18.7k
    int mask=(1<<(log2n-1))-1,i,j;
79
18.7k
    int msb=1<<(log2n-2);
80
3.07M
    for(i=0;i<n/8;i++){
81
3.06M
      int acc=0;
82
35.1M
      for(j=0;msb>>j;j++)
83
32.0M
        if((msb>>j)&i)acc|=1<<j;
84
3.06M
      bitrev[i*2]=((~acc)&mask)-1;
85
3.06M
      bitrev[i*2+1]=acc;
86
87
3.06M
    }
88
18.7k
  }
89
18.7k
  lookup->scale=FLOAT_CONV(4.f/n);
90
18.7k
}
91
92
/* 8 point butterfly (in place, 4 register) */
93
255M
STIN void mdct_butterfly_8(DATA_TYPE *x){
94
255M
  REG_TYPE r0   = x[6] + x[2];
95
255M
  REG_TYPE r1   = x[6] - x[2];
96
255M
  REG_TYPE r2   = x[4] + x[0];
97
255M
  REG_TYPE r3   = x[4] - x[0];
98
99
255M
           x[6] = r0   + r2;
100
255M
           x[4] = r0   - r2;
101
102
255M
           r0   = x[5] - x[1];
103
255M
           r2   = x[7] - x[3];
104
255M
           x[0] = r1   + r0;
105
255M
           x[2] = r1   - r0;
106
107
255M
           r0   = x[5] + x[1];
108
255M
           r1   = x[7] + x[3];
109
255M
           x[3] = r2   + r3;
110
255M
           x[1] = r2   - r3;
111
255M
           x[7] = r1   + r0;
112
255M
           x[5] = r1   - r0;
113
114
255M
}
115
116
/* 16 point butterfly (in place, 4 register) */
117
127M
STIN void mdct_butterfly_16(DATA_TYPE *x){
118
127M
  REG_TYPE r0     = x[1]  - x[9];
119
127M
  REG_TYPE r1     = x[0]  - x[8];
120
121
127M
           x[8]  += x[0];
122
127M
           x[9]  += x[1];
123
127M
           x[0]   = MULT_NORM((r0   + r1) * cPI2_8);
124
127M
           x[1]   = MULT_NORM((r0   - r1) * cPI2_8);
125
126
127M
           r0     = x[3]  - x[11];
127
127M
           r1     = x[10] - x[2];
128
127M
           x[10] += x[2];
129
127M
           x[11] += x[3];
130
127M
           x[2]   = r0;
131
127M
           x[3]   = r1;
132
133
127M
           r0     = x[12] - x[4];
134
127M
           r1     = x[13] - x[5];
135
127M
           x[12] += x[4];
136
127M
           x[13] += x[5];
137
127M
           x[4]   = MULT_NORM((r0   - r1) * cPI2_8);
138
127M
           x[5]   = MULT_NORM((r0   + r1) * cPI2_8);
139
140
127M
           r0     = x[14] - x[6];
141
127M
           r1     = x[15] - x[7];
142
127M
           x[14] += x[6];
143
127M
           x[15] += x[7];
144
127M
           x[6]  = r0;
145
127M
           x[7]  = r1;
146
147
127M
           mdct_butterfly_8(x);
148
127M
           mdct_butterfly_8(x+8);
149
127M
}
150
151
/* 32 point butterfly (in place, 4 register) */
152
63.9M
STIN void mdct_butterfly_32(DATA_TYPE *x){
153
63.9M
  REG_TYPE r0     = x[30] - x[14];
154
63.9M
  REG_TYPE r1     = x[31] - x[15];
155
156
63.9M
           x[30] +=         x[14];
157
63.9M
           x[31] +=         x[15];
158
63.9M
           x[14]  =         r0;
159
63.9M
           x[15]  =         r1;
160
161
63.9M
           r0     = x[28] - x[12];
162
63.9M
           r1     = x[29] - x[13];
163
63.9M
           x[28] +=         x[12];
164
63.9M
           x[29] +=         x[13];
165
63.9M
           x[12]  = MULT_NORM( r0 * cPI1_8  -  r1 * cPI3_8 );
166
63.9M
           x[13]  = MULT_NORM( r0 * cPI3_8  +  r1 * cPI1_8 );
167
168
63.9M
           r0     = x[26] - x[10];
169
63.9M
           r1     = x[27] - x[11];
170
63.9M
           x[26] +=         x[10];
171
63.9M
           x[27] +=         x[11];
172
63.9M
           x[10]  = MULT_NORM(( r0  - r1 ) * cPI2_8);
173
63.9M
           x[11]  = MULT_NORM(( r0  + r1 ) * cPI2_8);
174
175
63.9M
           r0     = x[24] - x[8];
176
63.9M
           r1     = x[25] - x[9];
177
63.9M
           x[24] += x[8];
178
63.9M
           x[25] += x[9];
179
63.9M
           x[8]   = MULT_NORM( r0 * cPI3_8  -  r1 * cPI1_8 );
180
63.9M
           x[9]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );
181
182
63.9M
           r0     = x[22] - x[6];
183
63.9M
           r1     = x[7]  - x[23];
184
63.9M
           x[22] += x[6];
185
63.9M
           x[23] += x[7];
186
63.9M
           x[6]   = r1;
187
63.9M
           x[7]   = r0;
188
189
63.9M
           r0     = x[4]  - x[20];
190
63.9M
           r1     = x[5]  - x[21];
191
63.9M
           x[20] += x[4];
192
63.9M
           x[21] += x[5];
193
63.9M
           x[4]   = MULT_NORM( r1 * cPI1_8  +  r0 * cPI3_8 );
194
63.9M
           x[5]   = MULT_NORM( r1 * cPI3_8  -  r0 * cPI1_8 );
195
196
63.9M
           r0     = x[2]  - x[18];
197
63.9M
           r1     = x[3]  - x[19];
198
63.9M
           x[18] += x[2];
199
63.9M
           x[19] += x[3];
200
63.9M
           x[2]   = MULT_NORM(( r1  + r0 ) * cPI2_8);
201
63.9M
           x[3]   = MULT_NORM(( r1  - r0 ) * cPI2_8);
202
203
63.9M
           r0     = x[0]  - x[16];
204
63.9M
           r1     = x[1]  - x[17];
205
63.9M
           x[16] += x[0];
206
63.9M
           x[17] += x[1];
207
63.9M
           x[0]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );
208
63.9M
           x[1]   = MULT_NORM( r1 * cPI1_8  -  r0 * cPI3_8 );
209
210
63.9M
           mdct_butterfly_16(x);
211
63.9M
           mdct_butterfly_16(x+16);
212
213
63.9M
}
214
215
/* N point first stage butterfly (in place, 2 register) */
216
STIN void mdct_butterfly_first(DATA_TYPE *T,
217
                                        DATA_TYPE *x,
218
2.43M
                                        int points){
219
220
2.43M
  DATA_TYPE *x1        = x          + points      - 8;
221
2.43M
  DATA_TYPE *x2        = x          + (points>>1) - 8;
222
2.43M
  REG_TYPE   r0;
223
2.43M
  REG_TYPE   r1;
224
225
126M
  do{
226
227
126M
               r0      = x1[6]      -  x2[6];
228
126M
               r1      = x1[7]      -  x2[7];
229
126M
               x1[6]  += x2[6];
230
126M
               x1[7]  += x2[7];
231
126M
               x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
232
126M
               x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
233
234
126M
               r0      = x1[4]      -  x2[4];
235
126M
               r1      = x1[5]      -  x2[5];
236
126M
               x1[4]  += x2[4];
237
126M
               x1[5]  += x2[5];
238
126M
               x2[4]   = MULT_NORM(r1 * T[5]  +  r0 * T[4]);
239
126M
               x2[5]   = MULT_NORM(r1 * T[4]  -  r0 * T[5]);
240
241
126M
               r0      = x1[2]      -  x2[2];
242
126M
               r1      = x1[3]      -  x2[3];
243
126M
               x1[2]  += x2[2];
244
126M
               x1[3]  += x2[3];
245
126M
               x2[2]   = MULT_NORM(r1 * T[9]  +  r0 * T[8]);
246
126M
               x2[3]   = MULT_NORM(r1 * T[8]  -  r0 * T[9]);
247
248
126M
               r0      = x1[0]      -  x2[0];
249
126M
               r1      = x1[1]      -  x2[1];
250
126M
               x1[0]  += x2[0];
251
126M
               x1[1]  += x2[1];
252
126M
               x2[0]   = MULT_NORM(r1 * T[13] +  r0 * T[12]);
253
126M
               x2[1]   = MULT_NORM(r1 * T[12] -  r0 * T[13]);
254
255
126M
    x1-=8;
256
126M
    x2-=8;
257
126M
    T+=16;
258
259
126M
  }while(x2>=x);
260
2.43M
}
261
262
/* N/stage point generic N stage butterfly (in place, 2 register) */
263
STIN void mdct_butterfly_generic(DATA_TYPE *T,
264
                                          DATA_TYPE *x,
265
                                          int points,
266
58.2M
                                          int trigint){
267
268
58.2M
  DATA_TYPE *x1        = x          + points      - 8;
269
58.2M
  DATA_TYPE *x2        = x          + (points>>1) - 8;
270
58.2M
  REG_TYPE   r0;
271
58.2M
  REG_TYPE   r1;
272
273
548M
  do{
274
275
548M
               r0      = x1[6]      -  x2[6];
276
548M
               r1      = x1[7]      -  x2[7];
277
548M
               x1[6]  += x2[6];
278
548M
               x1[7]  += x2[7];
279
548M
               x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
280
548M
               x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
281
282
548M
               T+=trigint;
283
284
548M
               r0      = x1[4]      -  x2[4];
285
548M
               r1      = x1[5]      -  x2[5];
286
548M
               x1[4]  += x2[4];
287
548M
               x1[5]  += x2[5];
288
548M
               x2[4]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
289
548M
               x2[5]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
290
291
548M
               T+=trigint;
292
293
548M
               r0      = x1[2]      -  x2[2];
294
548M
               r1      = x1[3]      -  x2[3];
295
548M
               x1[2]  += x2[2];
296
548M
               x1[3]  += x2[3];
297
548M
               x2[2]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
298
548M
               x2[3]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
299
300
548M
               T+=trigint;
301
302
548M
               r0      = x1[0]      -  x2[0];
303
548M
               r1      = x1[1]      -  x2[1];
304
548M
               x1[0]  += x2[0];
305
548M
               x1[1]  += x2[1];
306
548M
               x2[0]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
307
548M
               x2[1]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
308
309
548M
               T+=trigint;
310
548M
    x1-=8;
311
548M
    x2-=8;
312
313
548M
  }while(x2>=x);
314
58.2M
}
315
316
STIN void mdct_butterflies(mdct_lookup *init,
317
                             DATA_TYPE *x,
318
3.30M
                             int points){
319
320
3.30M
  DATA_TYPE *T=init->trig;
321
3.30M
  int stages=init->log2n-5;
322
3.30M
  int i,j;
323
324
3.30M
  if(--stages>0){
325
2.43M
    mdct_butterfly_first(T,x,points);
326
2.43M
  }
327
328
10.6M
  for(i=1;--stages>0;i++){
329
65.5M
    for(j=0;j<(1<<i);j++)
330
58.2M
      mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
331
7.30M
  }
332
333
67.2M
  for(j=0;j<points;j+=32)
334
63.9M
    mdct_butterfly_32(x+j);
335
336
3.30M
}
337
338
18.7k
void mdct_clear(mdct_lookup *l){
339
18.7k
  if(l){
340
18.7k
    if(l->trig)_ogg_free(l->trig);
341
18.7k
    if(l->bitrev)_ogg_free(l->bitrev);
342
18.7k
    memset(l,0,sizeof(*l));
343
18.7k
  }
344
18.7k
}
345
346
STIN void mdct_bitreverse(mdct_lookup *init,
347
3.30M
                            DATA_TYPE *x){
348
3.30M
  int        n       = init->n;
349
3.30M
  int       *bit     = init->bitrev;
350
3.30M
  DATA_TYPE *w0      = x;
351
3.30M
  DATA_TYPE *w1      = x = w0+(n>>1);
352
3.30M
  DATA_TYPE *T       = init->trig+n;
353
354
255M
  do{
355
255M
    DATA_TYPE *x0    = x+bit[0];
356
255M
    DATA_TYPE *x1    = x+bit[1];
357
358
255M
    REG_TYPE  r0     = x0[1]  - x1[1];
359
255M
    REG_TYPE  r1     = x0[0]  + x1[0];
360
255M
    REG_TYPE  r2     = MULT_NORM(r1     * T[0]   + r0 * T[1]);
361
255M
    REG_TYPE  r3     = MULT_NORM(r1     * T[1]   - r0 * T[0]);
362
363
255M
              w1    -= 4;
364
365
255M
              r0     = HALVE(x0[1] + x1[1]);
366
255M
              r1     = HALVE(x0[0] - x1[0]);
367
368
255M
              w0[0]  = r0     + r2;
369
255M
              w1[2]  = r0     - r2;
370
255M
              w0[1]  = r1     + r3;
371
255M
              w1[3]  = r3     - r1;
372
373
255M
              x0     = x+bit[2];
374
255M
              x1     = x+bit[3];
375
376
255M
              r0     = x0[1]  - x1[1];
377
255M
              r1     = x0[0]  + x1[0];
378
255M
              r2     = MULT_NORM(r1     * T[2]   + r0 * T[3]);
379
255M
              r3     = MULT_NORM(r1     * T[3]   - r0 * T[2]);
380
381
255M
              r0     = HALVE(x0[1] + x1[1]);
382
255M
              r1     = HALVE(x0[0] - x1[0]);
383
384
255M
              w0[2]  = r0     + r2;
385
255M
              w1[0]  = r0     - r2;
386
255M
              w0[3]  = r1     + r3;
387
255M
              w1[1]  = r3     - r1;
388
389
255M
              T     += 4;
390
255M
              bit   += 4;
391
255M
              w0    += 4;
392
393
255M
  }while(w0<w1);
394
3.30M
}
395
396
3.30M
void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
397
3.30M
  int n=init->n;
398
3.30M
  int n2=n>>1;
399
3.30M
  int n4=n>>2;
400
401
  /* rotate */
402
403
3.30M
  DATA_TYPE *iX = in+n2-7;
404
3.30M
  DATA_TYPE *oX = out+n2+n4;
405
3.30M
  DATA_TYPE *T  = init->trig+n4;
406
407
255M
  do{
408
255M
    oX         -= 4;
409
255M
    oX[0]       = MULT_NORM(-iX[2] * T[3] - iX[0]  * T[2]);
410
255M
    oX[1]       = MULT_NORM (iX[0] * T[3] - iX[2]  * T[2]);
411
255M
    oX[2]       = MULT_NORM(-iX[6] * T[1] - iX[4]  * T[0]);
412
255M
    oX[3]       = MULT_NORM (iX[4] * T[1] - iX[6]  * T[0]);
413
255M
    iX         -= 8;
414
255M
    T          += 4;
415
255M
  }while(iX>=in);
416
417
3.30M
  iX            = in+n2-8;
418
3.30M
  oX            = out+n2+n4;
419
3.30M
  T             = init->trig+n4;
420
421
255M
  do{
422
255M
    T          -= 4;
423
255M
    oX[0]       =  MULT_NORM (iX[4] * T[3] + iX[6] * T[2]);
424
255M
    oX[1]       =  MULT_NORM (iX[4] * T[2] - iX[6] * T[3]);
425
255M
    oX[2]       =  MULT_NORM (iX[0] * T[1] + iX[2] * T[0]);
426
255M
    oX[3]       =  MULT_NORM (iX[0] * T[0] - iX[2] * T[1]);
427
255M
    iX         -= 8;
428
255M
    oX         += 4;
429
255M
  }while(iX>=in);
430
431
3.30M
  mdct_butterflies(init,out+n2,n2);
432
3.30M
  mdct_bitreverse(init,out);
433
434
  /* roatate + window */
435
436
3.30M
  {
437
3.30M
    DATA_TYPE *oX1=out+n2+n4;
438
3.30M
    DATA_TYPE *oX2=out+n2+n4;
439
3.30M
    DATA_TYPE *iX =out;
440
3.30M
    T             =init->trig+n2;
441
442
255M
    do{
443
255M
      oX1-=4;
444
445
255M
      oX1[3]  =  MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
446
255M
      oX2[0]  = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
447
448
255M
      oX1[2]  =  MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
449
255M
      oX2[1]  = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
450
451
255M
      oX1[1]  =  MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
452
255M
      oX2[2]  = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
453
454
255M
      oX1[0]  =  MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
455
255M
      oX2[3]  = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
456
457
255M
      oX2+=4;
458
255M
      iX    +=   8;
459
255M
      T     +=   8;
460
255M
    }while(iX<oX1);
461
462
3.30M
    iX=out+n2+n4;
463
3.30M
    oX1=out+n4;
464
3.30M
    oX2=oX1;
465
466
255M
    do{
467
255M
      oX1-=4;
468
255M
      iX-=4;
469
470
255M
      oX2[0] = -(oX1[3] = iX[3]);
471
255M
      oX2[1] = -(oX1[2] = iX[2]);
472
255M
      oX2[2] = -(oX1[1] = iX[1]);
473
255M
      oX2[3] = -(oX1[0] = iX[0]);
474
475
255M
      oX2+=4;
476
255M
    }while(oX2<iX);
477
478
3.30M
    iX=out+n2+n4;
479
3.30M
    oX1=out+n2+n4;
480
3.30M
    oX2=out+n2;
481
255M
    do{
482
255M
      oX1-=4;
483
255M
      oX1[0]= iX[3];
484
255M
      oX1[1]= iX[2];
485
255M
      oX1[2]= iX[1];
486
255M
      oX1[3]= iX[0];
487
255M
      iX+=4;
488
255M
    }while(oX1>oX2);
489
3.30M
  }
490
3.30M
}
491
492
0
void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
493
0
  int n=init->n;
494
0
  int n2=n>>1;
495
0
  int n4=n>>2;
496
0
  int n8=n>>3;
497
0
  DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */
498
0
  DATA_TYPE *w2=w+n2;
499
500
  /* rotate */
501
502
  /* window + rotate + step 1 */
503
504
0
  REG_TYPE r0;
505
0
  REG_TYPE r1;
506
0
  DATA_TYPE *x0=in+n2+n4;
507
0
  DATA_TYPE *x1=x0+1;
508
0
  DATA_TYPE *T=init->trig+n2;
509
510
0
  int i=0;
511
512
0
  for(i=0;i<n8;i+=2){
513
0
    x0 -=4;
514
0
    T-=2;
515
0
    r0= x0[2] + x1[0];
516
0
    r1= x0[0] + x1[2];
517
0
    w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
518
0
    w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
519
0
    x1 +=4;
520
0
  }
521
522
0
  x1=in+1;
523
524
0
  for(;i<n2-n8;i+=2){
525
0
    T-=2;
526
0
    x0 -=4;
527
0
    r0= x0[2] - x1[0];
528
0
    r1= x0[0] - x1[2];
529
0
    w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
530
0
    w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
531
0
    x1 +=4;
532
0
  }
533
534
0
  x0=in+n;
535
536
0
  for(;i<n2;i+=2){
537
0
    T-=2;
538
0
    x0 -=4;
539
0
    r0= -x0[2] - x1[0];
540
0
    r1= -x0[0] - x1[2];
541
0
    w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
542
0
    w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
543
0
    x1 +=4;
544
0
  }
545
546
547
0
  mdct_butterflies(init,w+n2,n2);
548
0
  mdct_bitreverse(init,w);
549
550
  /* roatate + window */
551
552
0
  T=init->trig+n2;
553
0
  x0=out+n2;
554
555
0
  for(i=0;i<n4;i++){
556
0
    x0--;
557
0
    out[i] =MULT_NORM((w[0]*T[0]+w[1]*T[1])*init->scale);
558
0
    x0[0]  =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale);
559
0
    w+=2;
560
0
    T+=2;
561
0
  }
562
0
}