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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > 流媒体/Mpeg4/MP4 > 查看源码
psymodel.c
资源名称：NETVIDEO.rar　[点击查看]
上传用户：sun1608
上传日期：2007-02-02
资源大小：6116k
文件大小：33k
源码类别：

流媒体/Mpeg4/MP4

开发平台：

Visual C++

psymodel.c：源码内容
  1. /**********************************************************************
  2.  *   date   programmers         comment                               *
  3.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  4.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  5.  * 7/10/91  Earle Jennings      Ported to MsDos.                      *
  6.  *                              replace of floats with FLOAT          *
  7.  * 2/11/92  W. Joseph Carter    Fixed mem_alloc() arg for "absthr".   *
  8.  * 3/16/92  Masahiro Iwadare Modification for Layer III            *
  9.  * 17/4/93  Masahiro Iwadare    Updated for IS Modification           *
  10.  **********************************************************************/
  12. #include "util.h"
  13. #include "encoder.h"
  14. #include "psymodel.h"
  15. #include "l3side.h"
  16. #include <assert.h>
  17. #ifdef HAVEGTK
  18. #include "gtkanal.h"
  19. #endif
  20. #include "tables.h"
  21. #include "fft.h"
  23. #ifdef M_LN10
  24. #define LN_TO_LOG10 (M_LN10/10)
  25. #else
  26. #define         LN_TO_LOG10             0.2302585093
  27. #endif
  30. void L3para_read( FLOAT8 sfreq, int numlines[CBANDS],int numlines_s[CBANDS], int partition_l[HBLKSIZE],
  31.   FLOAT8 minval[CBANDS], FLOAT8 qthr_l[CBANDS], 
  32.   FLOAT8 s3_l[CBANDS + 1][CBANDS + 1],
  33.   FLOAT8 s3_s[CBANDS + 1][CBANDS + 1],
  34.                   FLOAT8 qthr_s[CBANDS],
  35.   FLOAT8 SNR_s[CBANDS],
  36.   int bu_l[SBPSY_l], int bo_l[SBPSY_l],
  37.   FLOAT8 w1_l[SBPSY_l], FLOAT8 w2_l[SBPSY_l],
  38.   int bu_s[SBPSY_s], int bo_s[SBPSY_s],
  39.   FLOAT8 w1_s[SBPSY_s], FLOAT8 w2_s[SBPSY_s] );
  48.  
  50. void L3psycho_anal( lame_global_flags *gfp,
  51.                     short int *buffer[2],int gr_out , 
  52.                     FLOAT8 *ms_ratio,
  53.                     FLOAT8 *ms_ratio_next,
  54.     FLOAT8 *ms_ener_ratio,
  55.     III_psy_ratio masking_ratio[2][2],
  56.     III_psy_ratio masking_MS_ratio[2][2],
  57.     FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2], 
  58.                     int blocktype_d[2])
  59. {
  61. /* to get a good cache performance, one has to think about
  62.  * the sequence, in which the variables are used
  63.  */
  64.   
  65. /* The static variables "r", "phi_sav", "new", "old" and "oldest" have    */
  66. /* to be remembered for the unpredictability measure.  For "r" and        */
  67. /* "phi_sav", the first index from the left is the channel select and     */
  68. /* the second index is the "age" of the data.                             */
  69.   static FLOAT8 minval[CBANDS],qthr_l[CBANDS];
  70.   static FLOAT8 qthr_s[CBANDS];
  71.   static FLOAT8 nb_1[4][CBANDS], nb_2[4][CBANDS];
  72.   static FLOAT8 s3_s[CBANDS + 1][CBANDS + 1];
  73.   static FLOAT8 s3_l[CBANDS + 1][CBANDS + 1];
  75.   static III_psy_xmin thm[4];
  76.   static III_psy_xmin en[4];
  77.   
  78.   /* unpredictability calculation
  79.    */
  80.   static int cw_upper_index;
  81.   static int cw_lower_index;
  82.   static FLOAT ax_sav[4][2][HBLKSIZE];
  83.   static FLOAT bx_sav[4][2][HBLKSIZE];
  84.   static FLOAT rx_sav[4][2][HBLKSIZE];
  85.   static FLOAT cw[HBLKSIZE];
  87.   /* fft and energy calculation
  88.    */
  89.   FLOAT (*wsamp_l)[BLKSIZE];
  90.   FLOAT (*wsamp_s)[3][BLKSIZE_s];
  91.   FLOAT tot_ener[4];
  92.   static FLOAT wsamp_L[2][BLKSIZE];
  93.   static FLOAT energy[HBLKSIZE];
  94.   static FLOAT wsamp_S[2][3][BLKSIZE_s];
  95.   static FLOAT energy_s[3][HBLKSIZE_s];
  97.   /* convolution
  98.    */
  99.   static FLOAT8 eb[CBANDS];
  100.   static FLOAT8 cb[CBANDS];
  101.   static FLOAT8 thr[CBANDS];
  102.   
  103.   /* Scale Factor Bands
  104.    */
  105.   static FLOAT8 w1_l[SBPSY_l], w2_l[SBPSY_l];
  106.   static FLOAT8 w1_s[SBPSY_s], w2_s[SBPSY_s];
  107.   static FLOAT8 mld_l[SBPSY_l],mld_s[SBPSY_s];
  108.   static int bu_l[SBPSY_l],bo_l[SBPSY_l] ;
  109.   static int bu_s[SBPSY_s],bo_s[SBPSY_s] ;
  110.   static int npart_l,npart_s;
  111.   static int npart_l_orig,npart_s_orig;
  112.   
  113.   static int s3ind[CBANDS][2];
  114.   static int s3ind_s[CBANDS][2];
  116.   static int numlines_s[CBANDS] ;
  117.   static int numlines_l[CBANDS];
  118.   static int partition_l[HBLKSIZE];
  119.   
  120.   /* frame analyzer 
  121.    */
  122. #ifdef HAVEGTK
  123.   static FLOAT energy_save[4][HBLKSIZE];
  124.   static FLOAT8 pe_save[4];
  125.   static FLOAT8 ers_save[4];
  126. #endif
  128.   /* ratios 
  129.    */
  130.   static FLOAT8 pe[4]={0,0,0,0};
  131.   static FLOAT8 ms_ratio_s_old=0,ms_ratio_l_old=0;
  132.   static FLOAT8 ms_ener_ratio_old=.25;
  133.   FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
  135.   /* block type 
  136.    */
  137.   static int blocktype_old[2];
  138.   int blocktype[2],uselongblock[2];
  139.   
  140.   /* usual variables like loop indices, etc..
  141.    */
  142.   int numchn, chn;
  143.   int   b, i, j, k;
  144.   int sb,sblock;
  145.   FLOAT cwlimit;
  148.   /* initialization of static variables
  149.    */
  150.   if((gfp->frameNum==0) && (gr_out==0)){
  151.     FLOAT8 SNR_s[CBANDS];
  152.     
  153.     blocktype_old[0]=STOP_TYPE;
  154.     blocktype_old[1]=STOP_TYPE;
  155.     i = gfp->out_samplerate;
  156.     switch(i){
  157.     case 32000: break;
  158.     case 44100: break;
  159.     case 48000: break;
  160.     case 16000: break;
  161.     case 22050: break;
  162.     case 24000: break;
  163.     default:    fprintf(stderr,"error, invalid sampling frequency: %d Hzn",i);
  164.       exit(-1);
  165.     }
  166.     
  167.     /* reset states used in unpredictability measure */
  168.     memset (rx_sav,0, sizeof(rx_sav));
  169.     memset (ax_sav,0, sizeof(ax_sav));
  170.     memset (bx_sav,0, sizeof(bx_sav));
  171.     memset (en,0, sizeof(en));
  172.     memset (thm,0, sizeof(thm));
  173.     
  175.     /*  gfp->cwlimit = sfreq*j/1024.0;  */
  176.     cw_lower_index=6;
  177.     if (gfp->cwlimit>0) 
  178.       cwlimit=gfp->cwlimit;
  179.     else
  180.       cwlimit=8.8717;
  181.     cw_upper_index = cwlimit*1000.0*1024.0/((FLOAT8) gfp->out_samplerate);
  182.     cw_upper_index=Min(HBLKSIZE-4,cw_upper_index);      /* j+3 < HBLKSIZE-1 */
  183.     cw_upper_index=Max(6,cw_upper_index);
  185.     for ( j = 0; j < HBLKSIZE; j++ )
  186.       cw[j] = 0.4;
  187.     
  188.     /* setup stereo demasking thresholds */
  189.     /* formula reverse enginerred from plot in paper */
  190.     for ( sb = 0; sb < SBPSY_s; sb++ ) {
  191.       FLOAT8 mld = 1.25*(1-cos(PI*sb/SBPSY_s))-2.5;
  192.       mld_s[sb] = pow(10.0,mld);
  193.     }
  194.     for ( sb = 0; sb < SBPSY_l; sb++ ) {
  195.       FLOAT8 mld = 1.25*(1-cos(PI*sb/SBPSY_l))-2.5;
  196.       mld_l[sb] = pow(10.0,mld);
  197.     }
  198.     
  199.     for (i=0;i<HBLKSIZE;i++) partition_l[i]=-1;
  201.     L3para_read( (FLOAT8) gfp->out_samplerate,numlines_l,numlines_s,partition_l,minval,qthr_l,s3_l,s3_s,
  202.  qthr_s,SNR_s,
  203.  bu_l,bo_l,w1_l,w2_l, bu_s,bo_s,w1_s,w2_s );
  204.     
  205.     
  206.     /* npart_l_orig   = number of partition bands before convolution */
  207.     /* npart_l  = number of partition bands after convolution */
  208.     npart_l_orig=0; npart_s_orig=0;
  209.     for (i=0;i<HBLKSIZE;i++) 
  210.       if (partition_l[i]>npart_l_orig) npart_l_orig=partition_l[i];
  211.     npart_l_orig++;
  213.     for (i=0;numlines_s[i]>=0;i++)
  214.       ;
  215.     npart_s_orig = i;
  216.     
  217.     npart_l=bo_l[SBPSY_l-1]+1;
  218.     npart_s=bo_s[SBPSY_s-1]+1;
  220.     /* MPEG2 tables are screwed up 
  221.      * the mapping from paritition bands to scalefactor bands will use
  222.      * more paritition bands than we have.  
  223.      * So we will not compute these fictitious partition bands by reducing
  224.      * npart_l below.  */
  225.     if (npart_l > npart_l_orig) {
  226.       npart_l=npart_l_orig;
  227.       bo_l[SBPSY_l-1]=npart_l-1;
  228.       w2_l[SBPSY_l-1]=1.0;
  229.     }
  230.     if (npart_s > npart_s_orig) {
  231.       npart_s=npart_s_orig;
  232.       bo_s[SBPSY_s-1]=npart_s-1;
  233.       w2_s[SBPSY_s-1]=1.0;
  234.     }
  235.     
  236.     
  237.     
  238.     for (i=0; i<npart_l; i++) {
  239.       for (j = 0; j < npart_l_orig; j++) {
  240. if (s3_l[i][j] != 0.0)
  241.   break;
  242.       }
  243.       s3ind[i][0] = j;
  244.       
  245.       for (j = npart_l_orig - 1; j > 0; j--) {
  246. if (s3_l[i][j] != 0.0)
  247.   break;
  248.       }
  249.       s3ind[i][1] = j;
  250.     }
  253.     for (i=0; i<npart_s; i++) {
  254.       for (j = 0; j < npart_s_orig; j++) {
  255. if (s3_s[i][j] != 0.0)
  256.   break;
  257.       }
  258.       s3ind_s[i][0] = j;
  259.       
  260.       for (j = npart_s_orig - 1; j > 0; j--) {
  261. if (s3_s[i][j] != 0.0)
  262.   break;
  263.       }
  264.       s3ind_s[i][1] = j;
  265.     }
  266.     
  267.     
  268.     /*  
  269.       #include "debugscalefac.c"
  270.     */
  271.     
  273. #define AACS3
  274. #define NEWS3XX
  276. #ifdef AACS3
  277.     /* AAC values, results in more masking over MP3 values */
  278. # define TMN 18
  279. # define NMT 6
  280. #else
  281.     /* MP3 values */
  282. # define TMN 29
  283. # define NMT 6
  284. #endif
  286. #define rpelev 2
  287. #define rpelev2 16
  289.     /* compute norm_l, norm_s instead of relying on table data */
  290.     for ( b = 0;b < npart_l; b++ ) {
  291.       FLOAT8 norm=0;
  292.       for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ ) {
  293. norm += s3_l[b][k];
  294.       }
  295.       for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ ) {
  296. s3_l[b][k] *= exp(-LN_TO_LOG10 * NMT) / norm;
  297.       }
  298.       /*printf("%i  norm=%f  norm_l=%f n",b,1/norm,norm_l[b]);*/
  299.     }
  301.     /* MPEG1 SNR_s data is given in db, convert to energy */
  302.     if (gfp->version == 1) {
  303.       for ( b = 0;b < npart_s; b++ ) {
  304. SNR_s[b]=exp( (FLOAT8) SNR_s[b] * LN_TO_LOG10 );
  305.       }
  306.     }
  308.     for ( b = 0;b < npart_s; b++ ) {
  309.       FLOAT8 norm=0;
  310.       for ( k = s3ind_s[b][0]; k <= s3ind_s[b][1]; k++ ) {
  311. norm += s3_s[b][k];
  312.       }
  313.       for ( k = s3ind_s[b][0]; k <= s3ind_s[b][1]; k++ ) {
  314. s3_s[b][k] *= SNR_s[b] / norm;
  315.       }
  316.       /*printf("%i  norm=%f  norm_s=%f n",b,1/norm,norm_l[b]);*/
  317.     }
  318.     
  319.     init_fft();
  320.   }
  321.   /************************* End of Initialization *****************************/
  322.   
  325.   
  326.   
  327.   numchn = gfp->stereo;
  328.   /* chn=2 and 3 = Mid and Side channels */
  329.   if (gfp->mode == MPG_MD_JOINT_STEREO) numchn=4;
  330.   for (chn=0; chn<numchn; chn++) {
  331.   
  332.     wsamp_s = wsamp_S+(chn & 1);
  333.     wsamp_l = wsamp_L+(chn & 1);
  336.     if (chn<2) {    
  337.       /**********************************************************************
  338.        *  compute FFTs
  339.        **********************************************************************/
  340.       fft_long ( *wsamp_l, chn, buffer);
  341.       fft_short( *wsamp_s, chn, buffer); 
  342.       
  343.       /* LR maskings  */
  344.       percep_entropy[chn] = pe[chn]; 
  345.       masking_ratio[gr_out][chn].thm = thm[chn];
  346.       masking_ratio[gr_out][chn].en = en[chn];
  347.     }else{
  348.       /* MS maskings  */
  349.       percep_MS_entropy[chn-2] = pe[chn]; 
  350.       masking_MS_ratio[gr_out][chn-2].en = en[chn];
  351.       masking_MS_ratio[gr_out][chn-2].thm = thm[chn];
  352.       
  353.       if (chn == 2)
  354.       {
  355.         for (j = BLKSIZE-1; j >=0 ; --j)
  356.         {
  357.           FLOAT l = wsamp_L[0][j];
  358.           FLOAT r = wsamp_L[1][j];
  359.           wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
  360.           wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5);
  361.         }
  362.         for (b = 2; b >= 0; --b)
  363.         {
  364.           for (j = BLKSIZE_s-1; j >= 0 ; --j)
  365.           {
  366.             FLOAT l = wsamp_S[0][b][j];
  367.             FLOAT r = wsamp_S[1][b][j];
  368.             wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
  369.             wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
  370.           }
  371.         }
  372.       }
  373.     }
  375.     /**********************************************************************
  376.      *  compute energies
  377.      **********************************************************************/
  378.     
  379.     
  380.     
  381.     energy[0]  = (*wsamp_l)[0];
  382.     energy[0] *= energy[0];
  383.     
  384.     tot_ener[chn] = energy[0]; /* sum total energy at nearly no extra cost */
  385.     
  386.     for (j=BLKSIZE/2-1; j >= 0; --j)
  387.     {
  388.       FLOAT re = (*wsamp_l)[BLKSIZE/2-j];
  389.       FLOAT im = (*wsamp_l)[BLKSIZE/2+j];
  390.       energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
  391.       
  392.       tot_ener[chn] += energy[BLKSIZE/2-j];
  393.     }
  394.     for (b = 2; b >= 0; --b)
  395.     {
  396.       energy_s[b][0]  = (*wsamp_s)[b][0];
  397.       energy_s[b][0] *=  energy_s [b][0];
  398.       for (j=BLKSIZE_s/2-1; j >= 0; --j)
  399.       {
  400.         FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
  401.         FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
  402.         energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
  403.       }
  404.     }
  407. #ifdef HAVEGTK
  408.   if(gfp->gtkflag) {
  409.     for (j=0; j<HBLKSIZE ; j++) {
  410.       pinfo->energy[gr_out][chn][j]=energy_save[chn][j];
  411.       energy_save[chn][j]=energy[j];
  412.     }
  413.   }
  414. #endif
  415.     
  416.     /**********************************************************************
  417.      *    compute unpredicatability of first six spectral lines            * 
  418.      **********************************************************************/
  419.     for ( j = 0; j < cw_lower_index; j++ )
  420.       {  /* calculate unpredictability measure cw */
  421. FLOAT an, a1, a2;
  422. FLOAT bn, b1, b2;
  423. FLOAT rn, r1, r2;
  424. FLOAT numre, numim, den;
  426. a2 = ax_sav[chn][1][j];
  427. b2 = bx_sav[chn][1][j];
  428. r2 = rx_sav[chn][1][j];
  429. a1 = ax_sav[chn][1][j] = ax_sav[chn][0][j];
  430. b1 = bx_sav[chn][1][j] = bx_sav[chn][0][j];
  431. r1 = rx_sav[chn][1][j] = rx_sav[chn][0][j];
  432. an = ax_sav[chn][0][j] = (*wsamp_l)[j];
  433. bn = bx_sav[chn][0][j] = j==0 ? (*wsamp_l)[0] : (*wsamp_l)[BLKSIZE-j];  
  434. rn = rx_sav[chn][0][j] = sqrt(energy[j]);
  436. { /* square (x1,y1) */
  437.   if( r1 != 0 ) {
  438.     numre = (a1*b1);
  439.     numim = (a1*a1-b1*b1)*(FLOAT)0.5;
  440.     den = r1*r1;
  441.   } else {
  442.     numre = 1;
  443.     numim = 0;
  444.     den = 1;
  445.   }
  446. }
  448. { /* multiply by (x2,-y2) */
  449.   if( r2 != 0 ) {
  450.     FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
  451.     FLOAT tmp1 = -a2*numre+tmp2;
  452.     numre =       -b2*numim+tmp2;
  453.     numim = tmp1;
  454.     den *= r2;
  455.   } else {
  456.     /* do nothing */
  457.   }
  458. }
  460. { /* r-prime factor */
  461.   FLOAT tmp = (2*r1-r2)/den;
  462.   numre *= tmp;
  463.   numim *= tmp;
  464. }
  465. den=rn+fabs(2*r1-r2);
  466. if( den != 0 ) {
  467.   numre = (an+bn)*(FLOAT)0.5-numre;
  468.   numim = (an-bn)*(FLOAT)0.5-numim;
  469.   den = sqrt(numre*numre+numim*numim)/den;
  470. }
  471. cw[j] = den;
  472.       }
  476.     /**********************************************************************
  477.      *     compute unpredicatibility of next 200 spectral lines            *
  478.      **********************************************************************/ 
  479.     for ( j = cw_lower_index; j < cw_upper_index; j += 4 )
  480.       {/* calculate unpredictability measure cw */
  481. FLOAT rn, r1, r2;
  482. FLOAT numre, numim, den;
  484. k = (j+2) / 4; 
  486. { /* square (x1,y1) */
  487.   r1 = energy_s[0][k];
  488.   if( r1 != 0 ) {
  489.     FLOAT a1 = (*wsamp_s)[0][k]; 
  490.     FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */
  491.     numre = (a1*b1);
  492.     numim = (a1*a1-b1*b1)*(FLOAT)0.5;
  493.     den = r1;
  494.     r1 = sqrt(r1);
  495.   } else {
  496.     numre = 1;
  497.     numim = 0;
  498.     den = 1;
  499.   }
  500. }
  503. { /* multiply by (x2,-y2) */
  504.   r2 = energy_s[2][k];
  505.   if( r2 != 0 ) {
  506.     FLOAT a2 = (*wsamp_s)[2][k]; 
  507.     FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k];
  508.     
  509.     
  510.     FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
  511.     FLOAT tmp1 = -a2*numre+tmp2;
  512.     numre =       -b2*numim+tmp2;
  513.     numim = tmp1;
  514.     
  515.     r2 = sqrt(r2);
  516.     den *= r2;
  517.   } else {
  518.     /* do nothing */
  519.   }
  520. }
  522. { /* r-prime factor */
  523.   FLOAT tmp = (2*r1-r2)/den;
  524.   numre *= tmp;
  525.   numim *= tmp;
  526. }
  528. rn = sqrt(energy_s[1][k]);
  529. den=rn+fabs(2*r1-r2);
  530. if( den != 0 ) {
  531.   FLOAT an = (*wsamp_s)[1][k]; 
  532.   FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k];
  533.   numre = (an+bn)*(FLOAT)0.5-numre;
  534.   numim = (an-bn)*(FLOAT)0.5-numim;
  535.   den = sqrt(numre*numre+numim*numim)/den;
  536. }
  538. cw[j+1] = cw[j+2] = cw[j+3] = cw[j] = den;
  539.       }
  540.     
  541. #if 0
  542.     for ( j = 14; j < HBLKSIZE-4; j += 4 )
  543.       {/* calculate energy from short ffts */
  544. FLOAT8 tot,ave;
  545. k = (j+2) / 4; 
  546. for (tot=0, sblock=0; sblock < 3; sblock++)
  547.   tot+=energy_s[sblock][k];
  548. ave = energy[j+1]+ energy[j+2]+ energy[j+3]+ energy[j];
  549. ave /= 4.;
  550. /*
  551.   printf("energy / tot %i %5.2f   %e  %en",j,ave/(tot*16./3.),
  552.   ave,tot*16./3.);
  553. */
  554. energy[j+1] = energy[j+2] = energy[j+3] =  energy[j]=tot;
  555.       }
  556. #endif
  557.     
  558.     
  559.     
  560.     
  561.     
  562.     
  563.     
  564.     
  565.     /**********************************************************************
  566.      *    Calculate the energy and the unpredictability in the threshold   *
  567.      *    calculation partitions                                           *
  568.      **********************************************************************/
  569. #if 0
  570.     for ( b = 0; b < CBANDS; b++ )
  571.       {
  572. eb[b] = 0;
  573. cb[b] = 0;
  574.       }
  575.     for ( j = 0; j < HBLKSIZE; j++ )
  576.       {
  577. int tp = partition_l[j];
  579. if ( tp >= 0 )
  580.   {
  581.     eb[tp] += energy[j];
  582.     cb[tp] += cw[j] * energy[j];
  583.   }
  584. assert(tp<npart_l_orig);
  585.       }
  586. #else
  587.     b = 0;
  588.     for (j = 0; j < cw_upper_index;)
  589.       {
  590. FLOAT8 ebb, cbb;
  591. int i;
  593. ebb = energy[j];
  594. cbb = energy[j] * cw[j];
  595. j++;
  597. for (i = numlines_l[b] - 1; i > 0; i--)
  598.   {
  599.     ebb += energy[j];
  600.     cbb += energy[j] * cw[j];
  601.     j++;
  602.   }
  603. eb[b] = ebb;
  604. cb[b] = cbb;
  605. b++;
  606.       }
  608.     for (; b < npart_l_orig; b++ )
  609.       {
  610. int i;
  611. FLOAT8 ebb = energy[j++];
  613. for (i = numlines_l[b] - 1; i > 0; i--)
  614.   {
  615.     ebb += energy[j++];
  616.   }
  617. eb[b] = ebb;
  618. cb[b] = ebb * 0.4;
  619.       }
  620. #endif
  622.     /**********************************************************************
  623.      *      convolve the partitioned energy and unpredictability           *
  624.      *      with the spreading function, s3_l[b][k]                        *
  625.      ******************************************************************** */
  626.     pe[chn] = 0; /*  calculate percetual entropy */
  627.     for ( b = 0;b < npart_l; b++ )
  628.       {
  629. FLOAT8 tbb,ecb,ctb;
  630. FLOAT8 temp_1; /* BUG of IS */
  632. ecb = 0;
  633. ctb = 0;
  634. for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ )
  635.   {
  636.     ecb += s3_l[b][k] * eb[k]; /* sprdngf for Layer III */
  637.     ctb += s3_l[b][k] * cb[k];
  638.   }
  640. /* calculate the tonality of each threshold calculation partition */
  641. /* calculate the SNR in each threshhold calculation partition */
  643. tbb = ecb;
  644. if (tbb != 0)
  645.   {
  646.     tbb = ctb / tbb;
  647.     if (tbb <= 0.04875584301)
  648.       {
  649. tbb = exp(-LN_TO_LOG10 * (TMN - NMT));
  650.       }
  651.     else if (tbb > 0.4989003827)
  652.       {
  653. tbb = 1;
  654.       }
  655.     else
  656.       {
  657. tbb = log(tbb);
  658. tbb = exp(((TMN - NMT)*(LN_TO_LOG10*0.299))
  659. + ((TMN - NMT)*(LN_TO_LOG10*0.43 ))*tbb);  /* conv1=-0.299, conv2=-0.43 */
  660.       }
  661.   }
  663. tbb = Min(minval[b], tbb);
  664. ecb *= tbb;
  666. /* pre-echo control */
  667. /* rpelev=2.0, rpelev2=16.0 */
  668. temp_1 = Min(ecb, Min(rpelev*nb_1[chn][b],rpelev2*nb_2[chn][b]) );
  669. thr[b] = Max( qthr_l[b], temp_1 ); 
  670. nb_2[chn][b] = nb_1[chn][b];
  671. nb_1[chn][b] = ecb;
  673. /* note: all surges in PE are because of the above pre-echo formula
  674.  * for temp_1.  it this is not used, PE is always around 600
  675.  */
  677. if (thr[b] < eb[b])
  678.   {
  679.     /* there's no non sound portition, because thr[b] is
  680.      maximum of qthr_l and temp_1 */
  681.     pe[chn] -= numlines_l[b] * log(thr[b] / eb[b]);
  682.   }
  683.       }
  686. #ifdef HAVEGTK
  687.     if (gfp->gtkflag) {
  688.       FLOAT mn,mx,ma=0,mb=0,mc=0;
  690.       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
  691.       {
  692.         ma += energy_s[0][j];
  693.         mb += energy_s[1][j];
  694.         mc += energy_s[2][j];
  695.       }
  696.       mn = Min(ma,mb);
  697.       mn = Min(mn,mc);
  698.       mx = Max(ma,mb);
  699.       mx = Max(mx,mc);
  701.       pinfo->ers[gr_out][chn]=ers_save[chn];
  702.       ers_save[chn]=mx/(1e-12+mn);
  703.       pinfo->pe[gr_out][chn]=pe_save[chn];
  704.       pe_save[chn]=pe[chn];
  705.     }
  706. #endif
  707.     
  708.     /*************************************************************** 
  709.      * determine the block type (window type) based on L & R channels
  710.      * 
  711.      ***************************************************************/
  712.     if (chn<2) {
  713.       if (gfp->no_short_blocks){
  714. uselongblock[chn]=1;
  715.       } else {
  716. /* tuned for t1.wav.  doesnt effect most other samples */
  717. if (pe[chn] > 3000) {
  718.   uselongblock[chn]=0;
  719. } else { 
  720.   FLOAT mn,mx,ma=0,mb=0,mc=0;
  722.   for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
  723.   {
  724.       ma += energy_s[0][j];
  725.       mb += energy_s[1][j];
  726.       mc += energy_s[2][j];
  727.   }
  728.   mn = Min(ma,mb);
  729.   mn = Min(mn,mc);
  730.   mx = Max(ma,mb);
  731.   mx = Max(mx,mc);
  733.   uselongblock[chn] = 1;
  734.   
  735.   if ( mx > 30*mn ) 
  736.   {/* big surge of energy - always use short blocks */
  737.     uselongblock[chn] = 0;
  738.   } 
  739.   else if ((mx > 10*mn) && (pe[chn] > 1000))
  740.   {/* medium surge, medium pe - use short blocks */
  741.     uselongblock[chn] = 0;
  742.   }
  743. } 
  744.       }
  745.     }
  749.     /*************************************************************** 
  750.      * compute masking thresholds for both short and long blocks
  751.      ***************************************************************/
  752.     /* longblock threshold calculation (part 2) */
  753.     for ( sb = 0; sb < SBPSY_l; sb++ )
  754.       {
  755. FLOAT8 enn = w1_l[sb] * eb[bu_l[sb]] + w2_l[sb] * eb[bo_l[sb]];
  756. FLOAT8 thmm = w1_l[sb] *thr[bu_l[sb]] + w2_l[sb] * thr[bo_l[sb]];
  757. for ( b = bu_l[sb]+1; b < bo_l[sb]; b++ )
  758.   {
  759.     enn  += eb[b];
  760.     thmm += thr[b];
  761.   }
  762. en[chn].l[sb] = enn;
  763. thm[chn].l[sb] = thmm;
  764.       }
  765.     
  766.     
  767.     /* threshold calculation for short blocks */
  768.     for ( sblock = 0; sblock < 3; sblock++ )
  769.       {
  770. j = 0;
  771. for ( b = 0; b < npart_s_orig; b++ )
  772.   {
  773.     int i;
  774.     FLOAT ecb = energy_s[sblock][j++];
  775.     for (i = numlines_s[b]; i > 0; i--)
  776.       {
  777. ecb += energy_s[sblock][j++];
  778.       }
  779.     eb[b] = ecb;
  780.   }
  782. for ( b = 0; b < npart_s; b++ )
  783.   {
  784.     FLOAT8 ecb = 0;
  785.     for ( k = s3ind_s[b][0]; k <= s3ind_s[b][1]; k++ )
  786.       {
  787. ecb += s3_s[b][k] * eb[k];
  788.       }
  789.     thr[b] = Max (qthr_s[b], ecb);
  790.   }
  792. for ( sb = 0; sb < SBPSY_s; sb++ )
  793.   {
  794.     FLOAT8 enn  = w1_s[sb] * eb[bu_s[sb]] + w2_s[sb] * eb[bo_s[sb]];
  795.     FLOAT8 thmm = w1_s[sb] *thr[bu_s[sb]] + w2_s[sb] * thr[bo_s[sb]];
  796.     for ( b = bu_s[sb]+1; b < bo_s[sb]; b++ )
  797.       {
  798. enn  += eb[b];
  799. thmm += thr[b];
  800.       }
  801.     en[chn].s[sb][sblock] = enn;
  802.     thm[chn].s[sb][sblock] = thmm;
  803.   }
  804.       }
  805.   } /* end loop over chn */
  808.   /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */
  809.   if ( numchn==4 /* mid/side and r/l */) {
  810.     FLOAT8 rside,rmid,mld;
  811.     int chmid=2,chside=3; 
  812.     
  813.     for ( sb = 0; sb < SBPSY_l; sb++ ) {
  814.       /* use this fix if L & R masking differs by 2db or less */
  815.       /* if db = 10*log10(x2/x1) < 2 */
  816.       /* if (x2 < 1.58*x1) { */
  817.       if (thm[0].l[sb] <= 1.58*thm[1].l[sb]
  818.   && thm[1].l[sb] <= 1.58*thm[0].l[sb]) {
  820. mld = mld_l[sb]*en[chside].l[sb];
  821. rmid = Max(thm[chmid].l[sb], Min(thm[chside].l[sb],mld));
  823. mld = mld_l[sb]*en[chmid].l[sb];
  824. rside = Max(thm[chside].l[sb],Min(thm[chmid].l[sb],mld));
  826. thm[chmid].l[sb]=rmid;
  827. thm[chside].l[sb]=rside;
  828.       }
  829.     }
  830.     for ( sb = 0; sb < SBPSY_s; sb++ ) {
  831.       for ( sblock = 0; sblock < 3; sblock++ ) {
  832. if (thm[0].s[sb][sblock] <= 1.58*thm[1].s[sb][sblock]
  833.     && thm[1].s[sb][sblock] <= 1.58*thm[0].s[sb][sblock]) {
  835.   mld = mld_s[sb]*en[chside].s[sb][sblock];
  836.   rmid = Max(thm[chmid].s[sb][sblock],Min(thm[chside].s[sb][sblock],mld));
  838.   mld = mld_s[sb]*en[chmid].s[sb][sblock];
  839.   rside = Max(thm[chside].s[sb][sblock],Min(thm[chmid].s[sb][sblock],mld));
  841.   thm[chmid].s[sb][sblock]=rmid;
  842.   thm[chside].s[sb][sblock]=rside;
  843. }
  844.       }
  845.     }
  846.   }
  849.   
  851.   
  852.   
  853.   if (gfp->mode == MPG_MD_JOINT_STEREO)  {
  854.     /* determin ms_ratio from masking thresholds*/
  855.     /* use ms_stereo (ms_ratio < .35) if average thresh. diff < 5 db */
  856.     FLOAT8 db,x1,x2,sidetot=0,tot=0;
  857.     for (sb= SBPSY_l/4 ; sb< SBPSY_l; sb ++ ) {
  858.       x1 = Min(thm[0].l[sb],thm[1].l[sb]);
  859.       x2 = Max(thm[0].l[sb],thm[1].l[sb]);
  860.       /* thresholds difference in db */
  861.       if (x2 >= 1000*x1)  db=3;
  862.       else db = log10(x2/x1);  
  863.       /*  printf("db = %f %e %e  n",db,thm[0].l[sb],thm[1].l[sb]);*/
  864.       sidetot += db;
  865.       tot++;
  866.     }
  867.     ms_ratio_l= (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
  868.     ms_ratio_l = Min(ms_ratio_l,0.5);
  869.     
  870.     sidetot=0; tot=0;
  871.     for ( sblock = 0; sblock < 3; sblock++ )
  872.       for ( sb = SBPSY_s/4; sb < SBPSY_s; sb++ ) {
  873. x1 = Min(thm[0].s[sb][sblock],thm[1].s[sb][sblock]);
  874. x2 = Max(thm[0].s[sb][sblock],thm[1].s[sb][sblock]);
  875. /* thresholds difference in db */
  876. if (x2 >= 1000*x1)  db=3;
  877. else db = log10(x2/x1);  
  878. sidetot += db;
  879. tot++;
  880.       }
  881.     ms_ratio_s = (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
  882.     ms_ratio_s = Min(ms_ratio_s,.5);
  883.   }
  885.   /*************************************************************** 
  886.    * determin final block type
  887.    ***************************************************************/
  889.   for (chn=0; chn<gfp->stereo; chn++) {
  890.     blocktype[chn] = NORM_TYPE;
  891.   }
  894.   if (gfp->stereo==2) {
  895.     if (!gfp->allow_diff_short || gfp->mode==MPG_MD_JOINT_STEREO) {
  896.       /* force both channels to use the same block type */
  897.       /* this is necessary if the frame is to be encoded in ms_stereo.  */
  898.       /* But even without ms_stereo, FhG  does this */
  899.       int bothlong= (uselongblock[0] && uselongblock[1]);
  900.       if (!bothlong) {
  901. uselongblock[0]=0;
  902. uselongblock[1]=0;
  903.       }
  904.     }
  905.   }
  907.   
  908.   
  909.   /* update the blocktype of the previous granule, since it depends on what
  910.    * happend in this granule */
  911.   for (chn=0; chn<gfp->stereo; chn++) {
  912.     if ( uselongblock[chn])
  913.       { /* no attack : use long blocks */
  914. switch( blocktype_old[chn] ) 
  915.   {
  916.   case NORM_TYPE:
  917.   case STOP_TYPE:
  918.     blocktype[chn] = NORM_TYPE;
  919.     break;
  920.   case SHORT_TYPE:
  921.     blocktype[chn] = STOP_TYPE; 
  922.     break;
  923.   case START_TYPE:
  924.     fprintf( stderr, "Error in block selectingn" );
  925.     abort();
  926.     break; /* problem */
  927.   }
  928.       } else   {
  929. /* attack : use short blocks */
  930. blocktype[chn] = SHORT_TYPE;
  931. if ( blocktype_old[chn] == NORM_TYPE ) {
  932.   blocktype_old[chn] = START_TYPE;
  933. }
  934. if ( blocktype_old[chn] == STOP_TYPE ) {
  935.   blocktype_old[chn] = SHORT_TYPE ;
  936. }
  937.       }
  938.     
  939.     blocktype_d[chn] = blocktype_old[chn];  /* value returned to calling program */
  940.     blocktype_old[chn] = blocktype[chn];    /* save for next call to l3psy_anal */
  941.   }
  942.   
  943.   if (blocktype_d[0]==2) 
  944.     *ms_ratio = ms_ratio_s_old;
  945.   else
  946.     *ms_ratio = ms_ratio_l_old;
  948.   ms_ratio_s_old = ms_ratio_s;
  949.   ms_ratio_l_old = ms_ratio_l;
  951.   /* we dont know the block type of this frame yet - assume long */
  952.   *ms_ratio_next = ms_ratio_l;
  956.   /*********************************************************************/
  957.   /* compute side_energy / (side+mid)_energy */
  958.   /* 0 = no energy in side channel */
  959.   /* .5 = half of total energy in side channel */
  960.   /*********************************************************************/
  961.   if (numchn==4)  {
  962.     FLOAT tmp = tot_ener[3]+tot_ener[2];
  963.     *ms_ener_ratio = ms_ener_ratio_old;
  964.     ms_ener_ratio_old=0;
  965.     if (tmp>0) ms_ener_ratio_old=tot_ener[3]/tmp;
  966.   } else
  967.     /* we didn't compute ms_ener_ratios */
  968.     *ms_ener_ratio = 0;
  969.  
  970. }
  977. void L3para_read(FLOAT8 sfreq, int *numlines_l,int *numlines_s, int *partition_l, FLOAT8 *minval,
  978. FLOAT8 *qthr_l, FLOAT8 s3_l[64][64], FLOAT8 s3_s[CBANDS + 1][CBANDS + 1],
  979. FLOAT8 *qthr_s, FLOAT8 *SNR, 
  980. int *bu_l, int *bo_l, FLOAT8 *w1_l, FLOAT8 *w2_l, 
  981. int *bu_s, int *bo_s, FLOAT8 *w1_s, FLOAT8 *w2_s)
  982. {
  983.   FLOAT8 freq_tp;
  984.   FLOAT8 bval_l[CBANDS], bval_s[CBANDS];
  985.   int   cbmax=0, cbmax_tp;
  986.   FLOAT8 *p = psy_data;
  988.   int  sbmax ;
  989.   int  i,j,k,k2,loop, part_max ;
  990.   int freq_scale=1;
  993.   /* use MPEG1 tables.  The MPEG2 tables in tables.c appear to be 
  994.    * junk.  MPEG2 doc claims data for these tables is the same as the
  995.    * MPEG1 data for 2x sampling frequency */
  996.   /*  if (sfreq<32000) freq_scale=2; */
  997.   
  1000.   /* Read long block data */
  1002.   for(loop=0;loop<6;loop++)
  1003.     {
  1004.       freq_tp = *p++;
  1005.       cbmax_tp = (int) *p++;
  1006.       cbmax_tp++;
  1008.       if (sfreq == freq_tp/freq_scale )
  1009. {
  1010.   cbmax = cbmax_tp;
  1011.   for(i=0,k2=0;i<cbmax_tp;i++)
  1012.     {
  1013.       j = (int) *p++;
  1014.       numlines_l[i] = (int) *p++;
  1015.       minval[i] = exp(-((*p++) - NMT) * LN_TO_LOG10);
  1016.       qthr_l[i] = *p++;
  1017.       /* norm_l[i] = *p++*/ p++;
  1018.       bval_l[i] = *p++;
  1019.       if (j!=i)
  1020. {
  1021.   fprintf(stderr,"1. please check "psy_data"");
  1022.   exit(-1);
  1023. }
  1024.       for(k=0;k<numlines_l[i];k++)
  1025. partition_l[k2++] = i ;
  1026.     }
  1027. }
  1028.       else
  1029. p += cbmax_tp * 6;
  1030.     }
  1032. #define NEWBARKXXX
  1033. #ifdef NEWBARK
  1034.   /* compute bark values of each critical band */
  1035.   j = 0;
  1036.   for(i=0;i<cbmax;i++)
  1037.     {
  1038.       FLOAT8 ji, freq, bark;
  1040.       ji = j + (numlines_l[i]-1)/2.0;
  1041.       freq = sfreq*ji/1024000.0;
  1042.       bark = 13*atan(.76*freq) + 3.5*atan(freq*freq/(7.5*7.5));
  1044.       printf("%i %i bval_l table=%f  f=%f  formaula=%f n",i,j,bval_l[i],freq,bark);
  1045.       bval_l[i]=bark;
  1046.       j += numlines_l[i];
  1047.     }
  1048. #endif
  1050.   /************************************************************************
  1051.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1052.    * ing function, centered at band j, for band i, store for later use    *
  1053.    ************************************************************************/
  1054.   /* i.e.: sum over j to spread into signal barkval=i  
  1055.      NOTE: i and j are used opposite as in the ISO docs */
  1056.   part_max = cbmax ;
  1057.   for(i=0;i<part_max;i++)
  1058.     {
  1059.       FLOAT8 tempx,x,tempy,temp;
  1060.       for(j=0;j<part_max;j++)
  1061. {
  1062.   /*tempx = (bval_l[i] - bval_l[j])*1.05;*/
  1063.   if (j>=i) tempx = (bval_l[i] - bval_l[j])*3.0;
  1064.   else    tempx = (bval_l[i] - bval_l[j])*1.5;
  1066. #ifdef AACS3
  1067.           if (i>=j) tempx = (bval_l[i] - bval_l[j])*3.0;
  1068.   else    tempx = (bval_l[i] - bval_l[j])*1.5; 
  1069. #endif
  1071.   if(tempx>=0.5 && tempx<=2.5)
  1072.     {
  1073.       temp = tempx - 0.5;
  1074.       x = 8.0 * (temp*temp - 2.0 * temp);
  1075.     }
  1076.   else x = 0.0;
  1077.   tempx += 0.474;
  1078.   tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1080. #ifdef NEWS3
  1081.   if (j>=i) tempy = (bval_l[j] - bval_l[i])*(-15);
  1082.   else    tempy = (bval_l[j] - bval_l[i])*25;
  1083.   x=0; 
  1084. #endif
  1085.   /*
  1086.   if ((i==part_max/2)  && (fabs(bval_l[j] - bval_l[i])) < 3) {
  1087.     printf("bark=%f   x+tempy = %f  n",bval_l[j] - bval_l[i],x+tempy);
  1088.   }
  1089.   */
  1091.   if (tempy <= -60.0) s3_l[i][j] = 0.0;
  1092.   else                s3_l[i][j] = exp( (x + tempy)*LN_TO_LOG10 ); 
  1093. }
  1094.     }
  1096.   /* Read short block data */
  1097.   for(loop=0;loop<6;loop++)
  1098.     {
  1099.       freq_tp = *p++;
  1100.       cbmax_tp = (int) *p++;
  1101.       cbmax_tp++;
  1103.       if (sfreq == freq_tp/freq_scale )
  1104. {
  1105.   cbmax = cbmax_tp;
  1106.   for(i=0,k2=0;i<cbmax_tp;i++)
  1107.     {
  1108.       j = (int) *p++;
  1109.       numlines_s[i] = (int) *p++;
  1110.       qthr_s[i] = *p++;         
  1111.       /* norm_s[i] =*p++ */ p++;         
  1112.       SNR[i] = *p++;            
  1113.       bval_s[i] = *p++;
  1114.       if (j!=i)
  1115. {
  1116.   fprintf(stderr,"3. please check "psy_data"");
  1117.   exit(-1);
  1118. }
  1119.       numlines_s[i]--;
  1120.     }
  1121.   numlines_s[i] = -1;
  1122. }
  1123.       else
  1124. p += cbmax_tp * 6;
  1125.     }
  1128. #ifdef NEWBARK
  1129.   /* compute bark values of each critical band */
  1130.   j = 0;
  1131.   for(i=0;i<cbmax;i++)
  1132.     {
  1133.       FLOAT8 ji, freq, bark;
  1134.       ji = (j * 2 + numlines_s[i]) / 2.0;
  1135.       freq = sfreq*ji/256000.0;
  1136.       bark = 13*atan(.76*freq) + 3.5*atan(freq*freq/(7.5*7.5));
  1137.       printf("%i %i bval_s = %f  %f  %f n",i,j,bval_s[i],freq,bark);
  1138.       bval_s[i]=bark;
  1139.       j += numlines_s[i] + 1;
  1140.     }
  1141. #endif
  1145.   /************************************************************************
  1146.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1147.    * ing function, centered at band j, for band i, store for later use    *
  1148.    ************************************************************************/
  1149.   part_max = cbmax ;
  1150.   for(i=0;i<part_max;i++)
  1151.     {
  1152.       FLOAT8 tempx,x,tempy,temp;
  1153.       for(j=0;j<part_max;j++)
  1154. {
  1155.   /* tempx = (bval_s[i] - bval_s[j])*1.05;*/
  1156.   if (j>=i) tempx = (bval_s[i] - bval_s[j])*3.0;
  1157.   else    tempx = (bval_s[i] - bval_s[j])*1.5;
  1158. #ifdef AACS3
  1159.           if (i>=j) tempx = (bval_s[i] - bval_s[j])*3.0;
  1160.   else    tempx = (bval_s[i] - bval_s[j])*1.5; 
  1161. #endif
  1162.   if(tempx>=0.5 && tempx<=2.5)
  1163.     {
  1164.       temp = tempx - 0.5;
  1165.       x = 8.0 * (temp*temp - 2.0 * temp);
  1166.     }
  1167.   else x = 0.0;
  1168.   tempx += 0.474;
  1169.   tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1170. #ifdef NEWS3
  1171.   if (j>=i) tempy = (bval_s[j] - bval_s[i])*(-15);
  1172.   else    tempy = (bval_s[j] - bval_s[i])*25;
  1173.   x=0; 
  1174. #endif
  1175.   if (tempy <= -60.0) s3_s[i][j] = 0.0;
  1176.   else                s3_s[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
  1177. }
  1178.     }
  1179.   /* Read long block data for converting threshold calculation 
  1180.      partitions to scale factor bands */
  1182.   for(loop=0;loop<6;loop++)
  1183.     {
  1184.       freq_tp = *p++;
  1185.       sbmax =  (int) *p++;
  1186.       sbmax++;
  1188.       if (sfreq == freq_tp/freq_scale)
  1189. {
  1190.   for(i=0;i<sbmax;i++)
  1191.     {
  1192.       j = (int) *p++;
  1193.       p++;             
  1194.       bu_l[i] = (int) *p++;
  1195.       bo_l[i] = (int) *p++;
  1196.       w1_l[i] = (FLOAT8) *p++;
  1197.       w2_l[i] = (FLOAT8) *p++;
  1198.       if (j!=i)
  1199. { fprintf(stderr,"30:please check "psy_data"n");
  1200. exit(-1);
  1201. }
  1203.       if (i!=0)
  1204. if ( (fabs(1.0-w1_l[i]-w2_l[i-1]) > 0.01 ) )
  1205.   {
  1206.     fprintf(stderr,"31l: please check "psy_data."n");
  1207.                   fprintf(stderr,"w1,w2: %f %f n",w1_l[i],w2_l[i-1]);
  1208.     exit(-1);
  1209.   }
  1210.     }
  1211. }
  1212.       else
  1213. p += sbmax * 6;
  1214.     }
  1216.   /* Read short block data for converting threshold calculation 
  1217.      partitions to scale factor bands */
  1219.   for(loop=0;loop<6;loop++)
  1220.     {
  1221.       freq_tp = *p++;
  1222.       sbmax = (int) *p++;
  1223.       sbmax++;
  1225.       if (sfreq == freq_tp/freq_scale)
  1226. {
  1227.   for(i=0;i<sbmax;i++)
  1228.     {
  1229.       j = (int) *p++;
  1230.       p++;
  1231.       bu_s[i] = (int) *p++;
  1232.       bo_s[i] = (int) *p++;
  1233.       w1_s[i] = *p++;
  1234.       w2_s[i] = *p++;
  1235.       if (j!=i)
  1236. { fprintf(stderr,"30:please check "psy_data"n");
  1237. exit(-1);
  1238. }
  1240.       if (i!=0)
  1241. if ( (fabs(1.0-w1_s[i]-w2_s[i-1]) > 0.01 ) )
  1242.   { 
  1243.                   fprintf(stderr,"31s: please check "psy_data."n");
  1244.                   fprintf(stderr,"w1,w2: %f %f n",w1_s[i],w2_s[i-1]);
  1245.   exit(-1);
  1246.   }
  1247.     }
  1248. }
  1249.       else
  1250. p += sbmax * 6;
  1251.     }
  1253. }