home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 13 / AACD13.ISO / AACD / Sound / LAME / src / psymodel.c < prev    next >
C/C++ Source or Header  |  2000-07-30  |  43KB  |  1,625 lines
  1. /*
  2.  *    psymodel.c
  3.  *
  4.  *    Copyright (c) 1999 Mark Taylor
  5.  *
  6.  * This library is free software; you can redistribute it and/or
  7.  * modify it under the terms of the GNU Library General Public
  8.  * License as published by the Free Software Foundation; either
  9.  * version 2 of the License, or (at your option) any later version.
  10.  *
  11.  * This library is distributed in the hope that it will be useful,
  12.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.     See the GNU
  14.  * Library General Public License for more details.
  15.  *
  16.  * You should have received a copy of the GNU Library General Public
  17.  * License along with this library; if not, write to the
  18.  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19.  * Boston, MA 02111-1307, USA.
  20.  */
  21.  
  22. #include "util.h"
  23. #include "encoder.h"
  24. #include "psymodel.h"
  25. #include "l3side.h"
  26. #include <assert.h>
  27. #include "tables.h"
  28. #include "fft.h"
  29.  
  30. #ifdef M_LN10
  31. #define        LN_TO_LOG10        (M_LN10/10)
  32. #else
  33. #define         LN_TO_LOG10             0.2302585093
  34. #endif
  35.  
  36.  
  37. int L3para_read( lame_global_flags *gfp,
  38.           FLOAT8 sfreq, int numlines[CBANDS],int numlines_s[CBANDS], 
  39.           FLOAT8 minval[CBANDS], 
  40.           FLOAT8 s3_l[CBANDS][CBANDS],
  41.           FLOAT8 s3_s[CBANDS][CBANDS],
  42.           FLOAT8 SNR_s[CBANDS],
  43.           int bu_l[SBPSY_l], int bo_l[SBPSY_l],
  44.           FLOAT8 w1_l[SBPSY_l], FLOAT8 w2_l[SBPSY_l],
  45.           int bu_s[SBPSY_s], int bo_s[SBPSY_s],
  46.           FLOAT8 w1_s[SBPSY_s], FLOAT8 w2_s[SBPSY_s],
  47.           int *, int *, int *, int *);
  48.  
  49. /* addition of simultaneous masking   Naoki Shibata 2000/7 */
  50. INLINE FLOAT8 mask_add(double m1,double m2,int b)
  51. {
  52.   static double table1[] = {
  53.     3.3246,3.23837,3.15437,3.00412,2.86103,2.65407,2.46209,2.284,
  54.     2.11879,1.96552,1.82335,1.69146,1.56911,1.46658,1.37074,1.31036,
  55.     1.25264,1.20648,1.16203,1.12765,1.09428,1.0659,1.03826,1.01895,
  56.     1
  57.   };
  58.  
  59.   static double table2[] = {
  60.     1.33352,1.35879,1.38454,1.39497,1.40548,1.3537,1.30382,1.22321,
  61.     1.14758
  62.   };
  63.  
  64.   static double table3[] = {
  65.     2.35364,2.29259,2.23313,2.12675,2.02545,1.87894,1.74303,1.61695,
  66.     1.49999,1.39148,1.29083,1.19746,1.11084,1.03826
  67.   };
  68.  
  69.  
  70.   int i;
  71.  
  72.   if (m1 == 0) return m2;
  73.  
  74.   if (b < 0) b = -b;
  75.  
  76.   i = 20*log10(m2 / m1)/20*16;
  77.   if (i < 0) i = -i;
  78.  
  79.   if (b <= 3) {
  80.     if (i > 8) return m1+m2;
  81.     if (m1+m2 < 200) {
  82.       double r;
  83.       r = (m1+m2)/200;
  84.       return (m1+m2)*(table2[i]*r+1*(1-r));
  85.     }
  86.     return (m1+m2)*table2[i];
  87.   }
  88.  
  89.   if (m1+m2 < 200) {
  90.     if (m1+m2 > 150) {
  91.       double f=1.0,r;
  92.       if (i > 24) return m1+m2;
  93.       if (i > 13) f = 1; else f = table3[i];
  94.       r = (m1+m2-150)/50;
  95.       return (m1+m2)*(table1[i]*r+f*(1-r));
  96.     }
  97.     if (i > 13) return m1+m2;
  98.     return (m1+m2)*table3[i];
  99.   }
  100.  
  101.   if (i > 24) return m1+m2;
  102.   return (m1+m2)*table1[i];
  103. }
  104.  
  105. int L3psycho_anal( lame_global_flags *gfp,
  106.                     short int *buffer[2],int gr_out , 
  107.                     FLOAT8 *ms_ratio,
  108.                     FLOAT8 *ms_ratio_next,
  109.             FLOAT8 *ms_ener_ratio,
  110.             III_psy_ratio masking_ratio[2][2],
  111.             III_psy_ratio masking_MS_ratio[2][2],
  112.             FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2], 
  113.                     int blocktype_d[2])
  114. {
  115.   lame_internal_flags *gfc=gfp->internal_flags;
  116.  
  117. /* to get a good cache performance, one has to think about
  118.  * the sequence, in which the variables are used.  
  119.  * (Note: these static variables have been moved to the gfc-> struct,
  120.  * and their order in memory is layed out in util.h)
  121.  */
  122.   
  123.  
  124.   /* fft and energy calculation   */
  125.   FLOAT (*wsamp_l)[BLKSIZE];
  126.   FLOAT (*wsamp_s)[3][BLKSIZE_s];
  127.   FLOAT tot_ener[4];
  128.  
  129.   /* convolution   */
  130.   FLOAT8 eb[CBANDS],eb2[CBANDS];
  131.   FLOAT8 cb[CBANDS];
  132.   FLOAT8 thr[CBANDS];
  133.   
  134.   FLOAT8 max[CBANDS],avg[CBANDS],tonality2[CBANDS];
  135.  
  136.   /* ratios    */
  137.   FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
  138.  
  139.   /* block type  */
  140.   int blocktype[2],uselongblock[2];
  141.   
  142.   /* usual variables like loop indices, etc..    */
  143.   int numchn, chn, samplerate;
  144.   int   b, i, j, k;
  145.   int    sb,sblock;
  146.   FLOAT cwlimit;
  147.  
  148.   /*  use a simplified spreading function: */
  149.     /*#define NEWS3  */ 
  150. #if 1
  151.     /* AAC values, results in more masking over MP3 values */
  152. # define TMN 18
  153. # define NMT 6
  154. #else
  155.     /* MP3 values */
  156. # define TMN 29
  157. # define NMT 6
  158. #endif
  159.  
  160.  
  161.  
  162.  
  163.  
  164.   if(gfc->psymodel_init==0) {
  165.     FLOAT8    SNR_s[CBANDS];
  166.     gfc->psymodel_init=1;
  167.  
  168.  
  169.     samplerate = gfp->out_samplerate;
  170.     switch(gfp->out_samplerate){
  171.     case 32000: break;
  172.     case 44100: break;
  173.     case 48000: break;
  174.     case 16000: break;
  175.     case 22050: break;
  176.     case 24000: break;
  177.     case  8000: samplerate *= 2; break;  /* kludge so mpeg2.5 uses mpeg2 tables  for now */
  178.     case 11025: samplerate *= 2; break;
  179.     case 12000: samplerate *= 2; break;
  180.     default:    ERRORF("error, invalid sampling frequency: %d Hz\n",
  181.             gfp->out_samplerate);
  182.     return -1;
  183.     }
  184.  
  185.     gfc->ms_ener_ratio_old=.25;
  186.     gfc->blocktype_old[0]=STOP_TYPE;
  187.     gfc->blocktype_old[1]=STOP_TYPE;
  188.     gfc->blocktype_old[0]=SHORT_TYPE;
  189.     gfc->blocktype_old[1]=SHORT_TYPE;
  190.  
  191.     for (i=0; i<4; ++i) {
  192.       for (j=0; j<CBANDS; ++j) {
  193.     gfc->nb_1[i][j]=1e20;
  194.     gfc->nb_2[i][j]=1e20;
  195.       }
  196.       for ( sb = 0; sb < SBPSY_l; sb++ ) {
  197.     gfc->en[i].l[sb] = 1e20;
  198.     gfc->thm[i].l[sb] = 1e20;
  199.       }
  200.       for (j=0; j<3; ++j) {
  201.     for ( sb = 0; sb < SBPSY_s; sb++ ) {
  202.       gfc->en[i].s[sb][j] = 1e20;
  203.       gfc->thm[i].s[sb][j] = 1e20;
  204.     }
  205.       }
  206.     }
  207.  
  208.  
  209.  
  210.     
  211.     /*  gfp->cwlimit = sfreq*j/1024.0;  */
  212.     gfc->cw_lower_index=6;
  213.     if (gfp->cwlimit>0) 
  214.       cwlimit=gfp->cwlimit;
  215.     else
  216.       cwlimit=8.8717;
  217.     gfc->cw_upper_index = cwlimit*1000.0*1024.0/((FLOAT8)samplerate);
  218.     gfc->cw_upper_index=Min(HBLKSIZE-4,gfc->cw_upper_index);      /* j+3 < HBLKSIZE-1 */
  219.     gfc->cw_upper_index=Max(6,gfc->cw_upper_index);
  220.  
  221.     for ( j = 0; j < HBLKSIZE; j++ )
  222.       gfc->cw[j] = 0.4;
  223.     
  224.     
  225.  
  226.     i=L3para_read( gfp,(FLOAT8) samplerate,gfc->numlines_l,gfc->numlines_s,
  227.           gfc->minval,gfc->s3_l,gfc->s3_s,SNR_s,gfc->bu_l,
  228.           gfc->bo_l,gfc->w1_l,gfc->w2_l, gfc->bu_s,gfc->bo_s,
  229.           gfc->w1_s,gfc->w2_s,&gfc->npart_l_orig,&gfc->npart_l,
  230.           &gfc->npart_s_orig,&gfc->npart_s );
  231.     if (i!=0) return -1;
  232.     
  233.  
  234.  
  235.     /* npart_l_orig   = number of partition bands before convolution */
  236.     /* npart_l  = number of partition bands after convolution */
  237.     
  238.     for (i=0; i<gfc->npart_l; i++) {
  239.       for (j = 0; j < gfc->npart_l_orig; j++) {
  240.     if (gfc->s3_l[i][j] != 0.0)
  241.       break;
  242.       }
  243.       gfc->s3ind[i][0] = j;
  244.       
  245.       for (j = gfc->npart_l_orig - 1; j > 0; j--) {
  246.     if (gfc->s3_l[i][j] != 0.0)
  247.       break;
  248.       }
  249.       gfc->s3ind[i][1] = j;
  250.     }
  251.  
  252.  
  253.     for (i=0; i<gfc->npart_s; i++) {
  254.       for (j = 0; j < gfc->npart_s_orig; j++) {
  255.     if (gfc->s3_s[i][j] != 0.0)
  256.       break;
  257.       }
  258.       gfc->s3ind_s[i][0] = j;
  259.       
  260.       for (j = gfc->npart_s_orig - 1; j > 0; j--) {
  261.     if (gfc->s3_s[i][j] != 0.0)
  262.       break;
  263.       }
  264.       gfc->s3ind_s[i][1] = j;
  265.     }
  266.     
  267.     
  268.     /*  
  269.       #include "debugscalefac.c"
  270.     */
  271.     
  272.  
  273.  
  274.  
  275. #define rpelev 2
  276. #define rpelev2 16
  277.  
  278.     /* compute norm_l, norm_s instead of relying on table data */
  279.     for ( b = 0;b < gfc->npart_l; b++ ) {
  280.       FLOAT8 norm=0;
  281.       for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) {
  282.     norm += gfc->s3_l[b][k];
  283.       }
  284.       if (gfp->exp_nspsytune) {
  285.     for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) {
  286.       gfc->s3_l[b][k] *= 0.7;
  287.     }
  288.       } else {
  289.     for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) {
  290.       gfc->s3_l[b][k] /=  norm;
  291.     }
  292.       }
  293.       /*DEBUGF("%i  norm=%f  norm_l=%f \n",b,1/norm,norm_l[b]);*/
  294.     }
  295.  
  296.  
  297.     for ( b = 0;b < gfc->npart_s; b++ ) {
  298.       FLOAT8 norm=0;
  299.       for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
  300.     norm += gfc->s3_s[b][k];
  301.       }
  302.       if (gfp->exp_nspsytune) {
  303.     for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
  304.       gfc->s3_s[b][k] *= 0.7;
  305.     }
  306.       } else {
  307.     for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
  308.       gfc->s3_s[b][k] *= SNR_s[b] / norm;
  309.     }
  310.       }
  311.       /*DEBUGF("%i  norm=%f  norm_s=%f \n",b,1/norm,norm_l[b]);*/
  312.     }
  313.     
  314.     init_fft(gfp);
  315.   }
  316.   /************************* End of Initialization *****************************/
  317.   
  318.  
  319.  
  320.   
  321.   
  322.   numchn = gfc->stereo;
  323.   /* chn=2 and 3 = Mid and Side channels */
  324.   if (gfp->mode == MPG_MD_JOINT_STEREO) numchn=4;
  325.   for (chn=0; chn<numchn; chn++) {
  326.  
  327.     /* there is a one granule delay.  Copy maskings computed last call
  328.      * into masking_ratio to return to calling program.
  329.      */
  330.     if (chn<2) {    
  331.       /* LR maskings  */
  332.       percep_entropy[chn] = gfc->pe[chn]; 
  333.       masking_ratio[gr_out][chn].thm = gfc->thm[chn];
  334.       masking_ratio[gr_out][chn].en = gfc->en[chn];
  335.     }else{
  336.       /* MS maskings  */
  337.       percep_MS_entropy[chn-2] = gfc->pe[chn]; 
  338.       masking_MS_ratio[gr_out][chn-2].en = gfc->en[chn];
  339.       masking_MS_ratio[gr_out][chn-2].thm = gfc->thm[chn];
  340.     }
  341.       
  342.  
  343.     /**********************************************************************
  344.      *  compute FFTs
  345.      **********************************************************************/
  346.     wsamp_s = gfc->wsamp_S+(chn & 1);
  347.     wsamp_l = gfc->wsamp_L+(chn & 1);
  348.     if (chn<2) {    
  349.       fft_long ( gfp, *wsamp_l, chn, buffer);
  350.       fft_short( gfp, *wsamp_s, chn, buffer);
  351.     } 
  352.     /* FFT data for mid and side channel is derived from L & R */
  353.     if (chn == 2)
  354.       {
  355.         for (j = BLKSIZE-1; j >=0 ; --j)
  356.         {
  357.           FLOAT l = gfc->wsamp_L[0][j];
  358.           FLOAT r = gfc->wsamp_L[1][j];
  359.           gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
  360.           gfc->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 = gfc->wsamp_S[0][b][j];
  367.             FLOAT r = gfc->wsamp_S[1][b][j];
  368.             gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
  369.             gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
  370.           }
  371.         }
  372.       }
  373.  
  374.  
  375.     /**********************************************************************
  376.      *  compute energies
  377.      **********************************************************************/
  378.     
  379.     
  380.     
  381.     gfc->energy[0]  = (*wsamp_l)[0];
  382.     gfc->energy[0] *= gfc->energy[0];
  383.     
  384.     tot_ener[chn] = gfc->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.       gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
  391.  
  392.       if (BLKSIZE/2-j > 10)
  393.     tot_ener[chn] += gfc->energy[BLKSIZE/2-j];
  394.     }
  395.     for (b = 2; b >= 0; --b)
  396.     {
  397.       gfc->energy_s[b][0]  = (*wsamp_s)[b][0];
  398.       gfc->energy_s[b][0] *=  gfc->energy_s [b][0];
  399.       for (j=BLKSIZE_s/2-1; j >= 0; --j)
  400.       {
  401.         FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
  402.         FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
  403.         gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
  404.       }
  405.     }
  406.  
  407.  
  408.   if(gfp->gtkflag) {
  409.     for (j=0; j<HBLKSIZE ; j++) {
  410.       gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j];
  411.       gfc->energy_save[chn][j]=gfc->energy[j];
  412.     }
  413.   }
  414.     
  415.     /**********************************************************************
  416.      *    compute unpredicatability of first six spectral lines            * 
  417.      **********************************************************************/
  418.     for ( j = 0; j < gfc->cw_lower_index; j++ )
  419.       {     /* calculate unpredictability measure cw */
  420.     FLOAT an, a1, a2;
  421.     FLOAT bn, b1, b2;
  422.     FLOAT rn, r1, r2;
  423.     FLOAT numre, numim, den;
  424.  
  425.     a2 = gfc-> ax_sav[chn][1][j];
  426.     b2 = gfc-> bx_sav[chn][1][j];
  427.     r2 = gfc-> rx_sav[chn][1][j];
  428.     a1 = gfc-> ax_sav[chn][1][j] = gfc-> ax_sav[chn][0][j];
  429.     b1 = gfc-> bx_sav[chn][1][j] = gfc-> bx_sav[chn][0][j];
  430.     r1 = gfc-> rx_sav[chn][1][j] = gfc-> rx_sav[chn][0][j];
  431.     an = gfc-> ax_sav[chn][0][j] = (*wsamp_l)[j];
  432.     bn = gfc-> bx_sav[chn][0][j] = j==0 ? (*wsamp_l)[0] : (*wsamp_l)[BLKSIZE-j];  
  433.     rn = gfc-> rx_sav[chn][0][j] = sqrt(gfc->energy[j]);
  434.  
  435.     { /* square (x1,y1) */
  436.       if( r1 != 0 ) {
  437.         numre = (a1*b1);
  438.         numim = (a1*a1-b1*b1)*(FLOAT)0.5;
  439.         den = r1*r1;
  440.       } else {
  441.         numre = 1;
  442.         numim = 0;
  443.         den = 1;
  444.       }
  445.     }
  446.     
  447.     { /* multiply by (x2,-y2) */
  448.       if( r2 != 0 ) {
  449.         FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
  450.         FLOAT tmp1 = -a2*numre+tmp2;
  451.         numre =       -b2*numim+tmp2;
  452.         numim = tmp1;
  453.         den *= r2;
  454.       } else {
  455.         /* do nothing */
  456.       }
  457.     }
  458.     
  459.     { /* r-prime factor */
  460.       FLOAT tmp = (2*r1-r2)/den;
  461.       numre *= tmp;
  462.       numim *= tmp;
  463.     }
  464.     den=rn+fabs(2*r1-r2);
  465.     if( den != 0 ) {
  466.       numre = (an+bn)*(FLOAT)0.5-numre;
  467.       numim = (an-bn)*(FLOAT)0.5-numim;
  468.       den = sqrt(numre*numre+numim*numim)/den;
  469.     }
  470.     gfc->cw[j] = den;
  471.       }
  472.  
  473.  
  474.  
  475.     /**********************************************************************
  476.      *     compute unpredicatibility of next 200 spectral lines            *
  477.      **********************************************************************/ 
  478.     for ( j = gfc->cw_lower_index; j < gfc->cw_upper_index; j += 4 )
  479.       {/* calculate unpredictability measure cw */
  480.     FLOAT rn, r1, r2;
  481.     FLOAT numre, numim, den;
  482.     
  483.     k = (j+2) / 4; 
  484.     
  485.     { /* square (x1,y1) */
  486.       r1 = gfc->energy_s[0][k];
  487.       if( r1 != 0 ) {
  488.         FLOAT a1 = (*wsamp_s)[0][k]; 
  489.         FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */
  490.         numre = (a1*b1);
  491.         numim = (a1*a1-b1*b1)*(FLOAT)0.5;
  492.         den = r1;
  493.         r1 = sqrt(r1);
  494.       } else {
  495.         numre = 1;
  496.         numim = 0;
  497.         den = 1;
  498.       }
  499.     }
  500.     
  501.     
  502.     { /* multiply by (x2,-y2) */
  503.       r2 = gfc->energy_s[2][k];
  504.       if( r2 != 0 ) {
  505.         FLOAT a2 = (*wsamp_s)[2][k]; 
  506.         FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k];
  507.         
  508.         
  509.         FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
  510.         FLOAT tmp1 = -a2*numre+tmp2;
  511.         numre =       -b2*numim+tmp2;
  512.         numim = tmp1;
  513.         
  514.         r2 = sqrt(r2);
  515.         den *= r2;
  516.       } else {
  517.         /* do nothing */
  518.       }
  519.     }
  520.     
  521.     { /* r-prime factor */
  522.       FLOAT tmp = (2*r1-r2)/den;
  523.       numre *= tmp;
  524.       numim *= tmp;
  525.     }
  526.     
  527.     rn = sqrt(gfc->energy_s[1][k]);
  528.     den=rn+fabs(2*r1-r2);
  529.     if( den != 0 ) {
  530.       FLOAT an = (*wsamp_s)[1][k]; 
  531.       FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k];
  532.       numre = (an+bn)*(FLOAT)0.5-numre;
  533.       numim = (an-bn)*(FLOAT)0.5-numim;
  534.       den = sqrt(numre*numre+numim*numim)/den;
  535.     }
  536.     
  537.     gfc->cw[j+1] = gfc->cw[j+2] = gfc->cw[j+3] = gfc->cw[j] = den;
  538.       }
  539.     
  540. #if 0
  541.     for ( j = 14; j < HBLKSIZE-4; j += 4 )
  542.       {/* calculate energy from short ffts */
  543.     FLOAT8 tot,ave;
  544.     k = (j+2) / 4; 
  545.     for (tot=0, sblock=0; sblock < 3; sblock++)
  546.       tot+=gfc->energy_s[sblock][k];
  547.     ave = gfc->energy[j+1]+ gfc->energy[j+2]+ gfc->energy[j+3]+ gfc->energy[j];
  548.     ave /= 4.;
  549.     /*
  550.       DEBUGF("energy / tot %i %5.2f   %e  %e\n",j,ave/(tot*16./3.),
  551.       ave,tot*16./3.);
  552.     */
  553.     gfc->energy[j+1] = gfc->energy[j+2] = gfc->energy[j+3] =  gfc->energy[j]=tot;
  554.       }
  555. #endif
  556.     
  557.     /**********************************************************************
  558.      *    Calculate the energy and the unpredictability in the threshold   *
  559.      *    calculation partitions                                           *
  560.      **********************************************************************/
  561.  
  562.     if (gfp->exp_nspsytune) {
  563.       b = 0;
  564.       for (j = 0; j < gfc->cw_upper_index && gfc->numlines_l[b] && b<gfc->npart_l_orig; )
  565.         {
  566.       FLOAT8 ebb, cbb,m,a;
  567.   
  568.       ebb = gfc->energy[j];
  569.       m = a = gfc->energy[j];
  570.       cbb = gfc->energy[j] * gfc->cw[j];
  571.       j++;
  572.   
  573.       for (i = gfc->numlines_l[b] - 1; i > 0; i--)
  574.         {
  575.           FLOAT8 el = gfc->energy[j];
  576.           ebb += gfc->energy[j];
  577.           cbb += gfc->energy[j] * gfc->cw[j];
  578.           a += el;
  579.           m = m < el ? el : m;
  580.           j++;
  581.         }
  582.       eb[b] = ebb;
  583.       cb[b] = cbb;
  584.       max[b] = m; avg[b] = a;
  585.       b++;
  586.         }
  587.   
  588.       for (; b < gfc->npart_l_orig; b++ )
  589.         {
  590.       FLOAT8 el = gfc->energy[j];
  591.       FLOAT8 m,a;
  592.       FLOAT8 ebb = gfc->energy[j++];
  593.       m = a = el;
  594.       assert(gfc->numlines_l[b]);
  595.       for (i = gfc->numlines_l[b] - 1; i > 0; i--)
  596.         {
  597.           el = gfc->energy[j];
  598.           ebb += gfc->energy[j++];
  599.           a += el;
  600.           m = m < el ? el : m;
  601.         }
  602.       eb[b] = ebb;
  603.       cb[b] = ebb * 0.4;
  604.  
  605.       max[b] = m; avg[b] = a;
  606.         }
  607.   
  608.       j = 0;
  609.       for (b=0; b < gfc->npart_l_orig; b++ )
  610.     {
  611.       int c1,c2;
  612.       FLOAT8 m,a;
  613.       tonality2[b] = 0;
  614.       c1 = c2 = 0;
  615.       m = a = 0;
  616.       for(k=b;k<=b;k++)
  617.         {
  618.           if (k >= 0 && k < gfc->npart_l_orig) {
  619.         c1++;
  620.         c2 += gfc->numlines_l[k];
  621.         a += avg[k];
  622.         m = m < max[k] ? max[k] : m;
  623.           }
  624.         }
  625.  
  626.       tonality2[b] = a == 0 ? 0 : (m / a)*c1;
  627.      }
  628.  
  629.       for (b=0; b < gfc->npart_l_orig; b++ )
  630.     {
  631.       FLOAT8 f=1;
  632.       eb2[b] = eb[b];
  633.  
  634.       if (tonality2[b]*b > 19) {
  635.         if (tonality2[b]*b > 21) f = 0.3;
  636.         else f = 1-(tonality2[b]*b-19)/2*0.7;
  637.       }
  638.       eb2[b] *= f;
  639.     }
  640.     } else {
  641.       b = 0;
  642.       for (j = 0; j < gfc->cw_upper_index && gfc->numlines_l[b] && b<gfc->npart_l_orig; )
  643.     {
  644.       FLOAT8 ebb, cbb;
  645.  
  646.       ebb = gfc->energy[j];
  647.       cbb = gfc->energy[j] * gfc->cw[j];
  648.       j++;
  649.  
  650.  
  651.       for (i = gfc->numlines_l[b] - 1; i > 0; i--)
  652.         {
  653.           ebb += gfc->energy[j];
  654.           cbb += gfc->energy[j] * gfc->cw[j];
  655.           j++;
  656.         }
  657.       eb[b] = ebb;
  658.       cb[b] = cbb;
  659.       b++;
  660.     }
  661.  
  662.       for (; b < gfc->npart_l_orig; b++ )
  663.     {
  664.       FLOAT8 ebb = gfc->energy[j++];
  665.       assert(gfc->numlines_l[b]);
  666.       for (i = gfc->numlines_l[b] - 1; i > 0; i--)
  667.         {
  668.           ebb += gfc->energy[j++];
  669.         }
  670.       eb[b] = ebb;
  671.       cb[b] = ebb * 0.4;
  672.     }
  673.     }
  674.  
  675.     /**********************************************************************
  676.      *      convolve the partitioned energy and unpredictability           *
  677.      *      with the spreading function, s3_l[b][k]                        *
  678.      ******************************************************************** */
  679.     gfc->pe[chn] = 0;        /*  calculate percetual entropy */
  680.     for ( b = 0;b < gfc->npart_l; b++ )
  681.       {
  682.     FLOAT8 tbb,ecb,ctb;
  683.  
  684.     ecb = 0;
  685.     ctb = 0;
  686.  
  687.     if (gfp->exp_nspsytune) {
  688.       for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
  689.         {
  690.           ecb = mask_add(ecb,gfc->s3_l[b][k] * eb2[k],k-b);
  691.           ctb += gfc->s3_l[b][k] * cb[k];
  692.         }
  693.     } else {
  694.       for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
  695.         {
  696.           ecb += gfc->s3_l[b][k] * eb[k];    /* sprdngf for Layer III */
  697.           ctb += gfc->s3_l[b][k] * cb[k];
  698.         }
  699.     }
  700.  
  701.     /* calculate the tonality of each threshold calculation partition 
  702.      * calculate the SNR in each threshhold calculation partition 
  703.      * tonality = -0.299 - .43*log(ctb/ecb);
  704.      * tonality = 0:           use NMT   (lots of masking)
  705.      * tonality = 1:           use TMN   (little masking)
  706.      */
  707.  
  708. /* ISO values */
  709. #define CONV1 (-.299)
  710. #define CONV2 (-.43)
  711.  
  712.     tbb = ecb;
  713.     if (tbb != 0)
  714.       {
  715.         tbb = ctb / tbb;
  716.         if (tbb <= exp((1-CONV1)/CONV2)) 
  717.           { /* tonality near 1 */
  718.         tbb = exp(-LN_TO_LOG10 * TMN);
  719.           }
  720.         else if (tbb >= exp((0-CONV1)/CONV2)) 
  721.           {/* tonality near 0 */
  722.         tbb = exp(-LN_TO_LOG10 * NMT);
  723.           }
  724.         else
  725.           {
  726.         /* convert to tonality index */
  727.         /* tonality small:   tbb=1 */
  728.         /* tonality large:   tbb=-.299 */
  729.         tbb = CONV1 + CONV2*log(tbb);
  730.         tbb = exp(-LN_TO_LOG10 * ( TMN*tbb + (1-tbb)*NMT) );
  731.           }
  732.       }
  733.     /* at this point, tbb represents the amount the spreading function
  734.      * will be reduced.  The smaller the value, the less masking.
  735.      * minval[] = 1 (0db)     says just use tbb.
  736.      * minval[]= .01 (-20db)  says reduce spreading function by 
  737.      *                        at least 20db.  
  738.      */
  739.  
  740.     if (gfp->exp_nspsytune) {
  741.       tbb = exp(-LN_TO_LOG10 * NMT);
  742.       ecb *= tbb;
  743.       ecb *= pow(10,4/20);   /* tuned by hearing tests */
  744.     } else {
  745.       tbb = Min(gfc->minval[b], tbb);
  746.       ecb *= tbb;
  747.     }
  748.  
  749.     /* long block pre-echo control.   */
  750.     /* rpelev=2.0, rpelev2=16.0 */
  751.     /* note: all surges in PE are because of this pre-echo formula
  752.      * for thr[b].  If it this is not used, PE is always around 600
  753.      */
  754.     /* dont use long block pre-echo control if previous granule was 
  755.      * a short block.  This is to avoid the situation:   
  756.      * frame0:  quiet (very low masking)  
  757.      * frame1:  surge  (triggers short blocks)
  758.      * frame2:  regular frame.  looks like pre-echo when compared to 
  759.      *          frame0, but all pre-echo was in frame1.
  760.      */
  761.     //    ecb = Max(ecb,gfc->ATH_partitionbands[b]);
  762.     if (gfc->blocktype_old[chn] == SHORT_TYPE )
  763.       thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */
  764.     else
  765.       thr[b] = Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b]) );
  766.  
  767.     gfc->nb_2[chn][b] = gfc->nb_1[chn][b];
  768.     gfc->nb_1[chn][b] = ecb;
  769.  
  770.     {
  771.       FLOAT8 thrpe;
  772.       thrpe = Max(thr[b],gfc->ATH_partitionbands[b]);
  773.       /*
  774.         printf("%i thr=%e   ATH=%e  \n",b,thr[b],gfc->ATH_partitionbands[b]);
  775.       */
  776.       if (thrpe < eb[b])
  777.         gfc->pe[chn] -= gfc->numlines_l[b] * log(thrpe / eb[b]);
  778.     }
  779.  
  780. #undef PRINTMASKING
  781. #ifdef PRINTMASKING
  782.      if (b != 0)
  783.       {
  784.         int i,t,e,a;
  785.  
  786.         t = tonality2[b]*b*50-950;
  787.         if (t<0) t = 0;
  788.         if (t>99) t = 99;
  789.         PRINTF1("%2d ",t);
  790.  
  791.         t = -82+20*log10(ecb); // -82+20*log10(thr[b]);
  792.         e = -82+20*log10(eb[b]);
  793.         a = -82+20*log10(gfc->ATH_partitionbands[b])+120;
  794.         for(i=0;i<a;i++) PRINTF1("a");
  795.         for(;i<t;i++) {
  796.           if (i >= e) PRINTF1("O");
  797.           else PRINTF1("@");
  798.         }
  799.         for(;i<e;i++) PRINTF1("-");
  800.         PRINTF1("\n");
  801.       }
  802. #endif
  803.       }
  804. #ifdef PRINTMASKING
  805.     PRINTF1("\n");
  806. #endif
  807.     
  808.     /*************************************************************** 
  809.      * determine the block type (window type) based on L & R channels
  810.      * 
  811.      ***************************************************************/
  812.     {  /* compute PE for all 4 channels */
  813.       FLOAT mn,mx,ma=0,mb=0,mc=0;
  814.       
  815.       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
  816.     {
  817.       ma += gfc->energy_s[0][j];
  818.       mb += gfc->energy_s[1][j];
  819.       mc += gfc->energy_s[2][j];
  820.     }
  821.       mn = Min(ma,mb);
  822.       mn = Min(mn,mc);
  823.       mx = Max(ma,mb);
  824.       mx = Max(mx,mc);
  825.       
  826.       /* bit allocation is based on pe.  */
  827.       if (mx>mn) {
  828.     FLOAT8 tmp = 400*log(mx/(1e-12+mn));
  829.     if (tmp>gfc->pe[chn]) gfc->pe[chn]=tmp;
  830.       }
  831.  
  832.       /* block type is based just on L or R channel */      
  833.       if (chn<2) {
  834.     uselongblock[chn] = 1;
  835.     
  836.     /* tuned for t1.wav.  doesnt effect most other samples */
  837.     if (gfc->pe[chn] > 3000) 
  838.       uselongblock[chn]=0;
  839.     
  840.     if ( mx > 30*mn ) 
  841.       {/* big surge of energy - always use short blocks */
  842.         uselongblock[chn] = 0;
  843.       } 
  844.     else if ((mx > 10*mn) && (gfc->pe[chn] > 1000))
  845.       {/* medium surge, medium pe - use short blocks */
  846.         uselongblock[chn] = 0;
  847.       }
  848.     
  849.     /* disable short blocks */
  850.     if (gfp->no_short_blocks)
  851.       uselongblock[chn]=1;
  852.       }
  853.     }
  854.  
  855.  
  856.     if (gfp->gtkflag) {
  857.       FLOAT mn,mx,ma=0,mb=0,mc=0;
  858.  
  859.       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
  860.       {
  861.         ma += gfc->energy_s[0][j];
  862.         mb += gfc->energy_s[1][j];
  863.         mc += gfc->energy_s[2][j];
  864.       }
  865.       mn = Min(ma,mb);
  866.       mn = Min(mn,mc);
  867.       mx = Max(ma,mb);
  868.       mx = Max(mx,mc);
  869.  
  870.       gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn];
  871.       gfc->ers_save[chn]=(mx/(1e-12+mn));
  872.       gfc->pinfo->pe[gr_out][chn]=gfc->pe_save[chn];
  873.       gfc->pe_save[chn]=gfc->pe[chn];
  874.     }
  875.  
  876.  
  877.     /*************************************************************** 
  878.      * compute masking thresholds for both short and long blocks
  879.      ***************************************************************/
  880.     /* longblock threshold calculation (part 2) */
  881.     for ( sb = 0; sb < SBPSY_l; sb++ )
  882.       {
  883.     FLOAT8 enn,thmm;
  884.     if (!gfp->exp_nspsytune) {
  885.       /* additive masking */
  886.       enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
  887.       thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]];
  888.       for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
  889.         {
  890.           enn  += eb[b];
  891.           thmm += thr[b];
  892.         }
  893.     } else {
  894.       enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
  895.       thmm = Min(thr[gfc->bu_l[sb]],thr[gfc->bo_l[sb]]);
  896.       for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
  897.         {
  898.           enn  += eb[b];
  899.           thmm = Min(thr[b],thmm);
  900.         }
  901.       thmm*=(1+gfc->bo_l[sb]-gfc->bu_l[sb]);
  902.     }
  903.  
  904.     gfc->en[chn].l[sb] = enn;
  905.     gfc->thm[chn].l[sb] = thmm;
  906.       }
  907.     
  908.  
  909.     /* threshold calculation for short blocks */
  910.     for ( sblock = 0; sblock < 3; sblock++ )
  911.       {
  912.     j = 0;
  913.     for ( b = 0; b < gfc->npart_s_orig; b++ )
  914.       {
  915.         FLOAT ecb = gfc->energy_s[sblock][j++];
  916.         for (i = 1 ; i<gfc->numlines_s[b]; ++i)
  917.           {
  918.         ecb += gfc->energy_s[sblock][j++];
  919.           }
  920.         eb[b] = ecb;
  921.       }
  922.  
  923.     for ( b = 0; b < gfc->npart_s; b++ )
  924.       {
  925.         FLOAT8 ecb = 0;
  926.         for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ )
  927.           {
  928.         ecb += gfc->s3_s[b][k] * eb[k];
  929.           }
  930.         thr[b] = Max (1e-6, ecb);
  931.       }
  932.  
  933.     for ( sb = 0; sb < SBPSY_s; sb++ )
  934.       {
  935.         FLOAT8 enn  = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]];
  936.         FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]];
  937.         for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ )
  938.           {
  939.         enn  += eb[b];
  940.         thmm += thr[b];
  941.           }
  942.         gfc->en[chn].s[sb][sblock] = enn;
  943.         gfc->thm[chn].s[sb][sblock] = thmm;
  944.       }
  945.       }
  946.  
  947.   } /* end loop over chn */
  948.  
  949.  
  950.   /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */
  951.   if ( numchn==4 /* mid/side and r/l */) {
  952.     FLOAT8 rside,rmid,mld;
  953.     int chmid=2,chside=3; 
  954.     
  955.     for ( sb = 0; sb < SBPSY_l; sb++ ) {
  956.       /* use this fix if L & R masking differs by 2db or less */
  957.       /* if db = 10*log10(x2/x1) < 2 */
  958.       /* if (x2 < 1.58*x1) { */
  959.       if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb]
  960.       && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) {
  961.  
  962.     mld = gfc->mld_l[sb]*gfc->en[chside].l[sb];
  963.     rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld));
  964.  
  965.     mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb];
  966.     rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld));
  967.  
  968.     gfc->thm[chmid].l[sb]=rmid;
  969.     gfc->thm[chside].l[sb]=rside;
  970.       }
  971.     }
  972.     for ( sb = 0; sb < SBPSY_s; sb++ ) {
  973.       for ( sblock = 0; sblock < 3; sblock++ ) {
  974.     if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock]
  975.         && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) {
  976.  
  977.       mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock];
  978.       rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld));
  979.  
  980.       mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock];
  981.       rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld));
  982.  
  983.       gfc->thm[chmid].s[sb][sblock]=rmid;
  984.       gfc->thm[chside].s[sb][sblock]=rside;
  985.     }
  986.       }
  987.     }
  988.   }
  989.  
  990.  
  991.   
  992.  
  993.   
  994.   
  995.   if (gfp->mode == MPG_MD_JOINT_STEREO)  {
  996.     /* determin ms_ratio from masking thresholds*/
  997.     /* use ms_stereo (ms_ratio < .35) if average thresh. diff < 5 db */
  998.     FLOAT8 db,x1,x2,sidetot=0,tot=0;
  999.     for (sb= SBPSY_l/4 ; sb< SBPSY_l; sb ++ ) {
  1000.       x1 = Min(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
  1001.       x2 = Max(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
  1002.       /* thresholds difference in db */
  1003.       if (x2 >= 1000*x1)  db=3;
  1004.       else db = log10(x2/x1);  
  1005.       /*  DEBUGF("db = %f %e %e  \n",db,gfc->thm[0].l[sb],gfc->thm[1].l[sb]);*/
  1006.       sidetot += db;
  1007.       tot++;
  1008.     }
  1009.     ms_ratio_l= (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
  1010.     ms_ratio_l = Min(ms_ratio_l,0.5);
  1011.     
  1012.     sidetot=0; tot=0;
  1013.     for ( sblock = 0; sblock < 3; sblock++ )
  1014.       for ( sb = SBPSY_s/4; sb < SBPSY_s; sb++ ) {
  1015.     x1 = Min(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
  1016.     x2 = Max(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
  1017.     /* thresholds difference in db */
  1018.     if (x2 >= 1000*x1)  db=3;
  1019.     else db = log10(x2/x1);  
  1020.     sidetot += db;
  1021.     tot++;
  1022.       }
  1023.     ms_ratio_s = (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
  1024.     ms_ratio_s = Min(ms_ratio_s,.5);
  1025.   }
  1026.  
  1027.   /*************************************************************** 
  1028.    * determin final block type
  1029.    ***************************************************************/
  1030.  
  1031.   for (chn=0; chn<gfc->stereo; chn++) {
  1032.     blocktype[chn] = NORM_TYPE;
  1033.   }
  1034.  
  1035.  
  1036.   if (gfc->stereo==2) {
  1037.     if (!gfp->allow_diff_short || gfp->mode==MPG_MD_JOINT_STEREO) {
  1038.       /* force both channels to use the same block type */
  1039.       /* this is necessary if the frame is to be encoded in ms_stereo.  */
  1040.       /* But even without ms_stereo, FhG  does this */
  1041.       int bothlong= (uselongblock[0] && uselongblock[1]);
  1042.       if (!bothlong) {
  1043.     uselongblock[0]=0;
  1044.     uselongblock[1]=0;
  1045.       }
  1046.     }
  1047.   }
  1048.  
  1049.   
  1050.   
  1051.   /* update the blocktype of the previous granule, since it depends on what
  1052.    * happend in this granule */
  1053.   for (chn=0; chn<gfc->stereo; chn++) {
  1054.     if ( uselongblock[chn])
  1055.       {                /* no attack : use long blocks */
  1056.     switch( gfc->blocktype_old[chn] ) 
  1057.       {
  1058.       case NORM_TYPE:
  1059.       case STOP_TYPE:
  1060.         blocktype[chn] = NORM_TYPE;
  1061.         break;
  1062.       case SHORT_TYPE:
  1063.         blocktype[chn] = STOP_TYPE; 
  1064.         break;
  1065.       case START_TYPE:
  1066.         ERRORF( "Error in block selecting\n" );
  1067.         LAME_ERROR_EXIT();
  1068.         break; /* problem */
  1069.       }
  1070.       } else   {
  1071.     /* attack : use short blocks */
  1072.     blocktype[chn] = SHORT_TYPE;
  1073.     if ( gfc->blocktype_old[chn] == NORM_TYPE ) {
  1074.       gfc->blocktype_old[chn] = START_TYPE;
  1075.     }
  1076.     if ( gfc->blocktype_old[chn] == STOP_TYPE ) {
  1077.       gfc->blocktype_old[chn] = SHORT_TYPE ;
  1078.     }
  1079.       }
  1080.     
  1081.     blocktype_d[chn] = gfc->blocktype_old[chn];  /* value returned to calling program */
  1082.     gfc->blocktype_old[chn] = blocktype[chn];    /* save for next call to l3psy_anal */
  1083.   }
  1084.   
  1085.   if (blocktype_d[0]==2) 
  1086.     *ms_ratio = gfc->ms_ratio_s_old;
  1087.   else
  1088.     *ms_ratio = gfc->ms_ratio_l_old;
  1089.  
  1090.   gfc->ms_ratio_s_old = ms_ratio_s;
  1091.   gfc->ms_ratio_l_old = ms_ratio_l;
  1092.  
  1093.   /* we dont know the block type of this frame yet - assume long */
  1094.   *ms_ratio_next = ms_ratio_l;
  1095.  
  1096.  
  1097.  
  1098.   /*********************************************************************/
  1099.   /* compute side_energy / (side+mid)_energy */
  1100.   /* 0 = no energy in side channel */
  1101.   /* .5 = half of total energy in side channel */
  1102.   /*********************************************************************/
  1103.   if (numchn==4)  {
  1104.     FLOAT tmp = tot_ener[3]+tot_ener[2];
  1105.     *ms_ener_ratio = gfc->ms_ener_ratio_old;
  1106.     gfc->ms_ener_ratio_old=0;
  1107.     if (tmp>0) gfc->ms_ener_ratio_old=tot_ener[3]/tmp;
  1108.   } else
  1109.     /* we didn't compute ms_ener_ratios */
  1110.     *ms_ener_ratio = 0;
  1111.  
  1112.   return 0;
  1113. }
  1114.  
  1115.  
  1116.  
  1117.  
  1118.  
  1119. int L3para_read(lame_global_flags *gfp, FLOAT8 sfreq, int *numlines_l,int *numlines_s, 
  1120. FLOAT8 *minval,
  1121. FLOAT8 s3_l[CBANDS][CBANDS], FLOAT8 s3_s[CBANDS][CBANDS],
  1122. FLOAT8 *SNR, 
  1123. int *bu_l, int *bo_l, FLOAT8 *w1_l, FLOAT8 *w2_l, 
  1124. int *bu_s, int *bo_s, FLOAT8 *w1_s, FLOAT8 *w2_s,
  1125. int *npart_l_orig,int *npart_l,int *npart_s_orig,int *npart_s)
  1126. {
  1127.   lame_internal_flags *gfc=gfp->internal_flags;
  1128.   FLOAT8 freq_tp;
  1129.   FLOAT8 bval_l[CBANDS], bval_s[CBANDS];
  1130.   int   cbmax=0, cbmax_tp;
  1131.   FLOAT8 *p = psy_data;
  1132.   int  sbmax ;
  1133.   int  i,j,k2,loop;
  1134.   int freq_scale=1;
  1135.   int partition[HBLKSIZE];
  1136.  
  1137.   /******************************************************************/
  1138.   /* Read long block data */
  1139.   /******************************************************************/
  1140.   for(loop=0;loop<6;loop++)
  1141.     {
  1142.       freq_tp = *p++;
  1143.       cbmax_tp = (int) *p++;
  1144.       cbmax_tp++;
  1145.  
  1146.       if (sfreq == freq_tp/freq_scale )
  1147.     {
  1148.       cbmax = cbmax_tp;
  1149.       for(i=0,k2=0;i<cbmax_tp;i++)
  1150.         {
  1151.           j = (int) *p++;
  1152.           numlines_l[i] = (int) *p++;
  1153.           minval[i] = exp(-((*p++) ) * LN_TO_LOG10);
  1154.           /* qthr_l[i] = *p++ */ p++;
  1155.           /* norm_l[i] = *p++*/ p++;
  1156.           /* bval_l[i] = *p++; */ p++;
  1157.           if (j!=i)
  1158.         {
  1159.           ERRORF("1. please check \"psy_data\"");
  1160.           return -1;
  1161.         }
  1162.         }
  1163.     }
  1164.       else
  1165.     p += cbmax_tp * 6;
  1166.     }
  1167.  
  1168.   *npart_l_orig = cbmax;
  1169.  
  1170.   /* Read short block data */
  1171.   for(loop=0;loop<6;loop++)
  1172.     {
  1173.       freq_tp = *p++;
  1174.       cbmax_tp = (int) *p++;
  1175.       cbmax_tp++;
  1176.       if (sfreq == freq_tp/freq_scale )
  1177.     {
  1178.       cbmax = cbmax_tp;
  1179.       for(i=0,k2=0;i<cbmax_tp;i++)
  1180.         {
  1181.           j = (int) *p++;
  1182.           numlines_s[i] = (int) *p++;
  1183.           /* qthr_s[i] = *p++*/  p++;         
  1184.           /* norm_s[i] =*p++ */ p++;         
  1185.           SNR[i] = *p++;            
  1186.           /* bval_s[i] = *p++ */ p++;
  1187.           if (j!=i)
  1188.         {
  1189.           ERRORF("3. please check \"psy_data\"");
  1190.           return -1;
  1191.         }
  1192.         }
  1193.     }
  1194.       else
  1195.     p += cbmax_tp * 6;
  1196.     }
  1197.   *npart_s_orig = cbmax;
  1198.  
  1199.   /* MPEG1 SNR_s data is given in db, convert to energy */
  1200.   if (gfp->version == 1) {
  1201.     for ( i = 0;i < *npart_s_orig; i++ ) {
  1202.       SNR[i]=exp( (FLOAT8) SNR[i] * LN_TO_LOG10 );
  1203.     }
  1204.   }
  1205.  
  1206.  
  1207.  
  1208.  
  1209.   /* Read long block data for converting threshold calculation 
  1210.      partitions to scale factor bands */
  1211.  
  1212.   for(loop=0;loop<6;loop++)
  1213.     {
  1214.       freq_tp = *p++;
  1215.       sbmax =  (int) *p++;
  1216.       sbmax++;
  1217.  
  1218.       if (sfreq == freq_tp/freq_scale)
  1219.     {
  1220.       for(i=0;i<sbmax;i++)
  1221.         {
  1222.           j = (int) *p++;
  1223.           p++;             
  1224.           bu_l[i] = (int) *p++;
  1225.           bo_l[i] = (int) *p++;
  1226.           w1_l[i] = (FLOAT8) *p++;
  1227.           w2_l[i] = (FLOAT8) *p++;
  1228.           if (j!=i)
  1229.         { ERRORF("30:please check \"psy_data\"\n");
  1230.         return -1;
  1231.         }
  1232.  
  1233.           if (i!=0)
  1234.         if ( (fabs(1.0-w1_l[i]-w2_l[i-1]) > 0.01 ) )
  1235.           {
  1236.             ERRORF("31l: please check \"psy_data.\"\n"
  1237.                            "w1,w2: %f %f \n",w1_l[i],w2_l[i-1]);
  1238.             return -1;
  1239.           }
  1240.         }
  1241.     }
  1242.       else
  1243.     p += sbmax * 6;
  1244.     }
  1245.  
  1246.   /* Read short block data for converting threshold calculation 
  1247.      partitions to scale factor bands */
  1248.  
  1249.   for(loop=0;loop<6;loop++)
  1250.     {
  1251.       freq_tp = *p++;
  1252.       sbmax = (int) *p++;
  1253.       sbmax++;
  1254.  
  1255.       if (sfreq == freq_tp/freq_scale)
  1256.     {
  1257.       for(i=0;i<sbmax;i++)
  1258.         {
  1259.           j = (int) *p++;
  1260.           p++;
  1261.           bu_s[i] = (int) *p++;
  1262.           bo_s[i] = (int) *p++;
  1263.           w1_s[i] = *p++;
  1264.           w2_s[i] = *p++;
  1265.           if (j!=i)
  1266.         { ERRORF("30:please check \"psy_data\"\n");
  1267.         return -1;
  1268.         }
  1269.  
  1270.           if (i!=0)
  1271.         if ( (fabs(1.0-w1_s[i]-w2_s[i-1]) > 0.01 ) )
  1272.           { 
  1273.                   ERRORF("31s: please check \"psy_data.\"\n"
  1274.                          "w1,w2: %f %f \n",w1_s[i],w2_s[i-1]);
  1275.           return -1;
  1276.           }
  1277.         }
  1278.     }
  1279.       else
  1280.     p += sbmax * 6;
  1281.     }
  1282.  
  1283.   /******************************************************************/
  1284.   /* done reading table data */
  1285.   /******************************************************************/
  1286.  
  1287.  
  1288.  
  1289.  
  1290. #undef NOTABLES
  1291. #ifdef NOTABLES
  1292.   /* compute numlines */
  1293.   j=0;
  1294.   for(i=0;i<CBANDS;i++)
  1295.     {
  1296.       FLOAT8 ji, bark1,bark2,delbark=.34;
  1297.       int k,j2;
  1298.  
  1299.       j2 = j;
  1300.       j2 = Min(j2,BLKSIZE/2);
  1301.       
  1302.       do {
  1303.     /* find smallest j2 >= j so that  (bark - bark_l[i-1]) < delbark */
  1304.     ji = j;
  1305.     bark1 = freq2bark(sfreq*ji/BLKSIZE);
  1306.     
  1307.     ++j2;
  1308.     ji = j2;
  1309.     bark2  = freq2bark(sfreq*ji/BLKSIZE);
  1310.       } while ((bark2 - bark1) < delbark  && j2<=BLKSIZE/2);
  1311.  
  1312.       /*
  1313.       DEBUGF("%i old n=%i  %f old numlines:  %i   new=%i (%i,%i) (%f,%f) \n",
  1314. i,*npart_l_orig,freq,numlines_l[i],j2-j,j,j2-1,bark1,bark2);
  1315.       */
  1316.       for (k=j; k<j2; ++k)
  1317.     partition[k]=i;
  1318.       numlines_l[i]=(j2-j);
  1319.       j = j2;
  1320.       if (j > BLKSIZE/2) break;
  1321.     }
  1322.   *npart_l_orig = i;
  1323.  
  1324.   /* compute which partition bands are in which scalefactor bands */
  1325.   { int i1,i2,sfb,start,end;
  1326.     FLOAT8 freq1,freq2;
  1327.     for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
  1328.       start = gfc->scalefac_band.l[ sfb ];
  1329.       end   = gfc->scalefac_band.l[ sfb+1 ];
  1330.       freq1 = sfreq*(start-.5)/(2*576);
  1331.       freq2 = sfreq*(end-1+.5)/(2*576);
  1332.              
  1333.       i1 = floor(.5 + BLKSIZE*freq1/sfreq);
  1334.       if (i1<0) i1=0;
  1335.       i2 = floor(.5 + BLKSIZE*freq2/sfreq);
  1336.       if (i2>BLKSIZE/2) i2=BLKSIZE/2;
  1337.  
  1338.       DEBUGF("old: (%i,%i)  new: (%i,%i) %i %i \n",bu_l[sfb],bo_l[sfb],
  1339.          partition[i1],partition[i2],i1,i2);
  1340.  
  1341.       w1_l[sfb]=.5;
  1342.       w2_l[sfb]=.5;
  1343.       bu_l[sfb]=partition[i1];
  1344.       bo_l[sfb]=partition[i2];
  1345.  
  1346.     }
  1347.   }
  1348.  
  1349. #endif
  1350.  
  1351.  
  1352.   /* compute bark value and ATH of each critical band */
  1353.   j=0;
  1354.   for(i=0;i<*npart_l_orig;i++)
  1355.     {
  1356.       FLOAT8 ji, freq, mval, bark;
  1357.       int k;
  1358.  
  1359.       ji = j;
  1360.       bark = freq2bark(sfreq*ji/BLKSIZE);
  1361.  
  1362.       ji = j + (numlines_l[i]-1);
  1363.       bark = .5*(bark + freq2bark(sfreq*ji/BLKSIZE));
  1364.       bval_l[i]=bark;
  1365.       //      j += numlines_l[i];
  1366.  
  1367.  
  1368.       gfc->ATH_partitionbands[i]=1e99;
  1369.       for (k=0; k < numlines_l[i]; ++k) {
  1370.     FLOAT8 freq = sfreq*j/(1000.0*BLKSIZE);
  1371.     assert( freq < 25 );
  1372.     //    freq = Min(.1,freq);    /* ignore ATH below 100hz */
  1373.     freq= ATHformula(freq);  
  1374.     freq += -20; /* scale to FFT units */
  1375.     freq = pow( 10.0, freq/10.0 );  /* convert from db -> energy */
  1376.     freq *= numlines_l[i];
  1377.     gfc->ATH_partitionbands[i]=Min(gfc->ATH_partitionbands[i],freq);
  1378.     ++j;
  1379.       }
  1380.  
  1381.  
  1382.       /*  minval formula needs work
  1383.       ji = j;
  1384.       freq = sfreq*ji/BLKSIZE;
  1385.       mval = Max(27.0 - freq/50.0, 0.0);
  1386.       //      mval = exp(-mval * LN_TO_LOG10);
  1387.  
  1388.       DEBUGF("%2i old minval=%f  new = %f  ",
  1389.          i,-log(minval[i])/LN_TO_LOG10,mval);
  1390.       if (mval < minval[i]) 
  1391.       DEBUGF("less masking than ISO tables \n");
  1392.       else
  1393.       DEBUGF("more masking than ISO tables \n");
  1394.       */
  1395.     }
  1396.  
  1397.  
  1398.  
  1399.   /************************************************************************
  1400.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1401.    * ing function, centered at band j, for band i, store for later use    *
  1402.    ************************************************************************/
  1403.   /* i.e.: sum over j to spread into signal barkval=i  
  1404.      NOTE: i and j are used opposite as in the ISO docs */
  1405.   for(i=0;i<*npart_l_orig;i++)
  1406.     {
  1407.       FLOAT8 tempx,x,tempy,temp;
  1408.       for(j=0;j<*npart_l_orig;j++)
  1409.     {
  1410.           if (i>=j) tempx = (bval_l[i] - bval_l[j])*3.0;
  1411.       else    tempx = (bval_l[i] - bval_l[j])*1.5; 
  1412.  
  1413.       if(tempx>=0.5 && tempx<=2.5)
  1414.         {
  1415.           temp = tempx - 0.5;
  1416.           x = 8.0 * (temp*temp - 2.0 * temp);
  1417.         }
  1418.       else x = 0.0;
  1419.       tempx += 0.474;
  1420.       tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1421.  
  1422. #ifdef NEWS3
  1423.       if (j>=i) tempy = (bval_l[j] - bval_l[i])*(-15);
  1424.       else    tempy = (bval_l[j] - bval_l[i])*25;
  1425.       x=0; 
  1426. #endif
  1427.       /*
  1428.       if ((i==cbmax/2)  && (fabs(bval_l[j] - bval_l[i])) < 3) {
  1429.         DEBUGF("bark=%f   x+tempy = %f  \n",bval_l[j] - bval_l[i],x+tempy);
  1430.       }
  1431.       */
  1432.  
  1433.       if (tempy <= -60.0) s3_l[i][j] = 0.0;
  1434.       else                s3_l[i][j] = exp( (x + tempy)*LN_TO_LOG10 ); 
  1435.     }
  1436.     }
  1437.  
  1438.  
  1439.   /************************************************************************/
  1440.   /* SHORT BLOCKS */
  1441.   /************************************************************************/
  1442.  
  1443. #ifdef NOTABLES
  1444.   /* compute numlines */
  1445.   j=0;
  1446.   for(i=0;i<CBANDS;i++)
  1447.     {
  1448.       FLOAT8 ji, bark1,bark2,delbark=.34;
  1449.       int k,j2;
  1450.  
  1451.       j2 = j;
  1452.       j2 = Min(j2,BLKSIZE_s/2);
  1453.       
  1454.       do {
  1455.     /* find smallest j2 >= j so that  (bark - bark_s[i-1]) < delbark */
  1456.     ji = j;
  1457.     bark1  = freq2bark(sfreq*ji/BLKSIZE_s);
  1458.     
  1459.     ++j2;
  1460.     ji = j2;
  1461.     bark2  = freq2bark(sfreq*ji/BLKSIZE_s);
  1462.  
  1463.       } while ((bark2 - bark1) < delbark  && j2<=BLKSIZE_s/2);
  1464.  
  1465.       /*
  1466.       DEBUGF("%i old n=%i  %f old numlines:  %i   new=%i (%i,%i) (%f,%f) \n",
  1467. i,*npart_s_orig,freq,numlines_s[i],j2-j,j,j2-1,bark1,bark2);
  1468.       */
  1469.       for (k=j; k<j2; ++k)
  1470.     partition[k]=i;
  1471.       numlines_s[i]=(j2-j);
  1472.       j = j2;
  1473.       if (j > BLKSIZE_s/2) break;
  1474.     }
  1475.   *npart_s_orig = i;
  1476.  
  1477.   /* compute which partition bands are in which scalefactor bands */
  1478.   { int i1,i2,sfb,start,end;
  1479.     FLOAT8 freq1,freq2;
  1480.     for ( sfb = 0; sfb < SBMAX_s; sfb++ ) {
  1481.       start = gfc->scalefac_band.s[ sfb ];
  1482.       end   = gfc->scalefac_band.s[ sfb+1 ];
  1483.       freq1 = sfreq*(start-.5)/(2*192);
  1484.       freq2 = sfreq*(end-1+.5)/(2*192);
  1485.              
  1486.       i1 = floor(.5 + BLKSIZE_s*freq1/sfreq);
  1487.       if (i1<0) i1=0;
  1488.       i2 = floor(.5 + BLKSIZE_s*freq2/sfreq);
  1489.       if (i2>BLKSIZE_s/2) i2=BLKSIZE_s/2;
  1490.  
  1491.       DEBUGF("old: (%i,%i)  new: (%i,%i) %i %i \n",bu_s[sfb],bo_s[sfb],
  1492.          partition[i1],partition[i2],i1,i2);
  1493.  
  1494.       w1_s[sfb]=.5;
  1495.       w2_s[sfb]=.5;
  1496.       bu_s[sfb]=partition[i1];
  1497.       bo_s[sfb]=partition[i2];
  1498.  
  1499.     }
  1500.   }
  1501.  
  1502. #endif
  1503.  
  1504.  
  1505.  
  1506.   /* compute bark values of each critical band */
  1507.   j = 0;
  1508.   for(i=0;i<*npart_s_orig;i++)
  1509.     {
  1510.       FLOAT8 ji, bark, freq, snr;
  1511.       ji = j;
  1512.       bark = freq2bark(sfreq*ji/BLKSIZE_s);
  1513.  
  1514.       ji = j + numlines_s[i] -1;
  1515.       bark = .5*(bark+freq2bark(sfreq*ji/BLKSIZE_s));
  1516.       /*
  1517.       DEBUGF("%i %i bval_s = %f  %f  numlines=%i  formula=%f \n",i,j,bval_s[i],freq,numlines_s[i],bark);
  1518.       */
  1519.       bval_s[i]=bark;
  1520.       j += numlines_s[i];
  1521.  
  1522.       /* SNR formula needs work 
  1523.       ji = j;
  1524.       freq = sfreq*ji/BLKSIZE;
  1525.       snr  = .14 + freq/80000;
  1526.       DEBUGF("%2i old SNR=%f  new = %f \n ",
  1527.          i,SNR[i],snr);
  1528.       */
  1529.     }
  1530.  
  1531.  
  1532.  
  1533.  
  1534.  
  1535.  
  1536.   /************************************************************************
  1537.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1538.    * ing function, centered at band j, for band i, store for later use    *
  1539.    ************************************************************************/
  1540.   for(i=0;i<*npart_s_orig;i++)
  1541.     {
  1542.       FLOAT8 tempx,x,tempy,temp;
  1543.       for(j=0;j<*npart_s_orig;j++)
  1544.     {
  1545.           if (i>=j) tempx = (bval_s[i] - bval_s[j])*3.0;
  1546.       else    tempx = (bval_s[i] - bval_s[j])*1.5; 
  1547.  
  1548.       if(tempx>=0.5 && tempx<=2.5)
  1549.         {
  1550.           temp = tempx - 0.5;
  1551.           x = 8.0 * (temp*temp - 2.0 * temp);
  1552.         }
  1553.       else x = 0.0;
  1554.       tempx += 0.474;
  1555.       tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1556. #ifdef NEWS3
  1557.       if (j>=i) tempy = (bval_s[j] - bval_s[i])*(-15);
  1558.       else    tempy = (bval_s[j] - bval_s[i])*25;
  1559.       x=0; 
  1560. #endif
  1561.       if (tempy <= -60.0) s3_s[i][j] = 0.0;
  1562.       else                s3_s[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
  1563.     }
  1564.     }
  1565.  
  1566.  
  1567.  
  1568.  
  1569.   /* compute: */
  1570.   /* npart_l_orig   = number of partition bands before convolution */
  1571.   /* npart_l  = number of partition bands after convolution */
  1572.  
  1573.   assert(*npart_l_orig <=CBANDS);
  1574.   assert(*npart_s_orig<=CBANDS);
  1575.  
  1576.   
  1577.   *npart_l=bo_l[SBPSY_l-1]+1;
  1578.   *npart_s=bo_s[SBPSY_s-1]+1;
  1579.   
  1580.   /* if npart_l = npart_l_orig + 1, we can fix that below.  else: */
  1581.   assert(*npart_l <= *npart_l_orig+1);
  1582.   assert(*npart_s <= *npart_s_orig+1);
  1583.  
  1584.  
  1585.   /* MPEG2 tables are screwed up 
  1586.    * the mapping from paritition bands to scalefactor bands will use
  1587.    * more paritition bands than we have.  
  1588.    * So we will not compute these fictitious partition bands by reducing
  1589.    * npart_l below.  */
  1590.   if (*npart_l > *npart_l_orig) {
  1591.     *npart_l=*npart_l_orig;
  1592.     bo_l[SBPSY_l-1]=(*npart_l)-1;
  1593.     w2_l[SBPSY_l-1]=1.0;
  1594.   }
  1595.  
  1596.   if (*npart_s > *npart_s_orig) {
  1597.     *npart_s=*npart_s_orig;
  1598.     bo_s[SBPSY_s-1]=(*npart_s)-1;
  1599.     w2_s[SBPSY_s-1]=1.0;
  1600.   }
  1601.  
  1602.  
  1603.     /* setup stereo demasking thresholds */
  1604.     /* formula reverse enginerred from plot in paper */
  1605.     for ( i = 0; i < SBPSY_s; i++ ) {
  1606.       FLOAT8 arg,mld;
  1607.       arg = freq2bark(sfreq*gfc->scalefac_band.s[i]/(2*192));
  1608.       arg = (Min(arg, 15.5)/15.5);
  1609.  
  1610.       mld = 1.25*(1-cos(PI*arg))-2.5;
  1611.       gfc->mld_s[i] = pow(10.0,mld);
  1612.     }
  1613.     for ( i = 0; i < SBPSY_l; i++ ) {
  1614.       FLOAT8 arg,mld;
  1615.       arg = freq2bark(sfreq*gfc->scalefac_band.l[i]/(2*576));
  1616.       arg = (Min(arg, 15.5)/15.5);
  1617.  
  1618.       mld = 1.25*(1-cos(PI*arg))-2.5;
  1619.       gfc->mld_l[i] = pow(10.0,mld);
  1620.     }
  1621.  
  1622.   
  1623.   return 0;
  1624. }
  1625.