several more bugfixes.
[swftools.git] / lib / bladeenc / mdct.c
1 /*
2                         (c) Copyright 1998-2000 - Tord Jansson
3                         ======================================
4
5                 This file is part of the BladeEnc MP3 Encoder, based on
6                 ISO's reference code for MPEG Layer 3 compression, and might
7                 contain smaller or larger sections that are directly taken
8                 from ISO's reference code.
9
10                 All changes to the ISO reference code herein are either
11                 copyrighted by Tord Jansson (tord.jansson@swipnet.se)
12                 or sublicensed to Tord Jansson by a third party.
13
14         BladeEnc is free software; you can redistribute this file
15         and/or modify it under the terms of the GNU Lesser General Public
16         License as published by the Free Software Foundation; either
17         version 2.1 of the License, or (at your option) any later version.
18
19
20
21         ------------    Changes    ------------
22
23         2000-11-04  Andre Piotrowski
24
25         -       some cosmetics
26         -       included "l3psy.h" to use the defined block type constants
27         -       negligible speed up of mdct_sub()
28         -       speed up of mdct()
29
30         2000-12-10  ap
31
32         -       speed up of mdct_sub() and mdct() - MDCT_CHANGE_LEVEL 5
33 */
34
35 /**********************************************************************
36  * ISO MPEG Audio Subgroup Software Simulation Group (1996)
37  * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension
38  *
39  * $Id: mdct.c,v 1.1 2002/01/10 17:30:01 kramm Exp $
40  *
41  * $Log: mdct.c,v $
42  * Revision 1.1  2002/01/10 17:30:01  kramm
43  * Version 0.94.1 of the bladeenc mp3 encoder
44  *
45  * Revision 1.1  1996/02/14 04:04:23  rowlands
46  * Initial revision
47  *
48  * Received from Mike Coleman
49  **********************************************************************/
50
51 #include "common.h"
52 #include <assert.h>
53 #include "l3side.h"
54 #include "l3psy.h"
55 #include "mdct.h"
56
57
58
59
60
61 /*
62         This is table B.9: coefficients for aliasing reduction
63 */
64 static  double                  c[8] = { -0.6, -0.535, -0.33, -0.185, -0.095, -0.041, -0.0142, -0.0037 };
65
66
67
68
69
70 int                                             gr_idx[3] = {0,1,2};
71
72
73
74
75
76 #if 0
77 static  double                  ca[8];
78 #endif
79 static  double                  cs[8];
80
81
82
83
84
85 int     fInit_mdct_sub;
86
87 void mdct_sub
88 (
89         L3SBS                                   (*sb_sample),
90         double                                  (*mdct_freq)[2][576],
91         int                                             stereo,
92         III_side_info_t                 *l3_side,
93         int                                             mode_gr
94 )
95 {
96         gr_info                                 *cod_info;
97 #if MDCT_CHANGE_LEVEL < 5
98         double                                  mdct_in[36];
99 #endif
100         int                                             gr, ch, band, k, j;
101         double                                  bu, bd;
102         int                                             block_type;
103         double                                  (*mdct_enc_long)[2][32][18] = (double (*)[2][32][18]) mdct_freq;
104
105         if (!fInit_mdct_sub)
106         {
107                 /* prepare the aliasing reduction butterflies */
108                 for (k = 0;  k < 8;  k++)
109                 {
110                         double sq = sqrt (1.0 + c[k]*c[k]);
111 #if 0
112                         ca[k] = c[k] / sq;
113 #endif
114                         cs[k] = 1.0 / sq;
115                 }
116                 fInit_mdct_sub++;
117         }
118
119
120         for (gr = 0;  gr < mode_gr;  gr++)
121         {
122                 int             pre_gr = gr_idx[gr  ];
123                 int             cur_gr = gr_idx[gr+1];
124
125                 for (ch = 0;  ch < stereo;  ch++)
126                 {
127                         double (*mdct_enc)[18] = mdct_enc_long[gr][ch];
128
129                         cod_info = &l3_side->gr[gr].ch[ch].tt;
130                         block_type = cod_info->block_type;
131
132                         /* Compensate for inversion in the analysis filter */
133                         for (k = 1;  k < 18;  k++)
134                                 if (k & 1)
135                                         for (band = 1;  band < 32;  band++)
136                                                 if (band & 1)
137                                                         (*sb_sample)[ch][cur_gr][k][band] *= -1.0;
138                         /*
139                                 Perform imdct of 18 previous subband samples
140                                 + 18 current subband samples
141                         */
142                         for (band = 0;  band < 32;  band++)
143                         {
144 #if MDCT_CHANGE_LEVEL < 5
145                                 for (k = 0;  k < 18;  k++)
146                                 {
147                                         mdct_in[k]    = (*sb_sample)[ch][pre_gr][k][band];
148                                         mdct_in[k+18] = (*sb_sample)[ch][cur_gr][k][band];
149                                 }
150                                 if (cod_info->mixed_block_flag  &&  (band < 2))
151                                         block_type = NORM_TYPE;   /* AND WHEN WILL THE BLOCK_TYPE BE SWITCHED BACK? */
152
153                                         /* &mdct_enc[gr][ch][band][0] */
154                                 mdct (mdct_in, mdct_enc/*[gr][ch]*/[band], block_type);
155 #else
156                                 mdct ((*sb_sample)[ch][pre_gr], (*sb_sample)[ch][cur_gr], band, mdct_enc[band], block_type);
157 #endif
158                         }
159
160                         /*
161                                 Perform aliasing reduction butterfly
162                                 on long blocks
163                         */
164
165                         if (block_type != SHORT_TYPE)
166                                 for (band = 0;  band < 31;  band++)
167                                         for (k = 0;  k < 8;  k++)
168                                         {
169 #if 1   /* This is faster because the calculation can be done more sequential */
170                                                 bu = (mdct_enc[band  ][17-k] + mdct_enc[band+1][   k] * c[k]) * cs[k];
171                                                 bd = (mdct_enc[band+1][   k] - mdct_enc[band  ][17-k] * c[k]) * cs[k];
172                                                 mdct_enc[band  ][17-k] = bu;
173                                                 mdct_enc[band+1][   k] = bd;
174 #else
175                                                 bu = mdct_enc[band  ][17-k] * cs[k] + mdct_enc[band+1][   k] * ca[k];
176                                                 bd = mdct_enc[band+1][   k] * cs[k] - mdct_enc[band  ][17-k] * ca[k];
177                                                 mdct_enc[band  ][17-k] = bu;
178                                                 mdct_enc[band+1][   k] = bd;
179 #endif
180                                         }
181                 }
182         }
183
184         /*
185                 Save latest granule's subband samples to be used in
186                 the next mdct call
187         */
188         j = gr_idx[mode_gr];
189         for (k = mode_gr;  k > 0;  k--)
190                 gr_idx[k] = gr_idx[k-1];
191         gr_idx[0] = j;
192 }
193
194
195
196
197
198 /*-------------------------------------------------------------------*/
199 /*                                                                   */
200 /*   Function: Calculation of the MDCT                               */
201 /*   In the case of long blocks ( block_type 0,1,3 ) there are       */
202 /*   36 coefficents in the time domain and 18 in the frequency       */
203 /*   domain.                                                         */
204 /*   In the case of short blocks (block_type 2 ) there are 3         */
205 /*   transformations with short length. This leads to 12 coefficents */
206 /*   in the time and 6 in the frequency domain. In this case the     */
207 /*   results are stored side by side in the vector out[].            */
208 /*                                                                   */
209 /*   New layer3                                                      */
210 /*                                                                   */
211 /*-------------------------------------------------------------------*/
212
213 int     fInit_mdct;
214
215
216
217 #if MDCT_CHANGE_LEVEL == 5
218
219 void mdct
220 (
221         double                                  inA[18][32],
222         double                                  inB[18][32],
223         int                                             band,
224         double                                  *out,
225         int                                             block_type
226 )
227 {
228         static  double                  cos_l[18][18];
229         static  double                  cos_s[ 6][ 6];
230
231         static  double                  winNorm [18];
232         static  double                  winShort[ 6];
233
234         int                                             k, m, N;
235         double                                  t1, t2, t3, t4, t5, t6, u1, u2;
236
237
238         if (!fInit_mdct)
239         {
240                 N = 36;
241                 for (k = 0;  k < N/2;  k++)  winNorm [k] = sin (PI/36 * (k+0.5));
242                 for (k = 0;  k < N/4;  k++)  for (m = 0;  m < N/2;  m++)  cos_l[k][m] = cos((PI/(2*N)) * (2* k     +1+N/2) * (2*m+1)) / (N/4);
243                 for (     ;  k < N/2;  k++)  for (m = 0;  m < N/2;  m++)  cos_l[k][m] = cos((PI/(2*N)) * (2*(k+N/4)+1+N/2) * (2*m+1)) / (N/4);
244
245                 N = 12;
246                 for (k = 0;  k < N/2;  k++)  winShort[k] = sin (PI/12 * (k+0.5));
247                 for (k = 0;  k < N/4;  k++)  for (m = 0;  m < N/2;  m++)  cos_s[k][m] = cos((PI/(2*N)) * (2*k      +1+N/2) * (2*m+1)) / (N/4);
248                 for (     ;  k < N/2;  k++)  for (m = 0;  m < N/2;  m++)  cos_s[k][m] = cos((PI/(2*N)) * (2*(k+N/4)+1+N/2) * (2*m+1)) / (N/4);
249
250                 fInit_mdct++;
251         }
252
253         for (m = 0;  m < 18;  m++)  out[m] = 0.0;
254         switch (block_type)
255         {
256                 case  NORM_TYPE:
257                         for (k = 0;  k < 9;  k++)
258                         {
259                                 t1 = winNorm [   k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
260                                 t2 = winNorm [17-k] * inB[k][band] + winNorm [   k] * inB[17-k][band];
261                                 for (m = 0;  m < 18;  m++)
262                                         out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
263                         }
264                         break;
265
266                 case START_TYPE:
267                         for (k = 0;  k < 6;  k++)
268                         {
269                                 t1 = winNorm [   k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
270                                 t2 =                  inB[k][band];
271                                 for (m = 0;  m < 18;  m++)
272                                         out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
273                         }
274                         for (k =  6;  k <  9;  k++)
275                         {
276                                 t1 = winNorm [   k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
277                                 t2 = winShort[11-k] * inB[k][band] + winShort[k- 6] * inB[17-k][band];
278                                 for (m = 0;  m < 18;  m++)
279                                         out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
280                         }
281                         break;
282
283                 case  STOP_TYPE:
284                         for (k =  0;  k <  6;  k++)
285                         {
286                                 t1 =                                                 -inA[17-k][band];
287                                 t2 = winNorm [17-k] * inB[k][band] + winNorm [   k] * inB[17-k][band];
288                                 for (m = 0;  m < 18;  m++)
289                                         out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
290                         }
291                         for (k =  6;  k <  9;  k++)
292                         {
293                                 t1 = winShort[k- 6] * inA[k][band] - winShort[11-k] * inA[17-k][band];
294                                 t2 = winNorm [17-k] * inB[k][band] + winNorm [   k] * inB[17-k][band];
295                                 for (m = 0;  m < 18;  m++)
296                                         out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
297                         }
298                         break;
299
300                 case SHORT_TYPE:
301                         for (k =  0;  k <  3;  k++)
302                         {
303                                 u1 = winShort[k];  u2 = winShort[5-k];
304                                 t1 = u1 * inA[k+ 6][band] - u2 * inA[11-k][band];
305                                 t2 = u2 * inA[k+12][band] + u1 * inA[17-k][band];
306                                 t3 = u1 * inA[k+12][band] - u2 * inA[17-k][band];
307                                 t4 = u2 * inB[k   ][band] + u1 * inB[ 5-k][band];
308                                 t5 = u1 * inB[k   ][band] - u2 * inB[ 5-k][band];
309                                 t6 = u2 * inB[k+ 6][band] + u1 * inB[11-k][band];
310                                 for (m = 0;  m <  6;  m++)
311                                 {
312                                         u1 = cos_s[k][m];  u2 = cos_s[k+3][m];
313                                         out[3*m  ] += t1 * u1 + t2 * u2;
314                                         out[3*m+1] += t3 * u1 + t4 * u2;
315                                         out[3*m+2] += t5 * u1 + t6 * u2;
316                                 }
317                         }
318         }
319 }
320
321 #elif MDCT_CHANGE_LEVEL == 4   /* reduce number of multiplications to nearly the half once more! (cos_x[9+k][m] = -cos_x[9-k][m] and cos_x[35+k][m] = cos_x[35-k][m]) */
322
323 void mdct
324 (
325         double                                  *in,
326         double                                  *out,
327         int                                             block_type
328 )
329 {
330         static  double                  cos_l[36][18];
331         static  double                  cos_s[12][ 6];
332
333         static  double                  winNorm [36];
334         static  double                  winShort[12];
335         static  double                  *winStart = winShort-18;
336         static  double                  *winStop  = winShort- 6;
337
338         int                                             l, k, m, N;
339         double                                  temp;
340
341
342         if (!fInit_mdct)
343         {
344                 N = 36;
345                 for (k = 0;  k < N;  k++)  winNorm [k] = sin (PI/36 * (k+0.5));
346                 for (k = 0;  k < N;  k++)  for (m = 0;  m < N/2;  m++)  cos_l[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
347
348                 N = 12;
349                 for (k = 0;  k < N;  k++)  winShort[k] = sin (PI/12 * (k+0.5));
350                 for (k = 0;  k < N;  k++)  for (m = 0;  m < N/2;  m++)  cos_s[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
351
352                 fInit_mdct++;
353         }
354
355         for (m = 0;  m < 18;  m++)  out[m] = 0.0;
356         switch (block_type)
357         {
358                 case  NORM_TYPE:
359                         for (k =  0;  k <  9;  k++)  {temp = (winNorm [k] * in[k] - winNorm [17-k] * in[17-k]);  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
360                         for (k = 18;  k < 27;  k++)  {temp = (winNorm [k] * in[k] + winNorm [53-k] * in[53-k]);  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
361                         break;
362
363                 case START_TYPE:
364                         for (k =  0;  k <  9;  k++)  {temp = (winNorm [k] * in[k] - winNorm [17-k] * in[17-k]);  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
365                         for (k = 18;  k < 24;  k++)  {temp =               in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
366                         for (k = 24;  k < 27;  k++)  {temp = (winStart[k] * in[k] + winStart[53-k] * in[53-k]);  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
367                         break;
368
369                 case  STOP_TYPE:
370                         for (k =  6;  k <  9;  k++)  {temp = (winStop [k] * in[k] - winStop [17-k] * in[17-k]);  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
371                         for (k = 12;  k < 18;  k++)  {temp =               in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
372                         for (k = 18;  k < 27;  k++)  {temp = (winNorm [k] * in[k] + winNorm [53-k] * in[53-k]);  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
373                         break;
374
375                 case SHORT_TYPE:
376                         for (l = 0;  l < 3;  l++)
377                         {
378                                 for (k =  0;  k <  3;  k++)  {temp = (winShort[k] * in[k+6+l*6] - winShort[ 5-k] * in[ 5-k+6+l*6]);  for (m = 0;  m <  6;  m++)  out[3*m+l] += temp * cos_s[k][m];}
379                                 for (k =  6;  k <  9;  k++)  {temp = (winShort[k] * in[k+6+l*6] + winShort[17-k] * in[17-k+6+l*6]);  for (m = 0;  m <  6;  m++)  out[3*m+l] += temp * cos_s[k][m];}
380                         }
381         }
382 }
383
384 #elif MDCT_CHANGE_LEVEL == 3   /* reduce number of multiplications to nearly the half (win_x[k]*in[k] is constant in m-loop)! flip cos_x components for faster access */
385
386 void mdct
387 (
388         double                                  *in,
389         double                                  *out,
390         int                                             block_type
391 )
392 {
393         static  double                  cos_l[36][18];
394         static  double                  cos_s[12][ 6];
395
396         static  double                  winNorm [36];
397         static  double                  winShort[12];
398         static  double                  *winStart = winShort-18;
399         static  double                  *winStop  = winShort- 6;
400
401         int                                             l, k, m, N;
402         double                                  temp;
403
404
405         if (!fInit_mdct)
406         {
407                 N = 36;
408                 for (k = 0;  k < N;  k++)  winNorm [k] = sin (PI/36 * (k+0.5));
409                 for (k = 0;  k < N;  k++)  for (m = 0;  m < N/2;  m++)  cos_l[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
410
411                 N = 12;
412                 for (k = 0;  k < N;  k++)  winShort[k] = sin (PI/12 * (k+0.5));
413                 for (k = 0;  k < N;  k++)  for (m = 0;  m < N/2;  m++)  cos_s[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
414
415                 fInit_mdct++;
416         }
417
418         for (m = 0;  m < 18;  m++)  out[m] = 0.0;
419         switch (block_type)
420         {
421                 case  NORM_TYPE:
422                         for (k =  0;  k < 36;  k++)  {temp = winNorm [k] * in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
423                         break;
424
425                 case START_TYPE:
426                         for (k =  0;  k < 18;  k++)  {temp = winNorm [k] * in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
427                         for (k = 18;  k < 24;  k++)  {temp =               in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
428                         for (k = 24;  k < 30;  k++)  {temp = winStart[k] * in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
429                         break;
430
431                 case  STOP_TYPE:
432                         for (k =  6;  k < 12;  k++)  {temp = winStop [k] * in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
433                         for (k = 12;  k < 18;  k++)  {temp =               in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
434                         for (k = 18;  k < 36;  k++)  {temp = winNorm [k] * in[k];  for (m = 0;  m < 18;  m++)  out[m] += temp * cos_l[k][m];}
435                         break;
436
437                 case SHORT_TYPE:
438                         for (l = 0;  l < 3;  l++)
439                                 for (k =  0;  k < 12;  k++)  {temp = winShort[k] * in[k+6+l*6];  for (m = 0;  m <  6;  m++)  out[3*m+l] += temp * cos_s[k][m];}
440         }
441 }
442
443 #elif MDCT_CHANGE_LEVEL == 2   /* avoid calculating of some redundant values; take care for some special values */
444
445 void mdct
446 (
447         double                                  *in,
448         double                                  *out,
449         int                                             block_type
450 )
451 {
452         static  double                  cos_l[18][36];
453         static  double                  cos_s[ 6][12];
454
455         static  double                  winNorm [36];
456         static  double                  winShort[12];
457         static  double                  *winStart = winShort-18;
458         static  double                  *winStop  = winShort- 6;
459
460         int                                             l, k, m, N;
461         double                                  sum;
462
463
464         if (!fInit_mdct)
465         {
466                 for (k =  0;  k < 36;  k++)  winNorm [k] = sin (PI/36 * (k+0.5));
467                 for (k =  0;  k < 12;  k++)  winShort[k] = sin (PI/12 * (k+0.5));
468
469                 N = 36;  for (m = 0;  m < N/2;  m++)  for (k = 0;  k < N;  k++)  cos_l[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
470                 N = 12;  for (m = 0;  m < N/2;  m++)  for (k = 0;  k < N;  k++)  cos_s[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
471
472                 fInit_mdct++;
473         }
474
475         switch (block_type)
476         {
477                 case  NORM_TYPE:
478                         for (m = 0;  m < 18;  m++)
479                         {
480                                 sum = 0.0;
481                                 for (k =  0;  k < 36;  k++) sum += winNorm [k] * in[k] * cos_l[m][k];
482                                 out[m] = sum;
483                         }
484                         break;
485
486                 case START_TYPE:
487                         for (m = 0;  m < 18;  m++)
488                         {
489                                 sum = 0.0;
490                                 for (k =  0;  k < 18;  k++) sum += winNorm [k] * in[k] * cos_l[m][k];
491                                 for (k = 18;  k < 24;  k++) sum +=               in[k] * cos_l[m][k];
492                                 for (k = 24;  k < 30;  k++) sum += winStart[k] * in[k] * cos_l[m][k];
493                                 out[m] = sum;
494                         }
495                         break;
496
497                 case  STOP_TYPE:
498                         for (m = 0;  m < 18;  m++)
499                         {
500                                 sum = 0.0;
501                                 for (k =  6;  k < 12;  k++) sum += winStop [k] * in[k] * cos_l[m][k];
502                                 for (k = 12;  k < 18;  k++) sum +=               in[k] * cos_l[m][k];
503                                 for (k = 18;  k < 36;  k++) sum += winNorm [k] * in[k] * cos_l[m][k];
504                                 out[m] = sum;
505                         }
506                         break;
507
508                 case SHORT_TYPE:
509                         for (l = 0;  l < 3;  l++)
510                         {
511                                 for (m = 0;  m <  6;  m++)
512                                 {
513                                         sum = 0.0;
514                                         for (k = 0;  k < 12;  k++) sum += winShort[k] * in[k+6+l*6] * cos_s[m][k];
515                                         out[3*m+l] = sum;
516                                 }
517                         }
518         }
519 }
520
521 #elif MDCT_CHANGE_LEVEL == 1   /* reformatted for better overview; use block type constants */
522
523 void mdct
524 (
525         double                                  *in,
526         double                                  *out,
527         int                                             block_type
528 )
529 {
530         static  double                  cos_l[18][36];
531         static  double                  cos_s[ 6][12];
532
533         static  double                  winNorm [36];
534         static  double                  winShort[12];
535         static  double                  winStart[36];
536         static  double                  winStop [36];
537
538         int                                             l, k, m, N;
539         double                                  sum;
540
541
542         if (!fInit_mdct)
543         {
544                 /* type 0 -- NORM_TYPE */
545                 for (k =  0;  k < 36;  k++)  winNorm [k] = sin (PI/36 * (k+0.5));
546                 /* type 1 -- START_TYPE */
547                 for (k =  0;  k < 18;  k++)  winStart[k] = sin (PI/36 * (k+0.5));
548                 for (k = 18;  k < 24;  k++)  winStart[k] = 1.0;
549                 for (k = 24;  k < 30;  k++)  winStart[k] = sin (PI/12 * (k+0.5 - 18));
550                 for (k = 30;  k < 36;  k++)  winStart[k] = 0.0;
551                 /* type 3 -- STOP_TYPE */
552                 for (k =  0;  k <  6;  k++)  winStop [k] = 0.0;
553                 for (k =  6;  k < 12;  k++)  winStop [k] = sin (PI/12 * (k+0.5 -  6));
554                 for (k = 12;  k < 18;  k++)  winStop [k] = 1.0;
555                 for (k = 18;  k < 36;  k++)  winStop [k] = sin (PI/36 * (k+0.5));
556                 /* type 2 -- SHORT_TYPE */
557                 for (k =  0;  k < 12;  k++)  winShort[k] = sin (PI/12 * (k+0.5));
558
559                 N = 12;  for (m = 0;  m < N/2;  m++)  for (k = 0;  k < N;  k++)  cos_s[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
560                 N = 36;  for (m = 0;  m < N/2;  m++)  for (k = 0;  k < N;  k++)  cos_l[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4);
561
562                 fInit_mdct++;
563         }
564
565         switch (block_type)
566         {
567         case  NORM_TYPE:  N = 36;                             for (m = 0;  m < N/2;  m++)  {sum = 0.0;  for (k = 0;  k < N;  k++) sum += winNorm [k] * in[k      ] * cos_l[m][k];  out[  m  ] = sum;}  break;
568         case START_TYPE:  N = 36;                             for (m = 0;  m < N/2;  m++)  {sum = 0.0;  for (k = 0;  k < N;  k++) sum += winStart[k] * in[k      ] * cos_l[m][k];  out[  m  ] = sum;}  break;
569         case  STOP_TYPE:  N = 36;                             for (m = 0;  m < N/2;  m++)  {sum = 0.0;  for (k = 0;  k < N;  k++) sum += winStop [k] * in[k      ] * cos_l[m][k];  out[  m  ] = sum;}  break;
570         case SHORT_TYPE:  N = 12;  for (l = 0;  l < 3;  l++) {for (m = 0;  m < N/2;  m++)  {sum = 0.0;  for (k = 0;  k < N;  k++) sum += winShort[k] * in[k+6+l*6] * cos_s[m][k];  out[3*m+l] = sum;}}
571         }
572 }
573
574 #endif /* MDCT_CHANGE_LEVEL */