MDCT.C
上传用户:njqiyou
上传日期:2007-01-08
资源大小:574k
文件大小:6k
源码类别:

mpeg/mp3

开发平台:

C/C++

  1. /**********************************************************************
  2.  * ISO MPEG Audio Subgroup Software Simulation Group (1996)
  3.  * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension
  4.  *
  5.  * $Id: mdct.c,v 1.1 1996/02/14 04:04:23 rowlands Exp $
  6.  *
  7.  * $Log: mdct.c,v $
  8.  * Revision 1.1  1996/02/14 04:04:23  rowlands
  9.  * Initial revision
  10.  *
  11.  * Received from Mike Coleman
  12.  **********************************************************************/
  13. #include "common.h"
  14. #include "l3side.h"
  15. #include "mdct.h"
  16. double ca[8], cs[8];
  17. /*
  18.   This is table B.9: coefficients for aliasing reduction
  19. */
  20. static double c[8] = { -0.6,-0.535,-0.33,-0.185,-0.095,-0.041,-0.0142, -0.0037 };
  21. void mdct_sub( L3SBS (*sb_sample), double (*mdct_freq)[2][576], int stereo, III_side_info_t *l3_side, int mode_gr )
  22. {
  23.     gr_info *cod_info;
  24.     double mdct_in[36];
  25.     int ch,gr,band,k,j;
  26.     double bu,bd;
  27.     static int init = 0;
  28.     int block_type;
  29.     double (*mdct_enc)[2][32][18] = (double (*)[2][32][18]) mdct_freq;
  30.     
  31.     if ( init == 0 )
  32.     {
  33. /* prepare the aliasing reduction butterflies */
  34. for ( k = 0; k < 8; k++ )
  35. {
  36.     double sq;
  37.     sq = sqrt( 1.0 + c[k] * c[k] );
  38.     ca[k] = c[k] / sq;
  39.     cs[k] = 1.0 / sq;
  40. }
  41. init++;
  42.     }
  43.     
  44.     for ( gr = 0; gr < mode_gr; gr++ )
  45. for ( ch = 0; ch < stereo; ch++ )
  46. {
  47.     cod_info = (gr_info *) &(l3_side->gr[gr].ch[ch]) ;
  48.     block_type = cod_info->block_type;
  49.     
  50.     /*
  51.       Compensate for inversion in the analysis filter
  52.     */
  53.     for ( band = 0; band < 32; band++ )
  54. for ( k = 0; k < 18; k++ )
  55.     if ( (band & 1) && (k & 1) )
  56. (*sb_sample)[ch][gr+1][k][band] *= -1.0;
  57.     
  58.     /*
  59.       Perform imdct of 18 previous subband samples
  60.       + 18 current subband samples
  61.     */
  62.     for ( band = 0; band < 32; band++ )
  63.     {
  64. for ( k = 0; k < 18; k++ )
  65. {
  66.     mdct_in[k]    = (*sb_sample)[ch][ gr ][k][band];
  67.     mdct_in[k+18] = (*sb_sample)[ch][gr+1][k][band];
  68. }
  69. if ( cod_info->mixed_block_flag && (band < 2) )
  70.     block_type = 0;
  71. mdct( mdct_in, &mdct_enc[gr][ch][band][0], block_type );
  72.     }
  73.     
  74.     /*
  75.       Perform aliasing reduction butterfly
  76.       on long blocks
  77.     */
  78.     if ( block_type != 2 )
  79. for ( band = 0; band < 31; band++ )
  80.     for ( k = 0; k < 8; k++ )
  81.     {
  82. bu = mdct_enc[gr][ch][band][17-k] * cs[k] + mdct_enc[gr][ch][band+1][k] * ca[k];
  83. bd = mdct_enc[gr][ch][band+1][k] * cs[k] - mdct_enc[gr][ch][band][17-k] * ca[k];
  84. mdct_enc[gr][ch][band][17-k] = bu;
  85. mdct_enc[gr][ch][band+1][k]  = bd;
  86.     }
  87.     
  88. }
  89.     
  90.     /*
  91.       Save latest granule's subband samples to be used in
  92.       the next mdct call
  93.     */
  94.     for ( ch = 0; ch < stereo; ch++ )
  95. for ( j = 0; j < 18; j++ )
  96.     for ( band = 0; band < 32; band++ )
  97. (*sb_sample)[ch][0][j][band] = (*sb_sample)[ch][mode_gr][j][band];
  98. }
  99. void mdct( double *in, double *out, int block_type )
  100. {
  101. /*-------------------------------------------------------------------*/
  102. /*                                                                   */
  103. /*   Function: Calculation of the MDCT                               */
  104. /*   In the case of long blocks ( block_type 0,1,3 ) there are       */
  105. /*   36 coefficents in the time domain and 18 in the frequency       */
  106. /*   domain.                                                         */
  107. /*   In the case of short blocks (block_type 2 ) there are 3         */
  108. /*   transformations with short length. This leads to 12 coefficents */
  109. /*   in the time and 6 in the frequency domain. In this case the     */
  110. /*   results are stored side by side in the vector out[].            */
  111. /*                                                                   */
  112. /*   New layer3                                                      */
  113. /*                                                                   */
  114. /*-------------------------------------------------------------------*/
  115.   int l,k,i,m,N;
  116.   double sum;
  117.   static double win[4][36];
  118.   static int init = 0;
  119.   static double cos_s[6][12], cos_l[18][36];
  120.   if ( init == 0 )
  121.   {
  122.     /* type 0 */
  123.     for ( i = 0; i < 36; i++ )
  124.       win[0][i] = sin( PI/36 * (i + 0.5) );
  125.     /* type 1*/
  126.     for ( i = 0; i < 18; i++ ) 
  127.       win[1][i] = sin( PI/36 * (i + 0.5) );
  128.     for ( i = 18; i < 24; i++ )
  129.       win[1][i] = 1.0;
  130.     for ( i = 24; i < 30; i++ )
  131.       win[1][i] = sin( PI/12 * ( i + 0.5 - 18) );
  132.     for ( i = 30; i < 36; i++ )
  133.       win[1][i] = 0.0;
  134.     /* type 3*/
  135.     for ( i = 0; i < 6; i++ )
  136.       win[3][i] = 0.0;
  137.     for ( i = 6; i < 12; i++ ) 
  138.       win[3][i] = sin( PI/12 * (i + 0.5 - 6) );
  139.     for ( i = 12; i < 18; i++ )
  140.       win[3][i] = 1.0;
  141.     for ( i = 18; i < 36; i++ )
  142.       win[3][i] = sin( PI/36 * (i + 0.5) );
  143.     /* type 2*/
  144.     for ( i = 0; i < 12; i++ )
  145.     win[2][i] = sin( PI/12 * (i + 0.5) );
  146.     for ( i = 12; i < 36; i++ )
  147.       win[2][i] = 0.0;
  148.     N = 12;
  149.     for ( m = 0; m < N / 2; m++ )
  150.       for ( k = 0; k < N; k++ )
  151.         cos_s[m][k] = cos( (PI /(2 * N)) * (2 * k + 1 + N / 2) *
  152.                      (2 * m + 1) ) / (N / 4);
  153.     N = 36;
  154.     for ( m = 0; m < N / 2; m++ )
  155.       for ( k = 0; k < N; k++ )
  156.         cos_l[m][k] = cos( (PI / (2 * N)) * (2 * k + 1 + N / 2) *
  157.                      (2 * m + 1) ) / (N / 4);
  158.     init++;
  159.   }
  160.   if ( block_type == 2 )
  161.   {
  162.     N = 12;
  163.     for ( l = 0; l < 3; l++ )
  164.     {
  165.       for ( m = 0; m < N / 2; m++ )
  166.       {
  167.         for ( sum = 0.0, k = 0; k < N; k++ )
  168.           sum += win[block_type][k] * in[k + 6 * l + 6] * cos_s[m][k];
  169.         out[ 3 * m + l] = sum;
  170.       }
  171.     }
  172.   }
  173.   else
  174.   {
  175.     N = 36;
  176.     for ( m = 0; m < N / 2; m++ )
  177.     {
  178.       for ( sum = 0.0, k = 0; k < N; k++ )
  179.         sum += win[block_type][k] * in[k] * cos_l[m][k];
  180.       out[m] = sum;
  181.     }
  182.   }
  183. }
  184. void
  185. delay( double (*xr)[2][576], int stereo )
  186. {
  187.     static double xr_buff[2][576];
  188.     double xr_buff2[2][576];
  189.     unsigned int i,j;
  190.     
  191.     for (i=0;i<stereo;i++)
  192.     {
  193. for (j=0;j<576;j++) xr_buff2[i][j] = xr_buff[i][j];
  194. for (j=0;j<576;j++) xr_buff[i][j]  = xr[1][i][j];
  195. for (j=0;j<576;j++) xr[1][i][j]    = xr[0][i][j];
  196. for (j=0;j<576;j++) xr[0][i][j]    = xr_buff2[i][j];
  197.     }
  198. }