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