/********************************************************************** * ISO MPEG Audio Subgroup Software Simulation Group (1996) * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension * * $Id: l3psy.c,v 1.2 1997/01/19 22:28:29 rowlands Exp $ * * $Log: l3psy.c,v $ * Revision 1.2  1997/01/19 22:28:29  rowlands * Layer 3 bug fixes from Seymour Shlien * * Revision 1.1  1996/02/14 04:04:23  rowlands * Initial revision * * Received from Mike Coleman **********************************************************************//********************************************************************** *   date   programmers         comment                               * * 2/25/91  Davis Pan           start of version 1.0 records          * * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         * * 7/10/91  Earle Jennings      Ported to MsDos.                      * *                              replace of floats with FLOAT          * * 2/11/92  W. Joseph Carter    Fixed mem_alloc() arg for "absthr".   * * 3/16/92  Masahiro Iwadare	Modification for Layer III            * * 17/4/93  Masahiro Iwadare    Updated for IS Modification           * **********************************************************************/#include "common.h"#include "encoder.h"#include "l3psy.h"#include "l3side.h"#include <assert.h>#define maximum(x,y) ( (x>y) ? x : y )#define minimum(x,y) ( (x<y) ? x : y )void L3para_read( double sfreq, int numlines[CBANDS], int partition_l[HBLKSIZE],		  double minval[CBANDS], double qthr_l[CBANDS], double norm_l[CBANDS],		  double s3_l[CBANDS][CBANDS], int partition_s[HBLKSIZE_s], double qthr_s[CBANDS_s],		  double norm_s[CBANDS_s], double SNR_s[CBANDS_s],		  int cbw_l[SBMAX_l], int bu_l[SBMAX_l], int bo_l[SBMAX_l],		  double w1_l[SBMAX_l], double w2_l[SBMAX_l],		  int cbw_s[SBMAX_s], int bu_s[SBMAX_s], int bo_s[SBMAX_s],		  double w1_s[SBMAX_s], double w2_s[SBMAX_s] );									void L3psycho_anal( short int *buffer, short int savebuf[1344], int chn, int lay, FLOAT snr32[32],		    double sfreq, double ratio_d[21], double ratio_ds[12][3],		    double *pe, gr_info *cod_info ){    static double ratio[2][21];    static double ratio_s[2][12][3];    int blocktype;    unsigned int   b, i, j, k;    double         r_prime, phi_prime; /* not FLOAT */    FLOAT          freq_mult, bval_lo, min_thres, sum_energy;    double         tb, temp1,temp2,temp3;    /*         nint(); Layer III */    double   thr[CBANDS];/* The static variables "r", "phi_sav", "new", "old" and "oldest" have    *//* to be remembered for the unpredictability measure.  For "r" and        *//* "phi_sav", the first index from the left is the channel select and     *//* the second index is the "age" of the data.                             */   static FLOAT window_s[BLKSIZE_s] ; static int     new = 0, old = 1, oldest = 0; static int     init = 0, flush, sync_flush, syncsize, sfreq_idx; static double 	cw[HBLKSIZE], eb[CBANDS]; static double 	ctb[CBANDS]; static double	SNR_l[CBANDS], SNR_s[CBANDS_s]; static int	init_L3; static double	minval[CBANDS],qthr_l[CBANDS],norm_l[CBANDS]; static double	qthr_s[CBANDS_s],norm_s[CBANDS_s]; static double	nb_1[2][CBANDS], nb_2[2][CBANDS]; static double	s3_l[CBANDS][CBANDS]; /* s3_s[CBANDS_s][CBANDS_s]; *//* Scale Factor Bands */ static int	cbw_l[SBMAX_l],bu_l[SBMAX_l],bo_l[SBMAX_l] ; static int	cbw_s[SBMAX_s],bu_s[SBMAX_s],bo_s[SBMAX_s] ; static double	w1_l[SBMAX_l], w2_l[SBMAX_l]; static double	w1_s[SBMAX_s], w2_s[SBMAX_s]; static double	en[SBMAX_l],   thm[SBMAX_l] ; static int	blocktype_old[2] ; int	sb,sblock; static int	partition_l[HBLKSIZE],partition_s[HBLKSIZE_s];/* The following static variables are constants.                           */ static double  nmt = 5.5; static FLOAT   crit_band[27] = {0,  100,  200, 300, 400, 510, 630,  770,                               920, 1080, 1270,1480,1720,2000,2320, 2700,                              3150, 3700, 4400,5300,6400,7700,9500,12000,                             15500,25000,30000}; static FLOAT   bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,                            10.0,  7.0,  4.4,  4.5,  4.5,  4.5,  4.5,                             4.5,  4.5,  4.5,  4.5,  4.5,  4.5,  4.5,                             4.5,  4.5,  4.5,  3.5,  3.5,  3.5};/* The following pointer variables point to large areas of memory         *//* dynamically allocated by the mem_alloc() function.  Dynamic memory     *//* allocation is used in order to avoid stack frame or data area          *//* overflow errors that otherwise would have occurred at compile time     *//* on the Macintosh computer.                                             */ FLOAT          *grouped_c, *grouped_e, *nb, *cb, *ecb, *bc; FLOAT          *wsamp_r, *wsamp_i, *phi, *energy; static FLOAT	energy_s[3][256]; static FLOAT phi_s[3][256] ; /* 256 samples not 129 */ FLOAT          *c, *fthr; F32            *snrtmp; static	int	*numlines ; static int     *partition; static FLOAT   *cbval, *rnorm; static FLOAT   *window; static FLOAT   *absthr; static double  *tmn; static FCB     *s; static FHBLK   *lthr; static F2HBLK  *r, *phi_sav;/* These dynamic memory allocations simulate "automatic" variables        *//* placed on the stack.  For each mem_alloc() call here, there must be    *//* a corresponding mem_free() call at the end of this function.           */ grouped_c = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_c"); grouped_e = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_e"); nb = (FLOAT *) mem_alloc(sizeof(FCB), "nb"); cb = (FLOAT *) mem_alloc(sizeof(FCB), "cb"); ecb = (FLOAT *) mem_alloc(sizeof(FCB), "ecb"); bc = (FLOAT *) mem_alloc(sizeof(FCB), "bc"); wsamp_r = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_r"); wsamp_i = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_i"); phi = (FLOAT *) mem_alloc(sizeof(FBLK), "phi"); energy = (FLOAT *) mem_alloc(sizeof(FBLK), "energy"); c = (FLOAT *) mem_alloc(sizeof(FHBLK), "c"); fthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "fthr"); snrtmp = (F32 *) mem_alloc(sizeof(F2_32), "snrtmp");    assert( lay == 3 ); if(init==0){/* These dynamic memory allocations simulate "static" variables placed    *//* in the data space.  Each mem_alloc() call here occurs only once at     *//* initialization time.  The mem_free() function must not be called.      */     numlines = (int *) mem_alloc(sizeof(ICB), "numlines");     partition = (int *) mem_alloc(sizeof(IHBLK), "partition");     cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");     rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");     window = (FLOAT *) mem_alloc(sizeof(FBLK), "window");     absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr");      tmn = (double *) mem_alloc(sizeof(DCB), "tmn");     s = (FCB *) mem_alloc(sizeof(FCBCB), "s");     lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");     r = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "r");     phi_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "phi_sav");/*#if 0 */     i = sfreq + 0.5;     switch(i){        case 32000: sfreq_idx = 0; break;        case 44100: sfreq_idx = 1; break;        case 48000: sfreq_idx = 2; break;        default:    printf("error, invalid sampling frequency: %d Hz\n",i);        exit(-1);     }     printf("absthr[][] sampling frequency index: %d\n",sfreq_idx);     read_absthr(absthr, sfreq_idx);     switch(lay){	case 1: sync_flush=576; flush=384; syncsize=1024; break;	case 2: sync_flush=480; flush=576; syncsize=1056; break;	case 3: sync_flush=768; flush=576; syncsize=1344; break;       default: printf("Bad lay value:(%d)",lay); exit(-1); break;     }/* #endif *//* calculate HANN window coefficients *//*   for(i=0;i<BLKSIZE;i++)  window[i]  =0.5*(1-cos(2.0*PI*i/(BLKSIZE-1.0)));*/     for(i=0;i<BLKSIZE;i++)  window[i]  =0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE));     for(i=0;i<BLKSIZE_s;i++)window_s[i]=0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE_s));/* reset states used in unpredictability measure */     for(i=0;i<HBLKSIZE;i++){        r[0][0][i]=r[1][0][i]=r[0][1][i]=r[1][1][i]=0;        phi_sav[0][0][i]=phi_sav[1][0][i]=0;        phi_sav[0][1][i]=phi_sav[1][1][i]=0;        lthr[0][i] = 60802371420160.0;        lthr[1][i] = 60802371420160.0;     }/***************************************************************************** * Initialization: Compute the following constants for use later             * *    partition[HBLKSIZE] = the partition number associated with each        * *                          frequency line                                   * *    cbval[CBANDS]       = the center (average) bark value of each          * *                          partition                                        * *    numlines[CBANDS]    = the number of frequency lines in each partition  * *    tmn[CBANDS]         = tone masking noise                               * *****************************************************************************//* compute fft frequency multiplicand */     freq_mult = sfreq/BLKSIZE; /* calculate fft frequency, then bval of each line (use fthr[] as tmp storage)*/     for(i=0;i<HBLKSIZE;i++){        temp1 = i*freq_mult;        j = 1;        while(temp1>crit_band[j])j++;        fthr[i]=j-1+(temp1-crit_band[j-1])/(crit_band[j]-crit_band[j-1]);     }     partition[0] = 0;/* temp2 is the counter of the number of frequency lines in each partition */     temp2 = 1;     cbval[0]=fthr[0];     bval_lo=fthr[0];     for(i=1;i<HBLKSIZE;i++){        if((fthr[i]-bval_lo)>0.33){           partition[i]=partition[i-1]+1;           cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;           cbval[partition[i]] = fthr[i];           bval_lo = fthr[i];           numlines[partition[i-1]] = temp2;           temp2 = 1;        }        else {           partition[i]=partition[i-1];           cbval[partition[i]] += fthr[i];           temp2++;        }     }     numlines[partition[i-1]] = temp2;     cbval[partition[i-1]] = cbval[partition[i-1]]/temp2; /************************************************************************ * Now compute the spreading function, s[j][i], the value of the spread-* * ing function, centered at band j, for band i, store for later use    * ************************************************************************/     for(j=0;j<CBANDS;j++){        for(i=0;i<CBANDS;i++){           temp1 = (cbval[i] - cbval[j])*1.05;           if(temp1>=0.5 && temp1<=2.5){              temp2 = temp1 - 0.5;              temp2 = 8.0 * (temp2*temp2 - 2.0 * temp2);           }           else temp2 = 0;           temp1 += 0.474;           temp3 = 15.811389+7.5*temp1-17.5*sqrt((double) (1.0+temp1*temp1));           if(temp3 <= -100) s[i][j] = 0;           else {              temp3 = (temp2 + temp3)*LN_TO_LOG10;              s[i][j] = exp(temp3);           }        }     }  /* Calculate Tone Masking Noise values */     for(j=0;j<CBANDS;j++){        temp1 = 15.5 + cbval[j];        tmn[j] = (temp1>24.5) ? temp1 : 24.5;  /* Calculate normalization factors for the net spreading functions */        rnorm[j] = 0;        for(i=0;i<CBANDS;i++){           rnorm[j] += s[j][i];        }     }     init++; } /************************* End of Initialization *****************************/ switch(lay) {  case 1:  case 2:	for ( i=0; i<lay; i++)  {/***************************************************************************** * Net offset is 480 samples (1056-576) for layer 2; this is because one must* * stagger input data by 256 samples to synchronize psychoacoustic model with* * filter bank outputs, then stagger so that center of 1024 FFT window lines * * up with center of 576 "new" audio samples.                                * *                                                                           * * For layer 1, the input data still needs to be staggered by 256 samples,   * * then it must be staggered again so that the 384 "new" samples are centered* * in the 1024 FFT window.  The net offset is then 576 and you need 448 "new"* * samples for each iteration to keep the 384 samples of interest centered   * *****************************************************************************/  for (j=0; j<syncsize; j++)  {    if (j < (sync_flush) )      savebuf[j] = savebuf[j+flush];    else      savebuf[j] = *buffer++;/**window data with HANN window***********************************************/    if (j<BLKSIZE)    {      wsamp_r[j] = window[j]*((FLOAT) savebuf[j]);       wsamp_i[j] = 0;    }  }/**Compute FFT****************************************************************/        fft(wsamp_r,wsamp_i,energy,phi,1024);/***************************************************************************** * calculate the unpredictability measure, given energy[f] and phi[f]        * *****************************************************************************/        for(j=0; j<HBLKSIZE; j++){           r_prime = 2.0 * r[chn][old][j] - r[chn][oldest][j];           phi_prime = 2.0 * phi_sav[chn][old][j] - phi_sav[chn][oldest][j];           r[chn][new][j] = sqrt((double) energy[j]);           phi_sav[chn][new][j] = phi[j];	   temp1 = r[chn][new][j] * cos((double) phi[j])		   - r_prime * cos(phi_prime);	   temp2=r[chn][new][j] * sin((double) phi[j])		   - r_prime * sin(phi_prime);           temp3=r[chn][new][j] + fabs(r_prime);           if(temp3 != 0)c[j]=sqrt(temp1*temp1+temp2*temp2)/temp3;           else c[j] = 0;        }/*only update data "age" pointers after you are done with the second channel *//*for layer 1 computations, for the layer 2 double computations, the pointers*//*are reset automatically on the second pass                                 */        if(lay==2 || chn==1){           if(new==0){new = 1; oldest = 1;}           else {new = 0; oldest = 0;}           if(old==0)old = 1; else old = 0;        }/***************************************************************************** * Calculate the grouped, energy-weighted, unpredictability measure,         * * grouped_c[], and the grouped energy. grouped_e[]                          * *****************************************************************************/        for(j=1;j<CBANDS;j++){           grouped_e[j] = 0;           grouped_c[j] = 0;        }        grouped_e[0] = energy[0];        grouped_c[0] = energy[0]*c[0];        for(j=1;j<HBLKSIZE;j++){           grouped_e[partition[j]] += energy[j];           grouped_c[partition[j]] += energy[j]*c[j];        }/***************************************************************************** * convolve the grouped energy-weighted unpredictability measure             * * and the grouped energy with the spreading function, s[j][k]               * *****************************************************************************/        for(j=0;j<CBANDS;j++){