home *** CD-ROM | disk | FTP | other *** search
/ PC Welt 2001 September / PC-WELT 9-2001.ISO / software / hw / brennen / flask_src.exe / Audio / MPEG / MPEG_imdct.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-06  |  9.9 KB  |  420 lines
  1. /* 
  2.  *  MPEG_imdct.cpp
  3.  *
  4.  *  Code from
  5.  *            NekoAmp 1.3 decoder by Avery Lee
  6.  *
  7.  *  FlasKMPEG
  8.  *    Copyright (C) Alberto Vigata - January 2000
  9.  *
  10.  *  This file is part of FlasKMPEG, a free MPEG to MPEG/AVI converter
  11.  *    
  12.  *  FlasKMPEG is free software; you can redistribute it and/or modify
  13.  *  it under the terms of the GNU General Public License as published by
  14.  *  the Free Software Foundation; either version 2, or (at your option)
  15.  *  any later version.
  16.  *   
  17.  *  FlasKMPEG is distributed in the hope that it will be useful,
  18.  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  19.  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20.  *  GNU General Public License for more details.
  21.  *   
  22.  *  You should have received a copy of the GNU General Public License
  23.  *  along with GNU Make; see the file COPYING.  If not, write to
  24.  *  the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
  25.  *
  26.  */
  27.  
  28.  
  29. #include <crtdbg.h>
  30. #include <math.h>
  31. #include <string.h>
  32.  
  33. #define M_PI (3.141592653589793238462643)
  34. #define PI (3.141592653589793238462643)
  35.  
  36. //#define MDCT_BLOCKTYPE_PROFILE
  37.  
  38. #ifdef MDCT_BLOCKTYPE_PROFILE
  39.     #include <windows.h>
  40. #endif
  41.  
  42. static double win[4][36];
  43.  
  44. static double w[18];
  45. static double w2[9];
  46. static double coef[9][4];
  47.  
  48. static double v[6];
  49. static double v2[3];
  50. static const double coef87 = 1.7320508075689;
  51.  
  52. void init_mdct()
  53. {
  54.     int k, p, n;
  55.     double t, pi;
  56.     
  57.     /*--- 18 point --*/
  58.     n = 18;
  59.     pi = 4.0 * atan(1.0);
  60.     t = pi / (4 * n);
  61.     for (p = 0; p < n; p++)
  62.         w[p] = (double) (2.0 * cos(t * (2 * p + 1)));
  63.     for (p = 0; p < 9; p++)
  64.         w2[p] = (double) 2.0 *cos(2 * t * (2 * p + 1));
  65.     
  66.     t = pi / (2 * n);
  67.     for (k = 0; k < 9; k++)
  68.     {
  69.         for (p = 0; p < 4; p++)
  70.             coef[k][p] = (double) cos(t * (2 * k) * (2 * p + 1));
  71.     }
  72.     
  73.     /*--- 6 point */
  74.     n = 6;
  75.     pi = 4.0 * atan(1.0);
  76.     t = pi / (4 * n);
  77.     for (p = 0; p < n; p++)
  78.         v[p] = (double) 2.0 *cos(t * (2 * p + 1));
  79.     
  80.     for (p = 0; p < 3; p++)
  81.         v2[p] = (double) 2.0 *cos(2 * t * (2 * p + 1));
  82.     
  83. #if 0
  84.     t = pi / (2 * n);
  85.     k = 1;
  86.     p = 0;
  87.     coef87 = (double) cos(t * (2 * k) * (2 * p + 1));
  88. #endif
  89.     
  90.     /* adjust scaling to save a few mults */
  91.     for (p = 0; p < 6; p++)
  92.         v[p] = v[p] / 2.0f;
  93.  
  94.     //coef87 = (double) 2.0 * coef87;
  95.     
  96.     int i;
  97.  
  98.     /* type 0  - invert for speed */
  99.       for(i=0;i<36;i++)
  100.          win[0][i] = -sin( PI/36 *(i+0.5) );
  101.  
  102.     /* type 1 - invert for speed*/
  103.       for(i=0;i<18;i++)
  104.          win[1][i] = -sin( PI/36 *(i+0.5) );
  105. //      for(i=18;i<24;i++)
  106. //         win[1][i] = 1.0;
  107.       for(i=24;i<30;i++)
  108.          win[1][i] = -sin( PI/12 *(i+0.5-18) );
  109. //      for(i=30;i<36;i++)
  110. //         win[1][i] = 0.0;
  111.  
  112.     /* type 3 - invert for speed*/
  113. //      for(i=0;i<6;i++)
  114. //         win[3][i] = 0.0;
  115.       for(i=6;i<12;i++)
  116.          win[3][i] = -sin( PI/12 *(i+0.5-6) );
  117. //      for(i=12;i<18;i++)
  118. //         win[3][i] =1.0;
  119.       for(i=18;i<36;i++)
  120.          win[3][i] = -sin( PI/36*(i+0.5) );
  121.  
  122.     /* type 2*/
  123.       for(i=0;i<12;i++)
  124.          win[2][i] = sin( PI/12*(i+0.5) ) ;
  125. //      for(i=12;i<36;i++)
  126. //         win[2][i] = 0.0 ;
  127.     
  128.     return;
  129. }
  130.  
  131. void imdct18(double f[18], float src[18])    /* 18 point */
  132. {
  133.    int p;
  134.    double a[9], b[9];
  135.    double ap, bp, a8p, b8p;
  136.    double g1, g2;
  137.  
  138.  
  139.    for (p = 0; p < 4; p++)
  140.    {
  141.       g1 = w[p] * src[p];
  142.       g2 = w[17 - p] * src[17 - p];
  143.       ap = g1 + g2;        // a[p]
  144.  
  145.       bp = w2[p] * (g1 - g2);    // b[p]
  146.  
  147.       g1 = w[8 - p] * src[8 - p];
  148.       g2 = w[9 + p] * src[9 + p];
  149.       a8p = g1 + g2;        // a[8-p]
  150.  
  151.       b8p = w2[8 - p] * (g1 - g2);    // b[8-p]
  152.  
  153.       a[p] = ap + a8p;
  154.       a[5 + p] = ap - a8p;
  155.       b[p] = bp + b8p;
  156.       b[5 + p] = bp - b8p;
  157.    }
  158.    g1 = w[p] * src[p];
  159.    g2 = w[17 - p] * src[17 - p];
  160.    a[p] = g1 + g2;
  161.    b[p] = w2[p] * (g1 - g2);
  162.  
  163.  
  164.    f[0] = 0.5f * (a[0] + a[1] + a[2] + a[3] + a[4]);
  165.    f[1] = 0.5f * (b[0] + b[1] + b[2] + b[3] + b[4]);
  166.  
  167.    f[2] = coef[1][0] * a[5] + coef[1][1] * a[6] + coef[1][2] * a[7] + coef[1][3] * a[8];
  168.    f[3] = coef[1][0] * b[5] + coef[1][1] * b[6] + coef[1][2] * b[7] + coef[1][3] * b[8] - f[1];
  169.    f[1] = f[1] - f[0];
  170.    f[2] = f[2] - f[1];
  171.  
  172.    f[4] = coef[2][0] * a[0] + coef[2][1] * a[1] + coef[2][2] * a[2] + coef[2][3] * a[3] - a[4];
  173.    f[5] = coef[2][0] * b[0] + coef[2][1] * b[1] + coef[2][2] * b[2] + coef[2][3] * b[3] - b[4] - f[3];
  174.    f[3] = f[3] - f[2];
  175.    f[4] = f[4] - f[3];
  176.  
  177.    f[6] = coef[3][0] * (a[5] - a[7] - a[8]);
  178.    f[7] = coef[3][0] * (b[5] - b[7] - b[8]) - f[5];
  179.    f[5] = f[5] - f[4];
  180.    f[6] = f[6] - f[5];
  181.  
  182.    f[8] = coef[4][0] * a[0] + coef[4][1] * a[1] + coef[4][2] * a[2] + coef[4][3] * a[3] + a[4];
  183.    f[9] = coef[4][0] * b[0] + coef[4][1] * b[1] + coef[4][2] * b[2] + coef[4][3] * b[3] + b[4] - f[7];
  184.    f[7] = f[7] - f[6];
  185.    f[8] = f[8] - f[7];
  186.  
  187.    f[10] = coef[5][0] * a[5] + coef[5][1] * a[6] + coef[5][2] * a[7] + coef[5][3] * a[8];
  188.    f[11] = coef[5][0] * b[5] + coef[5][1] * b[6] + coef[5][2] * b[7] + coef[5][3] * b[8] - f[9];
  189.    f[9] = f[9] - f[8];
  190.    f[10] = f[10] - f[9];
  191.  
  192.    f[12] = 0.5f * (a[0] + a[2] + a[3]) - a[1] - a[4];
  193.    f[13] = 0.5f * (b[0] + b[2] + b[3]) - b[1] - b[4] - f[11];
  194.    f[11] = f[11] - f[10];
  195.    f[12] = f[12] - f[11];
  196.  
  197.    f[14] = coef[7][0] * a[5] + coef[7][1] * a[6] + coef[7][2] * a[7] + coef[7][3] * a[8];
  198.    f[15] = coef[7][0] * b[5] + coef[7][1] * b[6] + coef[7][2] * b[7] + coef[7][3] * b[8] - f[13];
  199.    f[13] = f[13] - f[12];
  200.    f[14] = f[14] - f[13];
  201.  
  202.    f[16] = coef[8][0] * a[0] + coef[8][1] * a[1] + coef[8][2] * a[2] + coef[8][3] * a[3] + a[4];
  203.    f[17] = coef[8][0] * b[0] + coef[8][1] * b[1] + coef[8][2] * b[2] + coef[8][3] * b[3] + b[4] - f[15];
  204.    f[15] = f[15] - f[14];
  205.    f[16] = f[16] - f[15];
  206.    f[17] = f[17] - f[16];
  207.  
  208.  
  209.    return;
  210. }
  211. /*--------------------------------------------------------------------*/
  212. /* does 3, 6 pt dct.  changes order from f[i][window] c[window][i] */
  213. void imdct6_3(double f[], float src[])    /* 6 point */
  214. {
  215.    int w;
  216.    double buf[18];
  217.    double *a, *c;        // b[i] = a[3+i]
  218.  
  219.    double g1, g2;
  220.    double a02, b02;
  221.  
  222.    c = f;
  223.    a = buf;
  224.    for (w = 0; w < 3; w++)
  225.    {
  226.       g1 = v[0] * src[3 * 0];
  227.       g2 = v[5] * src[3 * 5];
  228.       a[0] = g1 + g2;
  229.       a[3 + 0] = v2[0] * (g1 - g2);
  230.  
  231.       g1 = v[1] * src[3 * 1];
  232.       g2 = v[4] * src[3 * 4];
  233.       a[1] = g1 + g2;
  234.       a[3 + 1] = v2[1] * (g1 - g2);
  235.  
  236.       g1 = v[2] * src[3 * 2];
  237.       g2 = v[3] * src[3 * 3];
  238.       a[2] = g1 + g2;
  239.       a[3 + 2] = v2[2] * (g1 - g2);
  240.  
  241.       a += 6;
  242.       src++;
  243.    }
  244.  
  245.    a = buf;
  246.    for (w = 0; w < 3; w++)
  247.    {
  248.       a02 = (a[0] + a[2]);
  249.       b02 = (a[3 + 0] + a[3 + 2]);
  250.       c[0] = a02 + a[1];
  251.       c[1] = b02 + a[3 + 1];
  252.       c[2] = coef87 * (a[0] - a[2]);
  253.       c[3] = coef87 * (a[3 + 0] - a[3 + 2]) - c[1];
  254.       c[1] = c[1] - c[0];
  255.       c[2] = c[2] - c[1];
  256.       c[4] = a02 - a[1] - a[1];
  257.       c[5] = b02 - a[3 + 1] - a[3 + 1] - c[3];
  258.       c[3] = c[3] - c[2];
  259.       c[4] = c[4] - c[3];
  260.       c[5] = c[5] - c[4];
  261.       a += 6;
  262.       c += 6;
  263.    }
  264.  
  265.    return;
  266. }
  267.  
  268. void inv_mdct(float *in0,
  269.               float out2[18][32],
  270.                 float *prevblock,
  271.                 int block_type) {
  272.     int i,m,p;
  273.  
  274. #ifdef MDCT_BLOCKTYPE_PROFILE
  275.     static int prof_bt[4]={0,0,0,0};
  276.     static int prof_bt_total=0;
  277.  
  278.     prof_bt[block_type]++;
  279.     prof_bt_total++;
  280.  
  281.     if (!(prof_bt_total&262143)) {
  282.         static char buf[256];
  283.  
  284.         wsprintf(buf, "%d MDCTs, %d%% type 0, %d%% type 1, %d%% type 2, %d%% type 3\n"
  285.             ,prof_bt_total
  286.             ,(prof_bt[0]*100)/prof_bt_total
  287.             ,(prof_bt[1]*100)/prof_bt_total
  288.             ,(prof_bt[2]*100)/prof_bt_total
  289.             ,(prof_bt[3]*100)/prof_bt_total);
  290.         OutputDebugString(buf);
  291.     }
  292. #endif
  293.  
  294.     double  sum;
  295.     double in[18], out[36];
  296.  
  297.     /*
  298.     for(p=0;p<36;p++) {
  299.         sum = 0.0;
  300.         for(m=0;m<18;m++)
  301.             sum += in[m] * cos( (PI/72) * (2*p+1+18)*(2*m+1));
  302.  
  303.         out[p] = sum;
  304.     }
  305.     */
  306.  
  307.     if (block_type == 2) {
  308.         int j;
  309.  
  310.         imdct6_3(in, in0);
  311.  
  312.         for(i=0;i<36;i++)
  313.             out[i]=0.0;
  314.  
  315.         for(i=0;i<3;i++){
  316.             for(p=0;p<3;p++) {
  317.                 out[6*i+6+p] += in[i*6+p+3] * win[2][p];
  318.                 out[6*i+11-p] += -in[i*6+p+3] * win[2][5-p];
  319.                 out[6*i+14-p] += -in[i*6+p] * win[2][3+p];
  320.                 out[6*i+15+p] += -in[i*6+p] * win[2][2-p];
  321.             }
  322.         }
  323.  
  324.         for(i=0; i<18; i++) {
  325.             out2[i][0] = out[i] + prevblock[i];
  326.             prevblock[i] = out[i+18];
  327.         }
  328.  
  329.     } else {
  330.  
  331.         imdct18(in, in0);
  332.  
  333.         switch(block_type) {
  334.         case 0:
  335. //            for(i=0; i<18; i++) {
  336. //                out[i   ] *= win[0][i];
  337. //                out[35-i] *= win[0][i];
  338. //            }
  339.  
  340. //            for(i=0; i<9; i++) {
  341. //                out[i] = in[i+9];
  342. //                out[17-i] = -in[i+9];
  343. //                out[26-i] = -in[i];
  344. //                out[27+i] = -in[i];
  345. //            }
  346.  
  347.             for(i=0; i<9; i++) {
  348.                 out2[i  ][0] = prevblock[i  ] - in[ 9+i]*win[0][i  ];
  349.                 out2[i+9][0] = prevblock[i+9] + in[17-i]*win[0][i+9];
  350.             }
  351.  
  352.             for(i=0; i<4; i++) {
  353.                 double in1 = in[8-i];
  354.                 double in2 = in[i];
  355.  
  356.                 prevblock[   i] = in1 * win[0][17-i];
  357.                 prevblock[ 9+i] = in2 * win[0][ 8-i];
  358.                 prevblock[ 8-i] = in2 * win[0][ 9+i];
  359.                 prevblock[17-i] = in1 * win[0][   i];
  360.             }
  361.  
  362.             prevblock[4  ] = in[4] * win[0][13];
  363.             prevblock[4+9] = in[4] * win[0][ 4];
  364.  
  365.             break;
  366.         case 1:
  367.             for(i=0; i<9; i++) {
  368.                 out[i] = -in[i+9];
  369.                 out[17-i] = in[i+9];
  370.                 out[26-i] = in[i];
  371.                 out[27+i] = in[i];
  372.             }
  373.  
  374.             for(i=0; i<18; i++)
  375.                 out[i] *= win[1][i];
  376.  
  377.             for(i=18; i<24; i++)
  378.                 out[i] = -out[i];
  379.  
  380.             for(i=24; i<30; i++)
  381.                 out[i] *= win[1][i];
  382.  
  383.             for(i=30; i<36; i++)
  384.                 out[i] = 0.0;
  385.  
  386.             for(i=0; i<18; i++) {
  387.                 out2[i][0] = out[i] + prevblock[i];
  388.                 prevblock[i] = out[i+18];
  389.             }
  390.             break;
  391.         case 3:
  392.             for(i=0; i<9; i++) {
  393.                 out[i] = -in[i+9];
  394.                 out[17-i] = in[i+9];
  395.                 out[26-i] = in[i];
  396.                 out[27+i] = in[i];
  397.             }
  398.  
  399.             for(i=0; i<6; i++)
  400.                 out[i] = 0.0;
  401.  
  402.             for(i=6; i<12; i++)
  403.                 out[i] *= win[3][i];
  404.  
  405.             for(i=12; i<18; i++)
  406.                 out[i] = -out[i];
  407.  
  408.             for(i=18; i<36; i++)
  409.                 out[i] *= win[3][i];
  410.  
  411.             for(i=0; i<18; i++) {
  412.                 out2[i][0] = out[i] + prevblock[i];
  413.                 prevblock[i] = out[i+18];
  414.             }
  415.             break;
  416.         }
  417.     }
  418. }
  419.  
  420.