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 */