home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 13 / AACD13.ISO / AACD / Sound / LAME / src / newmdct.c < prev    next >
C/C++ Source or Header  |  2000-07-26  |  21KB  |  722 lines
  1. /*
  2.  *    MP3 window subband -> subband filtering -> mdct routine
  3.  *
  4.  *    Copyright (c) 1999 Takehiro TOMINAGA
  5.  *
  6.  *
  7.  * This library is free software; you can redistribute it and/or
  8.  * modify it under the terms of the GNU Library General Public
  9.  * License as published by the Free Software Foundation; either
  10.  * version 2 of the License, or (at your option) any later version.
  11.  *
  12.  * This library is distributed in the hope that it will be useful,
  13.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.     See the GNU
  15.  * Library General Public License for more details.
  16.  *
  17.  * You should have received a copy of the GNU Library General Public
  18.  * License along with this library; if not, write to the
  19.  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  20.  * Boston, MA 02111-1307, USA.
  21.  */
  22.  
  23. /*
  24.  *         Special Thanks to Patrick De Smet for your advices.
  25.  */
  26.  
  27.  
  28. #include "util.h"
  29. #include "l3side.h"
  30. #include "newmdct.h"
  31.  
  32. #define SCALE 32768
  33.  
  34. #ifndef USE_GOGO_SUBBAND
  35. static const FLOAT8 enwindow[] = 
  36. {
  37.   3.5758972e-02, 3.401756e-03,  9.83715e-04,   9.9182e-05, /* 15*/
  38.       -4.77e-07,  1.03951e-04,  9.53674e-04, 2.841473e-03,
  39.      1.2398e-05,  1.91212e-04, 2.283096e-03,1.6994476e-02,
  40.   1.8756866e-02, 2.630711e-03,  2.47478e-04,   1.4782e-05,
  41.  
  42.   3.5694122e-02, 3.643036e-03,  9.91821e-04,   9.6321e-05, /* 14*/
  43.       -4.77e-07,  1.05858e-04,  9.30786e-04, 2.521515e-03,
  44.      1.1444e-05,  1.65462e-04, 2.110004e-03,1.6112804e-02,
  45.   1.9634247e-02, 2.803326e-03,  2.77042e-04,   1.6689e-05,
  46.  
  47.   -3.5586357e-02, -3.858566e-03,  -9.95159e-04,   -9.3460e-05, /* 13*/
  48.         4.77e-07,  -1.07288e-04,  -9.02653e-04, -2.174854e-03,
  49.      -1.0014e-05,  -1.40190e-04, -1.937389e-03,-1.5233517e-02,
  50.   -2.0506859e-02, -2.974033e-03,  -3.07560e-04,   -1.8120e-05,
  51.  
  52.   3.5435200e-02, 4.049301e-03,  9.94205e-04,   9.0599e-05, /* 12*/
  53.       -4.77e-07,  1.08242e-04,  8.68797e-04, 1.800537e-03,
  54.       9.060e-06,  1.16348e-04, 1.766682e-03,1.4358521e-02,
  55.   2.1372318e-02,  3.14188e-03,  3.39031e-04,   1.9550e-05,
  56.  
  57.   -3.5242081e-02, -4.215240e-03, - 9.89437e-04,-   8.7261e-05, /* 11*/
  58.         4.77e-07, - 1.08719e-04, - 8.29220e-04,- 1.399517e-03,
  59.   -    8.106e-06, -  9.3937e-05, -1.597881e-03,-1.3489246e-02,
  60.   -2.2228718e-02, -3.306866e-03, - 3.71456e-04,-   2.1458e-05,
  61.  
  62.   3.5007000e-02, 4.357815e-03,  9.80854e-04,   8.3923e-05, /* 10*/
  63.       -4.77e-07,  1.08719e-04,   7.8392e-04,  9.71317e-04,
  64.       7.629e-06,   7.2956e-05, 1.432419e-03,1.2627602e-02,
  65.   2.3074150e-02, 3.467083e-03,  4.04358e-04,   2.3365e-05,
  66.  
  67.   3.4730434e-02, 4.477024e-03,  9.68933e-04,   8.0585e-05, /* 9*/
  68.       -9.54e-07,  1.08242e-04,  7.31945e-04,  5.15938e-04,
  69.       6.676e-06,   5.2929e-05, 1.269817e-03,1.1775017e-02,
  70.   2.3907185e-02, 3.622532e-03,  4.38213e-04,   2.5272e-05,
  71.  
  72.   3.4412861e-02, 4.573822e-03,  9.54151e-04,   7.6771e-05, /* 8*/
  73.       -9.54e-07,  1.06812e-04,  6.74248e-04,   3.3379e-05,
  74.       6.199e-06,   3.4332e-05, 1.111031e-03,1.0933399e-02,
  75.   2.4725437e-02, 3.771782e-03,  4.72546e-04,   2.7657e-05,
  76.  
  77.   3.4055710e-02, 4.649162e-03,  9.35555e-04,   7.3433e-05, /* 7*/
  78.       -9.54e-07,  1.05381e-04,  6.10352e-04, -4.75883e-04,
  79.       5.245e-06,   1.7166e-05,  9.56535e-04,1.0103703e-02,
  80.   2.5527000e-02, 3.914356e-03,  5.07355e-04,   3.0041e-05,
  81.  
  82.   -3.3659935e-02, -4.703045e-03,  -9.15051e-04,-   7.0095e-05, /* 6*/
  83.         9.54e-07, - 1.02520e-04,  -5.39303e-04,  1.011848e-03,
  84.   -    4.768e-06, -    9.54e-07,  -8.06808e-04,- 9.287834e-03,
  85.   -2.6310921e-02, -4.048824e-03,  -5.42164e-04,-   3.2425e-05,
  86.  
  87.   3.3225536e-02, 4.737377e-03,  8.91685e-04,   6.6280e-05,  /* 5*/
  88.      -1.431e-06,   9.9182e-05,  4.62532e-04,-1.573563e-03,
  89.       4.292e-06,  -1.3828e-05,  6.61850e-04, 8.487225e-03,
  90.   2.7073860e-02, 4.174709e-03,  5.76973e-04,   3.4809e-05,
  91.  
  92.   3.2754898e-02, 4.752159e-03,  8.66413e-04,   6.2943e-05, /* 4*/
  93.      -1.431e-06,   9.5367e-05,  3.78609e-04,-2.161503e-03,
  94.       3.815e-06,   -2.718e-05,  5.22137e-04, 7.703304e-03,
  95.   2.7815342e-02, 4.290581e-03,  6.11782e-04,   3.7670e-05,
  96.  
  97.   -3.2248020e-02, -4.748821e-03,  -8.38757e-04,-   5.9605e-05, /* 3*/
  98.        1.907e-06, -  9.0122e-05,  -2.88486e-04,  2.774239e-03,
  99.   -    3.338e-06,    3.9577e-05,  -3.88145e-04,- 6.937027e-03,
  100.   -2.8532982e-02, -4.395962e-03,  -6.46591e-04,-   4.0531e-05,
  101.  
  102.   3.1706810e-02, 4.728317e-03,  8.09669e-04,    5.579e-05,
  103.      -1.907e-06,   8.4400e-05,  1.91689e-04,-3.411293e-03,
  104.       3.338e-06,  -5.0545e-05,  2.59876e-04, 6.189346e-03,
  105.   2.9224873e-02, 4.489899e-03,  6.80923e-04,   4.3392e-05,
  106.  
  107.   3.1132698e-02, 4.691124e-03,  7.79152e-04,   5.2929e-05,
  108.      -2.384e-06,   7.7724e-05,   8.8215e-05,-4.072189e-03,
  109.       2.861e-06,  -6.0558e-05,  1.37329e-04, 5.462170e-03,
  110.   2.9890060e-02, 4.570484e-03,  7.14302e-04,   4.6253e-05,
  111.  
  112.   3.5780907e-02,1.7876148e-02, 3.134727e-03, 2.457142e-03,
  113.     9.71317e-04,  2.18868e-04,  1.01566e-04,   1.3828e-05,
  114.  
  115.   3.0526638e-02, 4.638195e-03,  7.47204e-04,   4.9591e-05,
  116.    4.756451e-03,   2.1458e-05,  -6.9618e-05,    2.384e-06,
  117. };
  118. #endif
  119.  
  120.  
  121. #define NS 12
  122. #define NL 36
  123.  
  124. static FLOAT8 ca[8], cs[8];
  125. static FLOAT8 win[4][36];
  126. static FLOAT8 tantab_l[NL/4];
  127. static FLOAT8 tritab_s[NS/4*2];
  128.  
  129. /************************************************************************
  130. *
  131. * window_subband()
  132. *
  133. * PURPOSE:  Overlapping window on PCM samples
  134. *
  135. * SEMANTICS:
  136. * 32 16-bit pcm samples are scaled to fractional 2's complement and
  137. * concatenated to the end of the window buffer #x#. The updated window
  138. * buffer #x# is then windowed by the analysis window #c# to produce the
  139. * windowed sample #z#
  140. *
  141. ************************************************************************/
  142.  
  143. /*
  144.  *      new IDCT routine written by Naoki Shibata
  145.  */
  146. INLINE static void idct32(FLOAT8 a[])
  147. {
  148.     /* returns sum_j=0^31 a[j]*cos(PI*j*(k+1/2)/32), 0<=k<32 */
  149.  
  150.     static const FLOAT8 costab[] = {
  151.     0.54119610014619701222,
  152.     1.3065629648763763537,
  153.  
  154.     0.50979557910415917998,
  155.     2.5629154477415054814,
  156.     0.89997622313641556513,
  157.     0.60134488693504528634,
  158.  
  159.     0.5024192861881556782,
  160.     5.1011486186891552563,
  161.     0.78815462345125020249,
  162.     0.64682178335999007679,
  163.     0.56694403481635768927,
  164.     1.0606776859903470633,
  165.     1.7224470982383341955,
  166.     0.52249861493968885462,
  167.  
  168.     10.19000812354803287,
  169.     0.674808341455005678,
  170.     1.1694399334328846596,
  171.     0.53104259108978413284,    
  172.     2.0577810099534108446,
  173.     0.58293496820613388554,
  174.     0.83934964541552681272,
  175.     0.50547095989754364798,
  176.     3.4076084184687189804,
  177.     0.62250412303566482475,
  178.     0.97256823786196078263,
  179.     0.51544730992262455249,
  180.     1.4841646163141661852,
  181.     0.5531038960344445421,
  182.     0.74453627100229857749,
  183.     0.5006029982351962726,
  184.     };
  185.  
  186.     int i;
  187. #ifdef USE_GOGO_SUBBAND
  188.     a[31] += a[29]; a[29] += a[27];
  189.     a[27] += a[25]; a[25] += a[23];
  190.     a[23] += a[21]; a[21] += a[19];
  191.     a[19] += a[17]; a[17] += a[15];
  192.     a[15] += a[13]; a[13] += a[11];
  193.     a[11] += a[ 9]; a[ 9] += a[ 7];
  194.     a[ 7] += a[ 5]; a[ 5] += a[ 3];
  195.     a[ 3] += a[ 1];
  196.  
  197.     a[31] += a[27]; a[27] += a[23];
  198.     a[30] += a[26]; a[26] += a[22];
  199.     a[23] += a[19]; a[19] += a[15];
  200.     a[22] += a[18]; a[18] += a[14];
  201.     a[15] += a[11]; a[11] += a[ 7];
  202.     a[14] += a[10]; a[10] += a[ 6];
  203.     a[ 7] += a[ 3];
  204.     a[ 6] += a[ 2];
  205.  
  206.     a[31] += a[23]; a[23] += a[15];    a[15] += a[ 7];
  207.     a[30] += a[22]; a[22] += a[14];    a[14] += a[ 6];
  208.     a[29] += a[21]; a[21] += a[13];    a[13] += a[ 5];
  209.     a[28] += a[20]; a[20] += a[12];    a[12] += a[ 4];
  210.  
  211.     a[ 6] = -a[ 6];    a[22] = -a[22];
  212.     a[12] = -a[12];    a[28] = -a[28];
  213.     a[13] = -a[13];    a[29] = -a[29];
  214.     a[15] = -a[15];    a[31] = -a[31];
  215.     a[ 3] = -a[ 3];    a[19] = -a[19];
  216.     a[11] = -a[11];    a[27] = -a[27];
  217. #else
  218.     a[31] -= a[29]; a[29] -= a[27];
  219.     a[27] += a[25]; a[25] += a[23];
  220.     a[23] -= a[21]; a[21] += a[19];
  221.     a[19] -= a[17]; a[17] += a[15];
  222.     a[15] -= a[13]; a[13] += a[11];
  223.     a[11] -= a[ 9]; a[ 9] += a[ 7];
  224.     a[ 7] += a[ 5]; a[ 5] -= a[ 3];
  225.     a[ 3] -= a[ 1];
  226.  
  227.     a[31] += a[27]; a[27] = -(a[27]+a[23]);
  228.     a[30] -= a[26]; a[26] = a[22] - a[26];
  229.     a[23] -= a[19]; a[19] -= a[15];
  230.     a[22] += a[18]; a[18] += a[14];
  231.     a[15] -= a[11]; a[11] -= a[ 7];
  232.     a[14] += a[10]; a[10] -= a[ 6];
  233.     a[ 7] -= a[ 3];
  234.     a[ 6] -= a[ 2];
  235.  
  236.     a[31] = -(a[23]+a[31]); a[23] += a[15];    a[15] = -(a[15]+a[ 7]);
  237.     a[30] += a[22]; a[22] = -(a[14]+a[22]);    a[14] -= a[ 6];
  238.     a[29] += a[21]; a[21] = -(a[21]+a[13]);    a[13] -= a[ 5];
  239.     a[28] = -(a[28]+a[20]); a[20] += a[12];    a[12] = -(a[ 4]+a[12]);
  240. #endif
  241.     {
  242.     FLOAT8 const *xp = costab;
  243.     int p0,p1;
  244.     for (i = 0; i < 8; i++) {
  245.         FLOAT8 x1, x2, x3, x4;
  246.  
  247.         x3 = a[i+16] * (SQRT2*0.5);
  248.         x4 = a[i] - x3;
  249.         x3 = a[i] + x3;
  250.  
  251.         x2 = -(a[24+i] + a[i+8]) * (SQRT2*0.5);
  252.         x1 = (a[i+8] - x2) * xp[0];
  253.         x2 = (a[i+8] + x2) * xp[1];
  254.  
  255.         a[i   ] = x3 + x1;
  256.         a[i+ 8] = x4 - x2;
  257.         a[i+16] = x4 + x2;
  258.         a[i+24] = x3 - x1;
  259.     }
  260.  
  261.     xp += 2;
  262.  
  263.     for (i = 0; i < 4; i++) {
  264.         FLOAT8 xr;
  265.  
  266.         xr = a[i+28] * xp[0];
  267.         a[i+28] = (a[i] - xr);
  268.         a[i   ] = (a[i] + xr);
  269.  
  270.         xr = a[i+4] * xp[1];
  271.         a[i+ 4] = (a[i+24] - xr);
  272.         a[i+24] = (a[i+24] + xr);
  273.  
  274.         xr = a[i+20] * xp[2];
  275.         a[i+20] = (a[i+8] - xr);
  276.         a[i+ 8] = (a[i+8] + xr);
  277.  
  278.         xr = a[i+12] * xp[3];
  279.         a[i+12] = (a[i+16] - xr);
  280.         a[i+16] = (a[i+16] + xr);
  281.     }
  282.     xp += 4;
  283.  
  284.     for (i = 0; i < 4; i++) {
  285.         FLOAT8 xr;
  286.  
  287.         xr = a[30-i*4] * xp[0];
  288.         a[30-i*4] = (a[i*4] - xr);
  289.         a[   i*4] = (a[i*4] + xr);
  290.  
  291.         xr = a[ 2+i*4] * xp[1];
  292.         a[ 2+i*4] = (a[28-i*4] - xr);
  293.         a[28-i*4] = (a[28-i*4] + xr);
  294.  
  295.         xr = a[31-i*4] * xp[0];
  296.         a[31-i*4] = (a[1+i*4] - xr);
  297.         a[ 1+i*4] = (a[1+i*4] + xr);
  298.  
  299.         xr = a[ 3+i*4] * xp[1];
  300.         a[ 3+i*4] = (a[29-i*4] - xr);
  301.         a[29-i*4] = (a[29-i*4] + xr);
  302.  
  303.         xp += 2;
  304.     }
  305.  
  306.     p0 = 30;
  307.     p1 = 1;
  308.     do {
  309.         FLOAT8 xr = a[p1] * *xp++;
  310.         a[p1] = (a[p0] - xr);
  311.         a[p0] = (a[p0] + xr);
  312.         p0 -= 2; p1 += 2;
  313.     } while (p0 >= 0);
  314.     }
  315. }
  316.  
  317. INLINE static void window_subband(short *x1, FLOAT8 a[SBLIMIT])
  318. {
  319.     int i;
  320.     FLOAT8 const *wp = enwindow;
  321.  
  322.     short *x2 = &x1[238-14-286];
  323.  
  324.     for (i = -15; i < 0; i++) {
  325.     FLOAT8 w, s, t;
  326.  
  327.     w = wp[ 0]; s  = x2[  32] * w; t  = x1[- 32] * w;
  328.     w = wp[ 1]; s += x2[  96] * w; t += x1[- 96] * w;
  329.     w = wp[ 2]; s += x2[ 160] * w; t += x1[-160] * w;
  330.     w = wp[ 3]; s += x2[ 224] * w; t += x1[-224] * w;
  331.     w = wp[ 4]; s += x2[-224] * w; t += x1[ 224] * w;
  332.     w = wp[ 5]; s += x2[-160] * w; t += x1[ 160] * w;
  333.     w = wp[ 6]; s += x2[- 96] * w; t += x1[  96] * w;
  334.     w = wp[ 7]; s += x2[- 32] * w; t += x1[  32] * w;
  335.  
  336.     w = wp[ 8]; s += x1[-256] * w; t -= x2[ 256] * w;
  337.     w = wp[ 9]; s += x1[-192] * w; t -= x2[ 192] * w;
  338.     w = wp[10]; s += x1[-128] * w; t -= x2[ 128] * w;
  339.     w = wp[11]; s += x1[- 64] * w; t -= x2[  64] * w;
  340.     w = wp[12]; s -= x1[   0] * w; t += x2[   0] * w;
  341.     w = wp[13]; s -= x1[  64] * w; t += x2[- 64] * w;
  342.     w = wp[14]; s -= x1[ 128] * w; t += x2[-128] * w;
  343.     w = wp[15]; s -= x1[ 192] * w; t += x2[-192] * w;
  344.     wp += 16;
  345.  
  346.     *--a = t;
  347.     a[(i+16)*2] = s;
  348.     x1--;
  349.     x2++;
  350.     }
  351.     {
  352.     FLOAT8 s,t;
  353.     t  =  x1[- 16] * wp[0];               s  = x1[ -32] * wp[ 8];
  354.     t += (x1[- 48] - x1[ 16]) * wp[1];    s += x1[ -96] * wp[ 9];
  355.     t += (x1[- 80] + x1[ 48]) * wp[2];    s += x1[-160] * wp[10];
  356.     t += (x1[-112] - x1[ 80]) * wp[3];    s += x1[-224] * wp[11];
  357.     t += (x1[-144] + x1[112]) * wp[4];    s -= x1[  32] * wp[12];
  358.     t += (x1[-176] - x1[144]) * wp[5];    s -= x1[  96] * wp[13];
  359.     t += (x1[-208] + x1[176]) * wp[6];    s -= x1[ 160] * wp[14];
  360.     t += (x1[-240] - x1[208]) * wp[7];    s -= x1[ 224] * wp[15];
  361.  
  362.     a[15] = t;
  363.     a[-1] = s;
  364.     }
  365. }
  366.  
  367.  
  368. /*-------------------------------------------------------------------*/
  369. /*                                                                   */
  370. /*   Function: Calculation of the MDCT                               */
  371. /*   In the case of long blocks (type 0,1,3) there are               */
  372. /*   36 coefficents in the time domain and 18 in the frequency       */
  373. /*   domain.                                                         */
  374. /*   In the case of short blocks (type 2) there are 3                */
  375. /*   transformations with short length. This leads to 12 coefficents */
  376. /*   in the time and 6 in the frequency domain. In this case the     */
  377. /*   results are stored side by side in the vector out[].            */
  378. /*                                                                   */
  379. /*   New layer3                                                      */
  380. /*                                                                   */
  381. /*-------------------------------------------------------------------*/
  382.  
  383. INLINE static void mdct_short(FLOAT8 *out, FLOAT8 *in)
  384. {
  385.     int l;
  386.     for ( l = 0; l < 3; l++ ) {
  387.     FLOAT tc0,tc1,tc2,ts0,ts1,ts2;
  388.  
  389.     ts0 = in[5];
  390.     tc0 = in[3];
  391.     tc1 = ts0 + tc0;
  392.     tc2 = ts0 - tc0;
  393.  
  394.     ts0 = in[2];
  395.     tc0 = in[0];
  396.     ts1 = ts0 + tc0;
  397.     ts2 =-ts0 + tc0;
  398.  
  399.     tc0 = in[4];
  400.     ts0 = in[1];
  401.  
  402.     out[3*0] = tc1 + tc0;
  403.     out[3*5] =-ts1 + ts0;
  404.  
  405.     tc2 = tc2 * 0.86602540378443870761;
  406.     ts1 = ts1 * 0.5 + ts0;
  407.     out[3*1] = tc2-ts1;
  408.     out[3*2] = tc2+ts1;
  409.  
  410.     tc1 = tc1 * 0.5 - tc0;
  411.     ts2 = ts2 * 0.86602540378443870761;
  412.     out[3*3] = tc1+ts2;
  413.     out[3*4] = tc1-ts2;
  414.  
  415.     in += 6; out++;
  416.     }
  417. }
  418.  
  419. INLINE static void mdct_long(FLOAT8 *out, FLOAT8 *in)
  420. {
  421. #define inc(x) in[17-(x)]
  422. #define ins(x) in[8-(x)]
  423.  
  424.     const FLOAT8 c0=0.98480775301220802032, c1=0.64278760968653936292, c2=0.34202014332566882393;
  425.     const FLOAT8 c3=0.93969262078590842791, c4=-0.17364817766693030343, c5=-0.76604444311897790243;
  426.     FLOAT8 tc1 = inc(0)-inc(8),tc2 = (inc(1)-inc(7))*0.86602540378443870761, tc3 = inc(2)-inc(6), tc4 = inc(3)-inc(5);
  427.     FLOAT8 tc5 = inc(0)+inc(8),tc6 = (inc(1)+inc(7))*0.5,                    tc7 = inc(2)+inc(6), tc8 = inc(3)+inc(5);
  428.     FLOAT8 ts1 = ins(0)-ins(8),ts2 = (ins(1)-ins(7))*0.86602540378443870761, ts3 = ins(2)-ins(6), ts4 = ins(3)-ins(5);
  429.     FLOAT8 ts5 = ins(0)+ins(8),ts6 = (ins(1)+ins(7))*0.5,                    ts7 = ins(2)+ins(6), ts8 = ins(3)+ins(5);
  430.     FLOAT8 ct,st;
  431.  
  432.     ct = tc5+tc7+tc8+inc(1)+inc(4)+inc(7);
  433.     out[0] = ct;
  434.  
  435.     ct = tc1*c0 + tc2 + tc3*c1 + tc4*c2;
  436.     st = -ts5*c4 + ts6 - ts7*c5 + ts8*c3 + ins(4);
  437.     out[1] = ct+st;
  438.     out[2] = ct-st;
  439.  
  440.     ct =  tc5*c3 + tc6 + tc7*c4 + tc8*c5 - inc(4);
  441.     st = ts1*c2 + ts2 + ts3*c0 + ts4*c1;
  442.     out[3] = ct+st;
  443.     out[4] = ct-st;
  444.  
  445.     ct = (tc1-tc3-tc4)*0.86602540378443870761;
  446.     st = (ts5+ts7-ts8)*0.5+ins(1)-ins(4)+ins(7);
  447.     out[5] = ct+st;
  448.     out[6] = ct-st;
  449.  
  450.     ct = -tc5*c5 - tc6 - tc7*c3 - tc8*c4 + inc(4);
  451.     st = ts1*c1 + ts2 - ts3*c2 - ts4*c0;
  452.     out[7] = ct+st;
  453.     out[8] = ct-st;
  454.  
  455.     ct = tc1*c1 - tc2 - tc3*c2 + tc4*c0;
  456.     st = -ts5*c5 + ts6 - ts7*c3 + ts8*c4 + ins(4);
  457.     out[ 9] = ct+st;
  458.     out[10] = ct-st;
  459.  
  460.     ct = (tc5+tc7+tc8)*0.5-inc(1)-inc(4)-inc(7);
  461.     st = (ts1-ts3+ts4)*0.86602540378443870761;
  462.     out[11] = ct+st;
  463.     out[12] = ct-st;
  464.  
  465.     ct = tc1*c2 - tc2 + tc3*c0 - tc4*c1;
  466.     st =  ts5*c3 - ts6 + ts7*c4 - ts8*c5 - ins(4);
  467.     out[13] = ct+st;
  468.     out[14] = ct-st;
  469.  
  470.     ct = -tc5*c4 - tc6 - tc7*c5 - tc8*c3 + inc(4);
  471.     st = ts1*c0 - ts2 + ts3*c1 - ts4*c2;
  472.     out[15] = ct+st;
  473.     out[16] = ct-st;
  474.  
  475.     st = ts5+ts7-ts8-ins(1)+ins(4)-ins(7);
  476.     out[17] = st;
  477. }
  478.  
  479. static const int order[] = {
  480.     0,  16,  8, 24,  4,  20,  12,  28,
  481.     2,  18, 10, 26,  6,  22,  14,  30,
  482.     1,  17,  9, 25,  5,  21,  13,  29,
  483.     3,  19, 11, 27,  7,  23,  15,  31
  484. };
  485.  
  486.  
  487.  
  488. void mdct_init48(lame_global_flags *gfp)
  489. {
  490.     int i, k;
  491.     FLOAT8 sq;
  492.  
  493.     /* prepare the aliasing reduction butterflies */
  494.     for (k = 0; k < 8; k++) {
  495.     /*
  496.       This is table B.9: coefficients for aliasing reduction
  497.       */
  498.     static const FLOAT8 c[8] = {
  499.         -0.6,-0.535,-0.33,-0.185,-0.095,-0.041,-0.0142, -0.0037
  500.     };
  501.     sq = 1.0 + c[k] * c[k];
  502.     sq = sqrt(sq);
  503.     ca[k] = c[k] / sq;
  504.     cs[k] = 1.0 / sq;
  505.     }
  506.  
  507.     /* type 0*/
  508.     for (i = 0; i < 36; i++)
  509.     win[0][i] = sin(PI/36 * (i + 0.5));
  510.     /* type 1*/
  511.     for (i = 0; i < 18; i++) 
  512.     win[1][i] = win[0][i];
  513.     for (; i < 24; i++)
  514.     win[1][i] = 1.0;
  515.     for (; i < 30; i++)
  516.     win[1][i] = cos(PI/12 * (i + 0.5));
  517.     for (; i < 36; i++)
  518.     win[1][i] = 0.0;
  519.     /* type 3*/
  520.     for (i = 0; i < 36; i++)
  521.     win[3][i] = win[1][35 - i];
  522.  
  523.     /* swap window data*/
  524.     for (k = 0; k < 4; k++) {
  525.     FLOAT8 a;
  526.  
  527.     a = win[0][17-k];
  528.     win[0][17-k] = win[0][9+k];
  529.     win[0][9+k] = a;
  530.  
  531.     a = win[0][35-k];
  532.     win[0][35-k] = win[0][27+k];
  533.     win[0][27+k] = a;
  534.  
  535.     a = win[1][17-k];
  536.     win[1][17-k] = win[1][9+k];
  537.     win[1][9+k] = a;
  538.  
  539.     a = win[1][35-k];
  540.     win[1][35-k] = win[1][27+k];
  541.     win[1][27+k] = a;
  542.  
  543.     a = win[3][17-k];
  544.     win[3][17-k] = win[3][9+k];
  545.     win[3][9+k] = a;
  546.  
  547.     a = win[3][35-k];
  548.     win[3][35-k] = win[3][27+k];
  549.     win[3][27+k] = a;
  550.     }
  551.  
  552.     for (i = 0; i < NL; i++) {
  553.     win[0][i] *= cos((NL/4+0.5+i%9)*PI/NL) / SCALE/(NL/4);
  554.     win[1][i] *= cos((NL/4+0.5+i%9)*PI/NL) / SCALE/(NL/4);
  555.     win[3][i] *= cos((NL/4+0.5+i%9)*PI/NL) / SCALE/(NL/4);
  556.     }
  557.  
  558.     for (i = NL/2; i < NL; i++) {
  559.     win[0][i] *= -1;
  560.     win[1][i] *= -1;
  561.     win[3][i] *= -1;
  562.     }
  563.  
  564.     for (i = 0; i < NL/4; i++)
  565.     tantab_l[i] = tan((NL/4+0.5+i)*PI/NL);
  566.  
  567.     /* type 2(short)*/
  568.     for (i = 0; i < NS / 4; i++) {
  569.     FLOAT8 w2 = cos(PI / NS * (i + 0.5)) * (4.0/NS) / SCALE;
  570.     win[SHORT_TYPE][i] = tan(PI / NS * (i + 0.5));
  571.     tritab_s[i*2  ] = cos((0.5+2-i)*PI/NS) * w2;
  572.     tritab_s[i*2+1] = sin((0.5+2-i)*PI/NS) * w2;
  573.     }
  574. }
  575.  
  576. void mdct_sub48(lame_global_flags *gfp,
  577.     short *w0, short *w1,
  578.     FLOAT8 mdct_freq[2][2][576],
  579.     III_side_info_t *l3_side)
  580. {
  581.     int gr, k, ch;
  582.     short *wk;
  583.     static int init = 0;
  584.     lame_internal_flags *gfc=gfp->internal_flags;
  585.  
  586.     FLOAT8 work[18];
  587.  
  588.     if ( gfc->mdct_sub48_init == 0 ) {
  589.         gfc->mdct_sub48_init=1;
  590.     mdct_init48(gfp);
  591.     init++;
  592.     }
  593.  
  594.     wk = w0 + 286;
  595.     /* thinking cache performance, ch->gr loop is better than gr->ch loop */
  596.     for (ch = 0; ch < gfc->stereo; ch++) {
  597.     for (gr = 0; gr < gfc->mode_gr; gr++) {
  598.         int    band;
  599.         FLOAT8 *mdct_enc = &mdct_freq[gr][ch][0];
  600.         gr_info *gi = &(l3_side->gr[gr].ch[ch].tt);
  601.         FLOAT8 *samp = gfc->sb_sample[ch][1 - gr][0];
  602.  
  603.         for (k = 0; k < 18 / 2; k++) {
  604.         window_subband(wk, samp + 16);
  605.         idct32(samp);
  606.         window_subband(wk + 32, samp + 32+16);
  607.         idct32(samp+32);
  608.         samp += 64;
  609.         wk += 64;
  610.         /*
  611.          * Compensate for inversion in the analysis filter
  612.          */
  613.         for (band = 16-32; band < 0; band++)
  614.             samp[band] *= -1;
  615.         }
  616.  
  617.  
  618.         /* apply filters on the polyphase filterbank outputs */
  619.         /* bands <= gfc->highpass_band will be zeroed out below */
  620.         /* bands >= gfc->lowpass_band  will be zeroed out below */
  621.         if (gfc->filter_type==0) {
  622.               for (band=gfc->highpass_start_band;  band <= gfc->highpass_end_band; band++) { 
  623.           for (k=0; k<18; k++) 
  624.             gfc->sb_sample[ch][1-gr][k][order[band]]*=gfc->amp_highpass[band];
  625.           }
  626.               for (band=gfc->lowpass_start_band;  band <= gfc->lowpass_end_band; band++) { 
  627.           for (k=0; k<18; k++) 
  628.             gfc->sb_sample[ch][1-gr][k][order[band]]*=gfc->amp_lowpass[band];
  629.           }
  630.         }
  631.         
  632.  
  633.  
  634.         /*
  635.          * Perform imdct of 18 previous subband samples
  636.          * + 18 current subband samples
  637.          */
  638.         for (band = 0; band < 32; band++, mdct_enc += 18) 
  639.         {
  640.         int type = gi->block_type;
  641.         int band_swapped;
  642.         band_swapped = order[band];
  643. #ifdef ALLOW_MIXED
  644.         if (gi->mixed_block_flag && band < 2)
  645.             type = 0;
  646. #endif
  647.         if (band >= gfc->lowpass_band || band <= gfc->highpass_band) {
  648.             memset((char *)mdct_enc,0,18*sizeof(FLOAT8));
  649.         }else {
  650.           if (type == SHORT_TYPE) {
  651.             for (k = 2; k >= 0; --k) {
  652.             FLOAT8 win1 = win[SHORT_TYPE][k];
  653.             FLOAT8 a, b;
  654.  
  655.             a = gfc->sb_sample[ch][gr][k+6][band_swapped] * win1 -
  656.                 gfc->sb_sample[ch][gr][11-k][band_swapped];
  657.  
  658.             b = gfc->sb_sample[ch][gr][k+12][band_swapped] +
  659.                 gfc->sb_sample[ch][gr][17-k][band_swapped] * win1;
  660.  
  661.             work[k+3] = -b*tritab_s[k*2  ] + a * tritab_s[k*2+1];
  662.             work[k  ] =  b*tritab_s[k*2+1] + a * tritab_s[k*2  ];
  663.  
  664.             a = gfc->sb_sample[ch][gr][k+12][band_swapped] * win1 -
  665.                 gfc->sb_sample[ch][gr][17-k][band_swapped];
  666.  
  667.             b = gfc->sb_sample[ch][1-gr][k][band_swapped] +
  668.                 gfc->sb_sample[ch][1-gr][5-k][band_swapped] * win1;
  669.  
  670.             work[k+9] = -b*tritab_s[k*2  ] + a * tritab_s[k*2+1];
  671.             work[k+6] =  b*tritab_s[k*2+1] + a * tritab_s[k*2  ];
  672.  
  673.             a = gfc->sb_sample[ch][1-gr][k][band_swapped] * win1 -
  674.                 gfc->sb_sample[ch][1-gr][5-k][band_swapped];
  675.  
  676.             b = gfc->sb_sample[ch][1-gr][k+6][band_swapped] +
  677.                 gfc->sb_sample[ch][1-gr][11-k][band_swapped] * win1;
  678.  
  679.             work[k+15] = -b*tritab_s[k*2  ] + a * tritab_s[k*2+1];
  680.             work[k+12] =  b*tritab_s[k*2+1] + a * tritab_s[k*2  ];
  681.             }
  682.             mdct_short(mdct_enc, work);
  683.           } else {
  684.             for (k = -NL/4; k < 0; k++) {
  685.             FLOAT8 a, b;
  686.             a = win[type][k+27] * gfc->sb_sample[ch][1-gr][k+9][band_swapped]
  687.               + win[type][k+36] * gfc->sb_sample[ch][1-gr][8-k][band_swapped];
  688.             b = win[type][k+ 9] * gfc->sb_sample[ch][gr][k+9][band_swapped]
  689.               - win[type][k+18] * gfc->sb_sample[ch][gr][8-k][band_swapped];
  690.             work[k+ 9] = a - b*tantab_l[k+9];
  691.             work[k+18] = a*tantab_l[k+9] + b;
  692.             }
  693.  
  694.             mdct_long(mdct_enc, work);
  695.           }
  696.         }
  697.         
  698.         
  699.         /*
  700.           Perform aliasing reduction butterfly
  701.         */
  702.         if (type != SHORT_TYPE) {
  703.           if (band == 0)
  704.             continue;
  705.           for (k = 7; k >= 0; --k) {
  706.             FLOAT8 bu,bd;
  707.             bu = mdct_enc[k] * ca[k] + mdct_enc[-1-k] * cs[k];
  708.             bd = mdct_enc[k] * cs[k] - mdct_enc[-1-k] * ca[k];
  709.             
  710.             mdct_enc[-1-k] = bu;
  711.             mdct_enc[k]    = bd;
  712.           }
  713.         }
  714.           }
  715.     }
  716.     wk = w1 + 286;
  717.     if (gfc->mode_gr == 1) {
  718.         memcpy(gfc->sb_sample[ch][0], gfc->sb_sample[ch][1], 576 * sizeof(FLOAT8));
  719.     }
  720.     }
  721. }
  722.