开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > mpeg/mp3 > 查看源码
L3PSY.C
资源名称：mpeg.zip　[点击查看]
上传用户：njqiyou
上传日期：2007-01-08
资源大小：574k
文件大小：36k
源码类别：

mpeg/mp3

开发平台：

C/C++

L3PSY.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: l3psy.c,v 1.2 1997/01/19 22:28:29 rowlands Exp $
  6.  *
  7.  * $Log: l3psy.c,v $
  8.  * Revision 1.2  1997/01/19 22:28:29  rowlands
  9.  * Layer 3 bug fixes from Seymour Shlien
  10.  *
  11.  * Revision 1.1  1996/02/14 04:04:23  rowlands
  12.  * Initial revision
  13.  *
  14.  * Received from Mike Coleman
  15.  **********************************************************************/
  16. /**********************************************************************
  17.  *   date   programmers         comment                               *
  18.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  19.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  20.  * 7/10/91  Earle Jennings      Ported to MsDos.                      *
  21.  *                              replace of floats with FLOAT          *
  22.  * 2/11/92  W. Joseph Carter    Fixed mem_alloc() arg for "absthr".   *
  23.  * 3/16/92  Masahiro Iwadare Modification for Layer III            *
  24.  * 17/4/93  Masahiro Iwadare    Updated for IS Modification           *
  25.  **********************************************************************/
  27. #include "common.h"
  28. #include "encoder.h"
  29. #include "l3psy.h"
  30. #include "l3side.h"
  31. #include <assert.h>
  33. #define maximum(x,y) ( (x>y) ? x : y )
  34. #define minimum(x,y) ( (x<y) ? x : y )
  36. void L3para_read( double sfreq, int numlines[CBANDS], int partition_l[HBLKSIZE],
  37.   double minval[CBANDS], double qthr_l[CBANDS], double norm_l[CBANDS],
  38.   double s3_l[CBANDS][CBANDS], int partition_s[HBLKSIZE_s], double qthr_s[CBANDS_s],
  39.   double norm_s[CBANDS_s], double SNR_s[CBANDS_s],
  40.   int cbw_l[SBMAX_l], int bu_l[SBMAX_l], int bo_l[SBMAX_l],
  41.   double w1_l[SBMAX_l], double w2_l[SBMAX_l],
  42.   int cbw_s[SBMAX_s], int bu_s[SBMAX_s], int bo_s[SBMAX_s],
  43.   double w1_s[SBMAX_s], double w2_s[SBMAX_s] );
  45. void L3psycho_anal( short int *buffer, short int savebuf[1344], int chn, int lay, FLOAT snr32[32],
  46.     double sfreq, double ratio_d[21], double ratio_ds[12][3],
  47.     double *pe, gr_info *cod_info )
  48. {
  49.     static double ratio[2][21];
  50.     static double ratio_s[2][12][3];
  51.     int blocktype;
  52.     unsigned int   b, i, j, k;
  53.     double         r_prime, phi_prime; /* not FLOAT */
  54.     FLOAT          freq_mult, bval_lo, min_thres, sum_energy;
  55.     double         tb, temp1,temp2,temp3;
  57.     /*         nint(); Layer III */
  58.     double   thr[CBANDS];
  60. /* The static variables "r", "phi_sav", "new", "old" and "oldest" have    */
  61. /* to be remembered for the unpredictability measure.  For "r" and        */
  62. /* "phi_sav", the first index from the left is the channel select and     */
  63. /* the second index is the "age" of the data.                             */
  66.    static FLOAT window_s[BLKSIZE_s] ;
  67.  static int     new = 0, old = 1, oldest = 0;
  68.  static int     init = 0, flush, sync_flush, syncsize, sfreq_idx;
  69.  static double  cw[HBLKSIZE], eb[CBANDS];
  70.  static double  ctb[CBANDS];
  71.  static double SNR_l[CBANDS], SNR_s[CBANDS_s];
  72.  static int init_L3;
  73.  static double minval[CBANDS],qthr_l[CBANDS],norm_l[CBANDS];
  74.  static double qthr_s[CBANDS_s],norm_s[CBANDS_s];
  75.  static double nb_1[2][CBANDS], nb_2[2][CBANDS];
  76.  static double s3_l[CBANDS][CBANDS]; /* s3_s[CBANDS_s][CBANDS_s]; */
  78. /* Scale Factor Bands */
  79.  static int cbw_l[SBMAX_l],bu_l[SBMAX_l],bo_l[SBMAX_l] ;
  80.  static int cbw_s[SBMAX_s],bu_s[SBMAX_s],bo_s[SBMAX_s] ;
  81.  static double w1_l[SBMAX_l], w2_l[SBMAX_l];
  82.  static double w1_s[SBMAX_s], w2_s[SBMAX_s];
  83.  static double en[SBMAX_l],   thm[SBMAX_l] ;
  84.  static int blocktype_old[2] ;
  85.  int sb,sblock;
  86.  static int partition_l[HBLKSIZE],partition_s[HBLKSIZE_s];
  89. /* The following static variables are constants.                           */
  91.  static double  nmt = 5.5;
  93.  static FLOAT   crit_band[27] = {0,  100,  200, 300, 400, 510, 630,  770,
  94.                                920, 1080, 1270,1480,1720,2000,2320, 2700,
  95.                               3150, 3700, 4400,5300,6400,7700,9500,12000,
  96.                              15500,25000,30000};
  98.  static FLOAT   bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,
  99.                             10.0,  7.0,  4.4,  4.5,  4.5,  4.5,  4.5,
  100.                              4.5,  4.5,  4.5,  4.5,  4.5,  4.5,  4.5,
  101.                              4.5,  4.5,  4.5,  3.5,  3.5,  3.5};
  103. /* The following pointer variables point to large areas of memory         */
  104. /* dynamically allocated by the mem_alloc() function.  Dynamic memory     */
  105. /* allocation is used in order to avoid stack frame or data area          */
  106. /* overflow errors that otherwise would have occurred at compile time     */
  107. /* on the Macintosh computer.                                             */
  109.  FLOAT          *grouped_c, *grouped_e, *nb, *cb, *ecb, *bc;
  110.  FLOAT          *wsamp_r, *wsamp_i, *phi, *energy;
  111.  static FLOAT energy_s[3][256];
  112.  static FLOAT phi_s[3][256] ; /* 256 samples not 129 */
  113.  FLOAT          *c, *fthr;
  114.  F32            *snrtmp;
  116.  static int *numlines ;
  117.  static int     *partition;
  118.  static FLOAT   *cbval, *rnorm;
  119.  static FLOAT   *window;
  120.  static FLOAT   *absthr;
  121.  static double  *tmn;
  122.  static FCB     *s;
  123.  static FHBLK   *lthr;
  124.  static F2HBLK  *r, *phi_sav;
  126. /* These dynamic memory allocations simulate "automatic" variables        */
  127. /* placed on the stack.  For each mem_alloc() call here, there must be    */
  128. /* a corresponding mem_free() call at the end of this function.           */
  130.  grouped_c = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_c");
  131.  grouped_e = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_e");
  132.  nb = (FLOAT *) mem_alloc(sizeof(FCB), "nb");
  133.  cb = (FLOAT *) mem_alloc(sizeof(FCB), "cb");
  134.  ecb = (FLOAT *) mem_alloc(sizeof(FCB), "ecb");
  135.  bc = (FLOAT *) mem_alloc(sizeof(FCB), "bc");
  136.  wsamp_r = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_r");
  137.  wsamp_i = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_i");
  138.  phi = (FLOAT *) mem_alloc(sizeof(FBLK), "phi");
  139.  energy = (FLOAT *) mem_alloc(sizeof(FBLK), "energy");
  140.  c = (FLOAT *) mem_alloc(sizeof(FHBLK), "c");
  141.  fthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "fthr");
  142.  snrtmp = (F32 *) mem_alloc(sizeof(F2_32), "snrtmp");
  144.     assert( lay == 3 );
  145.  if(init==0){
  147. /* These dynamic memory allocations simulate "static" variables placed    */
  148. /* in the data space.  Each mem_alloc() call here occurs only once at     */
  149. /* initialization time.  The mem_free() function must not be called.      */
  150.      numlines = (int *) mem_alloc(sizeof(ICB), "numlines");
  151.      partition = (int *) mem_alloc(sizeof(IHBLK), "partition");
  152.      cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");
  153.      rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");
  154.      window = (FLOAT *) mem_alloc(sizeof(FBLK), "window");
  155.      absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr"); 
  156.      tmn = (double *) mem_alloc(sizeof(DCB), "tmn");
  157.      s = (FCB *) mem_alloc(sizeof(FCBCB), "s");
  158.      lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");
  159.      r = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "r");
  160.      phi_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "phi_sav");
  162. /*#if 0 */
  163.      i = sfreq + 0.5;
  164.      switch(i){
  165.         case 32000: sfreq_idx = 0; break;
  166.         case 44100: sfreq_idx = 1; break;
  167.         case 48000: sfreq_idx = 2; break;
  168.         default:    printf("error, invalid sampling frequency: %d Hzn",i);
  169.         exit(-1);
  170.      }
  171.      printf("absthr[][] sampling frequency index: %dn",sfreq_idx);
  172.      read_absthr(absthr, sfreq_idx);
  173.      switch(lay){
  174. case 1: sync_flush=576; flush=384; syncsize=1024; break;
  175. case 2: sync_flush=480; flush=576; syncsize=1056; break;
  176. case 3: sync_flush=768; flush=576; syncsize=1344; break;
  177.        default: printf("Bad lay value:(%d)",lay); exit(-1); break;
  178.      }
  179. /* #endif */
  181. /* calculate HANN window coefficients */
  182. /*   for(i=0;i<BLKSIZE;i++)  window[i]  =0.5*(1-cos(2.0*PI*i/(BLKSIZE-1.0)));*/
  183.      for(i=0;i<BLKSIZE;i++)  window[i]  =0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE));
  184.      for(i=0;i<BLKSIZE_s;i++)window_s[i]=0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE_s));
  185. /* reset states used in unpredictability measure */
  186.      for(i=0;i<HBLKSIZE;i++){
  187.         r[0][0][i]=r[1][0][i]=r[0][1][i]=r[1][1][i]=0;
  188.         phi_sav[0][0][i]=phi_sav[1][0][i]=0;
  189.         phi_sav[0][1][i]=phi_sav[1][1][i]=0;
  190.         lthr[0][i] = 60802371420160.0;
  191.         lthr[1][i] = 60802371420160.0;
  192.      }
  193. /*****************************************************************************
  194.  * Initialization: Compute the following constants for use later             *
  195.  *    partition[HBLKSIZE] = the partition number associated with each        *
  196.  *                          frequency line                                   *
  197.  *    cbval[CBANDS]       = the center (average) bark value of each          *
  198.  *                          partition                                        *
  199.  *    numlines[CBANDS]    = the number of frequency lines in each partition  *
  200.  *    tmn[CBANDS]         = tone masking noise                               *
  201.  *****************************************************************************/
  202. /* compute fft frequency multiplicand */
  203.      freq_mult = sfreq/BLKSIZE;
  204.  
  205. /* calculate fft frequency, then bval of each line (use fthr[] as tmp storage)*/
  206.      for(i=0;i<HBLKSIZE;i++){
  207.         temp1 = i*freq_mult;
  208.         j = 1;
  209.         while(temp1>crit_band[j])j++;
  210.         fthr[i]=j-1+(temp1-crit_band[j-1])/(crit_band[j]-crit_band[j-1]);
  211.      }
  212.      partition[0] = 0;
  213. /* temp2 is the counter of the number of frequency lines in each partition */
  214.      temp2 = 1;
  215.      cbval[0]=fthr[0];
  216.      bval_lo=fthr[0];
  217.      for(i=1;i<HBLKSIZE;i++){
  218.         if((fthr[i]-bval_lo)>0.33){
  219.            partition[i]=partition[i-1]+1;
  220.            cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  221.            cbval[partition[i]] = fthr[i];
  222.            bval_lo = fthr[i];
  223.            numlines[partition[i-1]] = temp2;
  224.            temp2 = 1;
  225.         }
  226.         else {
  227.            partition[i]=partition[i-1];
  228.            cbval[partition[i]] += fthr[i];
  229.            temp2++;
  230.         }
  231.      }
  232.      numlines[partition[i-1]] = temp2;
  233.      cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  234.  
  235. /************************************************************************
  236.  * Now compute the spreading function, s[j][i], the value of the spread-*
  237.  * ing function, centered at band j, for band i, store for later use    *
  238.  ************************************************************************/
  239.      for(j=0;j<CBANDS;j++){
  240.         for(i=0;i<CBANDS;i++){
  241.            temp1 = (cbval[i] - cbval[j])*1.05;
  242.            if(temp1>=0.5 && temp1<=2.5){
  243.               temp2 = temp1 - 0.5;
  244.               temp2 = 8.0 * (temp2*temp2 - 2.0 * temp2);
  245.            }
  246.            else temp2 = 0;
  247.            temp1 += 0.474;
  248.            temp3 = 15.811389+7.5*temp1-17.5*sqrt((double) (1.0+temp1*temp1));
  249.            if(temp3 <= -100) s[i][j] = 0;
  250.            else {
  251.               temp3 = (temp2 + temp3)*LN_TO_LOG10;
  252.               s[i][j] = exp(temp3);
  253.            }
  254.         }
  255.      }
  257.   /* Calculate Tone Masking Noise values */
  258.      for(j=0;j<CBANDS;j++){
  259.         temp1 = 15.5 + cbval[j];
  260.         tmn[j] = (temp1>24.5) ? temp1 : 24.5;
  261.   /* Calculate normalization factors for the net spreading functions */
  262.         rnorm[j] = 0;
  263.         for(i=0;i<CBANDS;i++){
  264.            rnorm[j] += s[j][i];
  265.         }
  266.      }
  267.      init++;
  268.  }
  269.  
  270. /************************* End of Initialization *****************************/
  271.  switch(lay) {
  272.   case 1:
  273.   case 2:
  274. for ( i=0; i<lay; i++)
  275.   {
  276. /*****************************************************************************
  277.  * Net offset is 480 samples (1056-576) for layer 2; this is because one must*
  278.  * stagger input data by 256 samples to synchronize psychoacoustic model with*
  279.  * filter bank outputs, then stagger so that center of 1024 FFT window lines *
  280.  * up with center of 576 "new" audio samples.                                *
  281.  *                                                                           *
  282.  * For layer 1, the input data still needs to be staggered by 256 samples,   *
  283.  * then it must be staggered again so that the 384 "new" samples are centered*
  284.  * in the 1024 FFT window.  The net offset is then 576 and you need 448 "new"*
  285.  * samples for each iteration to keep the 384 samples of interest centered   *
  286.  *****************************************************************************/
  287.   for (j=0; j<syncsize; j++)
  288.   {
  289.     if (j < (sync_flush) )
  290.       savebuf[j] = savebuf[j+flush];
  291.     else
  292.       savebuf[j] = *buffer++;
  294. /**window data with HANN window***********************************************/
  295.     if (j<BLKSIZE)
  296.     {
  297.       wsamp_r[j] = window[j]*((FLOAT) savebuf[j]); 
  298.       wsamp_i[j] = 0;
  299.     }
  300.   }
  301. /**Compute FFT****************************************************************/
  302.         fft(wsamp_r,wsamp_i,energy,phi,1024);
  303. /*****************************************************************************
  304.  * calculate the unpredictability measure, given energy[f] and phi[f]        *
  305.  *****************************************************************************/
  306.         for(j=0; j<HBLKSIZE; j++){
  307.            r_prime = 2.0 * r[chn][old][j] - r[chn][oldest][j];
  308.            phi_prime = 2.0 * phi_sav[chn][old][j] - phi_sav[chn][oldest][j];
  309.            r[chn][new][j] = sqrt((double) energy[j]);
  310.            phi_sav[chn][new][j] = phi[j];
  311.    temp1 = r[chn][new][j] * cos((double) phi[j])
  312.    - r_prime * cos(phi_prime);
  313.    temp2=r[chn][new][j] * sin((double) phi[j])
  314.    - r_prime * sin(phi_prime);
  315.            temp3=r[chn][new][j] + fabs(r_prime);
  316.            if(temp3 != 0)c[j]=sqrt(temp1*temp1+temp2*temp2)/temp3;
  317.            else c[j] = 0;
  318.         }
  319. /*only update data "age" pointers after you are done with the second channel */
  320. /*for layer 1 computations, for the layer 2 double computations, the pointers*/
  321. /*are reset automatically on the second pass                                 */
  322.         if(lay==2 || chn==1){
  323.            if(new==0){new = 1; oldest = 1;}
  324.            else {new = 0; oldest = 0;}
  325.            if(old==0)old = 1; else old = 0;
  326.         }
  327. /*****************************************************************************
  328.  * Calculate the grouped, energy-weighted, unpredictability measure,         *
  329.  * grouped_c[], and the grouped energy. grouped_e[]                          *
  330.  *****************************************************************************/
  331.         for(j=1;j<CBANDS;j++){
  332.            grouped_e[j] = 0;
  333.            grouped_c[j] = 0;
  334.         }
  335.         grouped_e[0] = energy[0];
  336.         grouped_c[0] = energy[0]*c[0];
  337.         for(j=1;j<HBLKSIZE;j++){
  338.            grouped_e[partition[j]] += energy[j];
  339.            grouped_c[partition[j]] += energy[j]*c[j];
  340.         }
  341. /*****************************************************************************
  342.  * convolve the grouped energy-weighted unpredictability measure             *
  343.  * and the grouped energy with the spreading function, s[j][k]               *
  344.  *****************************************************************************/
  345.         for(j=0;j<CBANDS;j++){
  346.            ecb[j] = 0;
  347.            cb[j] = 0;
  348.            for(k=0;k<CBANDS;k++){
  349.               if(s[j][k] != 0.0){
  350.                  ecb[j] += s[j][k]*grouped_e[k];
  351.                  cb[j] += s[j][k]*grouped_c[k];
  352.               }
  353.            }
  354.            if(ecb[j] !=0)cb[j] = cb[j]/ecb[j];
  355.            else cb[j] = 0;
  356.         }
  357. /*****************************************************************************
  358.  * Calculate the required SNR for each of the frequency partitions           *
  359.  *         this whole section can be accomplished by a table lookup          *
  360.  *****************************************************************************/
  361.         for(j=0;j<CBANDS;j++){
  362.            if(cb[j]<.05)cb[j]=0.05;
  363.            else if(cb[j]>.5)cb[j]=0.5;
  364. /*         tb = -0.434294482*log((double) cb[j])-0.301029996; */
  365.            tb = -0.43 *log((double) cb[j]) - 0.29 ;
  366.            if(tb<0.0) tb=0.0; else if(tb>1.0) tb=1.0;
  367.            bc[j] = tmn[j]*tb + nmt*(1.0-tb);
  368.            k = cbval[j] + 0.5;
  369.            bc[j] = (bc[j] > bmax[k]) ? bc[j] : bmax[k];
  370.            bc[j] = exp((double) -bc[j]*LN_TO_LOG10);
  371.         }
  372. /*****************************************************************************
  373.  * Calculate the permissible noise energy level in each of the frequency     *
  374.  * partitions. Include absolute threshold and pre-echo controls              *
  375.  *         this whole section can be accomplished by a table lookup          *
  376.  *****************************************************************************/
  377.         for(j=0;j<CBANDS;j++)
  378.            if(rnorm[j] && numlines[j])
  379.               nb[j] = ecb[j]*bc[j]/(rnorm[j]*numlines[j]);
  380.            else nb[j] = 0;
  381.         for(j=0;j<HBLKSIZE;j++){
  382.            temp1=nb[partition[j]];  /* preliminary threshold */
  383.            temp1=(temp1>absthr[j])?temp1:absthr[j];
  384. /*do not use pre-echo control for layer 2 because it may do bad things to the*/
  385. /*  MUSICAM bit allocation algorithm                                         */
  386.            if(lay==1){
  387.               fthr[j] = (temp1 < lthr[chn][j]) ? temp1 : lthr[chn][j];
  388.               temp2 = temp1 * 0.00316;
  389.               fthr[j] = (temp2 > fthr[j]) ? temp2 : fthr[j];
  390.            }
  391.            else fthr[j] = temp1;
  392.            lthr[chn][j] = LXMIN*temp1;
  393.         }
  394. /*****************************************************************************
  395.  * Translate the 512 threshold values to the 32 filter bands of the coder    *
  396.  *****************************************************************************/
  397.         for(j=0;j<193;j += 16){
  398.            min_thres = fthr[j];
  399.            sum_energy = energy[j];
  400.            for(k=1;k<17;k++){
  401.               if(min_thres>fthr[j+k]) min_thres = fthr[j+k];
  402.               sum_energy += energy[j+k];
  403.            }
  404.            snrtmp[i][j/16] = sum_energy/(min_thres * 17.0);
  405.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  406.         }
  407.         for(j=208;j<(HBLKSIZE-1);j += 16){
  408.            min_thres = 0.0;
  409.            sum_energy = 0.0;
  410.            for(k=0;k<17;k++){
  411.               min_thres += fthr[j+k];
  412.               sum_energy += energy[j+k];
  413.            }
  414.            snrtmp[i][j/16] = sum_energy/min_thres;
  415.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  416.         }
  417. /*****************************************************************************
  418.  * End of Psychoacuostic calculation loop                                    *
  419.  *****************************************************************************/
  420.      }
  421.      for(i=0; i<32; i++){
  422.         if(lay==2) /* if(lay==2 && chn==2) MI */
  423.            snr32[i]=(snrtmp[0][i]>snrtmp[1][i])?snrtmp[0][i]:snrtmp[1][i];
  424.         else snr32[i]=snrtmp[0][i];
  425.      }
  426.      break;
  428. /*************************************************************************/
  429. /**       LAYER 3                                                        */
  430. /*************************************************************************/
  432.    case 3:
  433. if ( init_L3 == 0 )
  434. {
  435.     L3para_read( sfreq,numlines,partition_l,minval,qthr_l,norm_l,s3_l,
  436.  partition_s,qthr_s,norm_s,SNR_s,
  437.  cbw_l,bu_l,bo_l,w1_l,w2_l, cbw_s,bu_s,bo_s,w1_s,w2_s );
  438.     init_L3 ++ ;
  439. }
  441. for ( j = 0; j < 21; j++ )
  442.     ratio_d[j] = ratio[chn][j];
  443. for ( j = 0; j < 12; j++ )
  444.     for ( i = 0; i < 3; i++ )
  445. ratio_ds[j][i] = ratio_s[chn][j][i];
  447. if ( chn == 0 )
  448.     if ( new == 0 )
  449.     {
  450. new = 1;
  451. old = 0;
  452. oldest = 1;
  453.     }
  454.     else
  455.     {
  456. new = 0;
  457. old = 1;
  458. oldest = 0;
  459.     }
  462. /**********************************************************************
  463. *  Delay signal by sync_flush=768 samples                             *
  464. **********************************************************************/
  466. for ( j = 0; j < sync_flush; j++ ) /* for long window samples */
  467.     savebuf[j] = savebuf[j+flush];
  469. for ( j = sync_flush; j < syncsize; j++ )
  470.     savebuf[j] = *buffer++;
  472. for ( j = 0; j < BLKSIZE; j++ )
  473. { /**window data with HANN window**/
  474.     wsamp_r[j] = window[j] * savebuf[j];  
  475.     wsamp_i[j] = 0.0;
  476. }
  479. /**********************************************************************
  480. *    compute unpredicatability of first six spectral lines            * 
  481. **********************************************************************/
  483. fft( wsamp_r, wsamp_i, energy, phi, 1024 ); /**long FFT**/
  484. for ( j = 0; j < 6; j++ )
  485. {  /* calculate unpredictability measure cw */
  486.     r_prime = 2.0 * r[chn][old][j] - r[chn][oldest][j];
  487.     phi_prime = 2.0 * phi_sav[chn][old][j]-phi_sav[chn][oldest][j];
  488.     r[chn][new][j] = sqrt((double) energy[j]);
  489.     phi_sav[chn][new][j] = phi[j];
  490.     temp1 = r[chn][new][j] * cos((double) phi[j])
  491. - r_prime * cos(phi_prime);
  492.     temp2 = r[chn][new][j] * sin((double) phi[j])
  493. - r_prime * sin(phi_prime);
  494.     temp3=r[chn][new][j] + fabs(r_prime);
  495.     
  496.     if ( temp3 != 0.0 )
  497. cw[j] = sqrt( temp1*temp1+temp2*temp2 ) / temp3;
  498.     else
  499. cw[j] = 0;
  500. }
  503. /**********************************************************************
  504. *     compute unpredicatibility of next 200 spectral lines            *
  505. **********************************************************************/ 
  506. for ( sblock = 0; sblock < 3; sblock++ )
  507. { /**window data with HANN window**/
  508.     for ( j = 0, k = 128 * (2 + sblock); j < 256; j++, k++ )
  509.     {
  510. wsamp_r[j] = window_s[j]* savebuf[k]; 
  511. wsamp_i[j] = 0.0;
  512.     } /* short FFT*/
  513.     
  514.     fft( wsamp_r, wsamp_i, &energy_s[sblock][0], &phi_s[sblock][0], 256 );
  515.         }
  516.  
  517.         sblock = 1;
  519. for ( j = 6; j < 206; j += 4 )
  520. {/* calculate unpredictability measure cw */
  521.     double r2, phi2, temp1, temp2, temp3;
  522.     
  523.     k = (j+2) / 4; 
  524.     r_prime = 2.0 * sqrt((double) energy_s[0][k])
  525. - sqrt((double) energy_s[2][k]);
  526.     phi_prime = 2.0 * phi_s[0][k] - phi_s[2][k];
  527.     r2 = sqrt((double) energy_s[1][k]);
  528.     phi2 = phi_s[1][k];
  529.     temp1 = r2 * cos( phi2 ) - r_prime * cos( phi_prime );
  530.     temp2 = r2 * sin( phi2 ) - r_prime * sin( phi_prime );
  531.     temp3 = r2 + fabs( r_prime );
  532.     if ( temp3 != 0.0 )
  533. cw[j] = sqrt( temp1 * temp1 + temp2 * temp2 ) / temp3;
  534.     else
  535. cw[j] = 0.0;
  536.     cw[j+1] = cw[j+2] = cw[j+3] = cw[j];
  537. }
  540. /**********************************************************************
  541. *    Set unpredicatiblility of remaining spectral lines to 0.4        *
  542. **********************************************************************/
  543. for ( j = 206; j < HBLKSIZE; j++ )
  544.     cw[j] = 0.4;
  548. /**********************************************************************
  549. *    Calculate the energy and the unpredictability in the threshold   *
  550. *    calculation partitions                                           *
  551. **********************************************************************/
  553. for ( b = 0; b < CBANDS; b++ )
  554. {
  555.     eb[b] = 0.0;
  556.     cb[b] = 0.0;
  557. }
  558. for ( j = 0; j < HBLKSIZE; j++ )
  559. {
  560.     int tp = partition_l[j];
  561.     if ( tp >= 0 )
  562.     {
  563. eb[tp] += energy[j];
  564. cb[tp] += cw[j] * energy[j];
  565.     }
  566. }
  569. /**********************************************************************
  570. *      convolve the partitioned energy and unpredictability           *
  571. *      with the spreading function, s3_l[b][k]                        *
  572. ******************************************************************** */
  574. for ( b = 0; b < CBANDS; b++ )
  575. {
  576.     ecb[b] = 0.0;
  577.     ctb[b] = 0.0;
  578. }
  579. for ( b = 0;b < CBANDS; b++ )
  580. {
  581.     for ( k = 0; k < CBANDS; k++ )
  582.     {
  583. ecb[b] += s3_l[b][k] * eb[k]; /* sprdngf for Layer III */
  584. ctb[b] += s3_l[b][k] * cb[k];
  585.     }
  586. }
  588. /* calculate the tonality of each threshold calculation partition */
  589. /* calculate the SNR in each threshhold calculation partition */
  591. for ( b = 0; b < CBANDS; b++ )
  592. {
  593.     double cbb,tbb;
  594.     if (ecb[b] != 0.0 )
  595.                 {
  596. cbb = ctb[b]/ecb[b];
  597.                 if (cbb <0.01) cbb = 0.01;
  598. cbb = log( cbb);
  599.                 }
  600.     else
  601. cbb = 0.0 ;
  602.     tbb = -0.299 - 0.43*cbb;  /* conv1=-0.299, conv2=-0.43 */
  603.     tbb = minimum( 1.0, maximum( 0.0, tbb) ) ;  /* 0<tbb<1 */
  604.     SNR_l[b] = maximum( minval[b], 29.0*tbb+6.0*(1.0-tbb) );
  605. } /* TMN=29.0,NMT=6.0 for all calculation partitions */
  607. for ( b = 0; b < CBANDS; b++ ) /* calculate the threshold for each partition */
  608.     nb[b] = ecb[b] * norm_l[b] * exp( -SNR_l[b] * LN_TO_LOG10 );
  610. for ( b = 0; b < CBANDS; b++ )
  611. { /* pre-echo control */
  612.     double temp_1; /* BUG of IS */
  613.     temp_1 = minimum( nb[b], minimum(2.0*nb_1[chn][b],16.0*nb_2[chn][b]) );
  614.     thr[b] = maximum( qthr_l[b], temp_1 );/* rpelev=2.0, rpelev2=16.0 */
  615.     nb_2[chn][b] = nb_1[chn][b];
  616.     nb_1[chn][b] = nb[b];
  617. }
  620. *pe = 0.0; /*  calculate percetual entropy */
  621. for ( b = 0; b < CBANDS; b++ )
  622. {
  623.     double tp ;
  624.     tp = minimum( 0.0, log((thr[b]+1.0) / (eb[b]+1.0) ) ) ; /*not log*/
  625.     *pe -= numlines[b] * tp ;
  626. } /* thr[b] -> thr[b]+1.0 : for non sound portition */
  628. #define switch_pe  1800
  629.         blocktype = NORM_TYPE;
  632. if ( *pe < switch_pe )
  633. { /* no attack : use long blocks */
  634.     switch( blocktype_old[chn] ) 
  635.     {
  636.       case NORM_TYPE:
  637.       case STOP_TYPE:
  638. blocktype = NORM_TYPE;
  639. break;
  640.     
  641.       case SHORT_TYPE:
  642. blocktype = STOP_TYPE;
  643. break;
  644.     
  645.       case START_TYPE:
  646. fprintf( stderr, "Error in block selectingn" );
  647. abort();
  648. break; /* problem */
  649.     }
  651.     /* threshold calculation (part 2) */
  652.     for ( sb = 0; sb < SBMAX_l; sb++ )
  653.     {
  654. en[sb] = w1_l[sb] * eb[bu_l[sb]] + w2_l[sb] * eb[bo_l[sb]];
  655. thm[sb] = w1_l[sb] *thr[bu_l[sb]] + w2_l[sb] * thr[bo_l[sb]];
  656. for ( b = bu_l[sb]+1; b < bo_l[sb]; b++ )
  657. {
  658.     en[sb]  += eb[b];
  659.     thm[sb] += thr[b];
  660. }
  661. if ( en[sb] != 0.0 )
  662.     ratio[chn][sb] = thm[sb]/en[sb];
  663. else
  664.     ratio[chn][sb] = 0.0;
  665.     }
  666. }
  667. else 
  668. {
  669.     /* attack : use short blocks */
  670.     blocktype = SHORT_TYPE;
  671.     
  672.     if ( blocktype_old[chn] == NORM_TYPE ) 
  673. blocktype_old[chn] = START_TYPE;
  674.     if ( blocktype_old[chn] == STOP_TYPE )
  675. blocktype_old[chn] = SHORT_TYPE ;
  676.     
  677.     /* threshold calculation for short blocks */
  678.     
  679.     for ( sblock = 0; sblock < 3; sblock++ )
  680.     {
  681. for ( b = 0; b < CBANDS_s; b++ )
  682. {
  683.     eb[b] = 0.0;
  684.     ecb[b] = 0.0;
  685. }
  686. for ( j = 0; j < HBLKSIZE_s; j++ )
  687.     eb[partition_s[j]] += energy_s[sblock][j];
  688. for ( b = 0; b < CBANDS_s; b++ )
  689.     for ( k = 0; k < CBANDS_s; k++ )
  690. ecb[b] += s3_l[b][k] * eb[k];
  691. for ( b = 0; b < CBANDS_s; b++ )
  692. {
  693.     nb[b] = ecb[b] * norm_l[b] * exp( (double) SNR_s[b] * LN_TO_LOG10 );
  694.     thr[b] = maximum (qthr_s[b],nb[b]);
  695. }
  696. for ( sb = 0; sb < SBMAX_s; sb++ )
  697. {
  698.     en[sb] = w1_s[sb] * eb[bu_s[sb]] + w2_s[sb] * eb[bo_s[sb]];
  699.     thm[sb] = w1_s[sb] *thr[bu_s[sb]] + w2_s[sb] * thr[bo_s[sb]];
  700.     for ( b = bu_s[sb]+1; b < bo_s[sb]; b++ )
  701.     {
  702. en[sb] += eb[b];
  703. thm[sb] += thr[b];
  704.     }
  705.     if ( en[sb] != 0.0 )
  706. ratio_s[chn][sb][sblock] = thm[sb]/en[sb];
  707.     else
  708. ratio_s[chn][sb][sblock] = 0.0;
  709. }
  710.     }
  711. } 
  713. cod_info->block_type = blocktype_old[chn];
  714. blocktype_old[chn] = blocktype;
  716. if ( cod_info->block_type == NORM_TYPE )
  717.     cod_info->window_switching_flag = 0;
  718. else
  719.     cod_info->window_switching_flag = 1;
  720. cod_info->mixed_block_flag = 0;
  721. break;
  724.   default:
  725.      printf("error, invalid MPEG/audio coding layer: %dn",lay);
  726.  }
  728. /* These mem_free() calls must correspond with the mem_alloc() calls     */
  729. /* used at the beginning of this function to simulate "automatic"        */
  730. /* variables placed on the stack.                                        */
  732.  mem_free((void **) &grouped_c);
  733.  mem_free((void **) &grouped_e);
  734.  mem_free((void **) &nb);
  735.  mem_free((void **) &cb);
  736.  mem_free((void **) &ecb);
  737.  mem_free((void **) &bc);
  738.  mem_free((void **) &wsamp_r);
  739.  mem_free((void **) &wsamp_i);
  740.  mem_free((void **) &phi);
  741.  mem_free((void **) &energy);
  742.  mem_free((void **) &c);
  743.  mem_free((void **) &fthr);
  744.  mem_free((void **) &snrtmp);
  745. }
  746. #ifdef DEBUG
  747. #undef DEBUG
  748. #endif
  751. void L3para_read(double sfreq, int *numlines, int *partition_l, double *minval, double *qthr_l, double *norm_l, double (*s3_l)[63], int *partition_s, double *qthr_s, double *norm_s, double *SNR, int *cbw_l, int *bu_l, int *bo_l, double *w1_l, double *w2_l, int *cbw_s, int *bu_s, int *bo_s, double *w1_s, double *w2_s)
  752. {
  753.    double freq_tp;
  754.    static double bval_l[CBANDS], bval_s[CBANDS];
  755.    int   cbmax, cbmax_tp;
  756.    static double s3_s[CBANDS][CBANDS];
  758.    FILE *fin;
  759.    char tp[256];
  760.    int  sbmax ;
  761.    int  i,j,k,k2,loop, part_max ;
  763.    fin = OpenTableFile( "psy_data" );
  764.    if (fin == NULL)
  765.        exit( 1 );
  767. /* Read long block data */
  769.       for(loop=0;loop<6;loop++)
  770.       {
  771. fscanf(fin,"freq = %lf partition = %dn",&freq_tp,&cbmax_tp);
  772. cbmax_tp++;
  773. #ifdef DEBUG
  774. printf("freq = %f partition = %dn",freq_tp,cbmax);
  775. #endif
  777. if (sfreq == freq_tp )
  778.   {
  779.      cbmax = cbmax_tp;
  780.      for(i=0,k2=0;i<cbmax_tp;i++)
  781.        {
  782. fscanf(fin,
  783.   "No=%d #lines=%d minval=%lf qthr=%lf norm=%lf bval=%lfn",
  784.   &j,&numlines[i],&minval[i],&qthr_l[i],&norm_l[i],&bval_l[i]);
  785.         if (j!=i)
  786.          { printf("please check "psy_data"");
  787.    exit(-1);
  788.          }
  789. for(k=0;k<numlines[i];k++)
  790.   partition_l[k2++] = i ;
  791. #ifdef DEBUG
  792.      printf("No=%2d #lines=%2d minval=%4.1f qthr=%8.3f norm=%5.3f bval=%8.3fn",
  793.      i,numlines[i],minval[i],qthr_l[i],norm_l[i],bval_l[i] );
  794. #endif
  795. }
  796.    }
  797.    else
  798.    {
  799.      for(j=0;j<cbmax_tp;j++)
  800.        {
  801. fgets(tp,255,fin);
  802.         sscanf(tp,"No=%d %sn",&i,tp);
  803.         if (j!=i)
  804.          { printf("please check "psy_data."n");
  805.    exit(-1);
  806.          }
  807.        }
  808.    }
  809.        }
  811. /************************************************************************
  812.  * Now compute the spreading function, s[j][i], the value of the spread-*
  813.  * ing function, centered at band j, for band i, store for later use    *
  814.  ************************************************************************/
  815. #ifdef DEBUG
  816. printf("freq = %fn",sfreq);
  817. #endif
  818.   part_max = cbmax ;
  819.           for(i=0;i<part_max;i++)
  820.   {
  821.   double tempx,x,tempy,temp;
  822.             for(j=0;j<part_max;j++)
  823.     {
  824.              tempx = (bval_l[i] - bval_l[j])*1.05;
  825.              if (j>=i) tempx = (bval_l[i] - bval_l[j])*3.0;
  826.                else    tempx = (bval_l[i] - bval_l[j])*1.5;
  827. /*             if (j>=i) tempx = (bval_l[j] - bval_l[i])*3.0;
  828.                else    tempx = (bval_l[j] - bval_l[i])*1.5; */
  829.              if(tempx>=0.5 && tempx<=2.5)
  830.      {
  831.                temp = tempx - 0.5;
  832.                x = 8.0 * (temp*temp - 2.0 * temp);
  833.              }
  834.              else x = 0.0;
  835.              tempx += 0.474;
  836.              tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  837.              if (tempy <= -60.0) s3_l[i][j] = 0.0;
  838.              else                s3_l[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
  839. #ifdef DEBUG_S3
  840.      printf("s[%d][%d]=%fn",i,j,s3_l[i][j]);
  841. #endif
  843. #ifdef DEBUGP
  844.      printf("j=%d i=%d tempy=%f s[i][j]=%f n",i,j,tempy,s[i][j]);
  845.    minval[i] = bmax[j-1];
  846.    printf("minval[%d] = %f, j-1=%d %fn",i, minval[i] , j,fthr[i]) ;
  847. #endif
  848.             }
  849.           }
  852. /* Read short block data */
  854.       for(loop=0;loop<6;loop++)
  855.       {
  856. fscanf(fin,"freq = %lf partition = %dn",&freq_tp,&cbmax_tp);
  857. cbmax_tp++;
  858. #ifdef DEBUG
  859. printf("freq = %f partition = %dn",freq_tp,cbmax);
  860. #endif
  862. if (sfreq == freq_tp )
  863.   {
  864.      cbmax = cbmax_tp;
  865.      for(i=0,k2=0;i<cbmax_tp;i++)
  866.        {
  867. fscanf(fin,
  868.   "No=%d #lines=%d qthr=%lf norm=%lf SNR=%lf bval=%lfn",
  869.    &j,&numlines[i],&qthr_s[i],&norm_s[i],&SNR[i],&bval_s[i]);
  870.         if (j!=i)
  871.          { printf("please check "psy_data"");
  872.    exit(-1);
  873.          }
  874. for(k=0;k<numlines[i];k++)
  875.   partition_s[k2++] = i ;
  876. #ifdef DEBUG
  877.       printf("No=%2d #lines=%2d qthr=%8.3f norm=%5.3f SNR=%6.3f bval=%8.3fn",
  878.       i,numlines[i],qthr_s[i],norm_s[i],SNR[i],bval_s[i] );
  879. #endif
  880. }
  881.    }
  882.    else
  883.    {
  884.      for(j=0;j<cbmax_tp;j++)
  885.        {
  886. fgets(tp,255,fin);
  887.         sscanf(tp,"No=%d %sn",&i,tp);
  888.         if (j!=i)
  889.          { printf("please check "psy_data."n");
  890.    exit(-1);
  891.          }
  892.        }
  893.    }
  894.        }
  896. /************************************************************************
  897.  * Now compute the spreading function, s[j][i], the value of the spread-*
  898.  * ing function, centered at band j, for band i, store for later use    *
  899.  ************************************************************************/
  900. #ifdef DEBUG_S3
  901. fpp=fopen("s3_s","w");
  902. #endif
  903.   part_max = cbmax ;
  904.           for(i=0;i<part_max;i++)
  905.   {
  906.   double tempx,x,tempy,temp;
  907.             for(j=0;j<part_max;j++)
  908.     {
  909.              tempx = (bval_s[i] - bval_s[j])*1.05;
  910.              if (j>=i) tempx = (bval_s[i] - bval_s[j])*3.0;
  911.                else    tempx = (bval_s[i] - bval_s[j])*1.5;
  912.              if(tempx>=0.5 && tempx<=2.5)
  913.      {
  914.                temp = tempx - 0.5;
  915.                x = 8.0 * (temp*temp - 2.0 * temp);
  916.              }
  917.              else x = 0.0;
  918.              tempx += 0.474;
  919.              tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  920.              if (tempy <= -60.0) s3_s[i][j] = 0.0;
  921.              else                s3_s[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
  922. #ifdef DEBUG_S3
  923.      fprintf(fpp,"s3_s[%d][%d]=%fn",i,j,s3_s[i][j]);
  924. #endif
  925. #ifdef DEBUGP
  926.      printf("j=%d i=%d tempy=%f s[i][j]=%f n",i,j,tempy,s[i][j]);
  927.    minval[i] = bmax[j-1];
  928.    printf("minval[%d] = %f, j-1=%d %fn",i, minval[i] , j,fthr[i]) ;
  929. #endif
  930.             }
  931.           }
  932. #ifdef DEBUG_S3
  933. fclose(fpp);
  934. #endif
  935. /* Read long block data for converting threshold calculation 
  936.    partitions to scale factor bands */
  938.       for(loop=0;loop<6;loop++)
  939.       {
  940. fscanf(fin,"freq=%lf sb=%dn",&freq_tp,&sbmax);
  941. sbmax++;
  942. #ifdef DEBUG
  943. printf("freq = %f sb = %dn",freq_tp,sbmax);
  944. #endif
  946. if (sfreq == freq_tp)
  947.   {
  948.      for(i=0;i<sbmax;i++)
  949.       {
  950. fscanf(fin,
  951.   "sb=%d cbw=%d bu=%d bo=%d w1=%lf w2=%lfn",
  952.   &j,&cbw_l[i],&bu_l[i],&bo_l[i],&w1_l[i],&w2_l[i]);
  953.         if (j!=i)
  954.          { printf("30:please check "psy_data"n");
  955.    exit(-1);
  956.          }
  957. #ifdef DEBUG
  958. printf(
  959.   "sb=%2d cbw=%1d bu=%2d bo=%2d w1=%5.3f w2=%5.3fn",
  960.   j,cbw_l[i],bu_l[i],bo_l[i],w1_l[i],w2_l[i]);
  961. #endif
  962.         if (i!=0)
  963.  if ( (bo_l[i] != (bu_l[i]+cbw_l[i])) ||
  964.  (fabs(1.0-w1_l[i]-w2_l[i-1]) > 0.01 ) )
  965.          { printf("31:please check "psy_data."n");
  966.    exit(-1);
  967.          }
  968.       }
  969.    }
  970.    else
  971.    {
  972.      for(j=0;j<sbmax;j++)
  973.        {
  974. fgets(tp,255,fin);
  975.         sscanf(tp,"sb=%d %sn",&i,tp);
  976.         if (j!=i)
  977.          { printf("please check "psy_data."n");
  978.    exit(-1);
  979.          }
  980.        }
  981.    }
  982.        }
  984. /* Read short block data for converting threshold calculation 
  985.    partitions to scale factor bands */
  987.       for(loop=0;loop<6;loop++)
  988.       {
  989. fscanf(fin,"freq=%lf sb=%dn",&freq_tp,&sbmax);
  990. sbmax++;
  991. #ifdef DEBUG
  992. printf("freq = %f sb = %dn",freq_tp,sbmax);
  993. #endif
  995. if (sfreq == freq_tp)
  996.   {
  997.      for(i=0;i<sbmax;i++)
  998.       {
  999. fscanf(fin,
  1000.   "sb=%d cbw=%d bu=%d bo=%d w1=%lf w2=%lfn",
  1001.   &j,&cbw_s[i],&bu_s[i],&bo_s[i],&w1_s[i],&w2_s[i]);
  1002.         if (j!=i)
  1003.          { printf("30:please check "psy_data"n");
  1004.    exit(-1);
  1005.          }
  1006. #ifdef DEBUG
  1007. printf(
  1008.   "sb=%2d cbw=%1d bu=%2d bo=%2d w1=%5.3f w2=%5.3fn",
  1009.   j,cbw_s[i],bu_s[i],bo_s[i],w1_s[i],w2_s[i]);
  1010. #endif
  1011.         if (i!=0)
  1012.  if ( (bo_s[i] != (bu_s[i]+cbw_s[i])) ||
  1013.  (fabs(1.0-w1_s[i]-w2_s[i-1]) > 0.01 ) )
  1014.          { printf("31:please check "psy_data."n");
  1015.    exit(-1);
  1016.          }
  1017.       }
  1018.    }
  1019.    else
  1020.    {
  1021.      for(j=0;j<sbmax;j++)
  1022.        {
  1023. fgets(tp,255,fin);
  1024.         sscanf(tp,"sb=%d %sn",&i,tp);
  1025.         if (j!=i)
  1026.          { printf("please check "psy_data."n");
  1027.    exit(-1);
  1028.          }
  1029.        }
  1030.    }
  1031.        }
  1033. }