/********************************************************************** * ISO MPEG Audio Subgroup Software Simulation Group (1996) * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension * * $Id: psy.c,v 1.1 1996/02/14 04:04:23 rowlands Exp $ * * $Log: psy.c,v $ * 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".   * * 7/24/92  M. Iwadare          HANN window coefficients modified.    * * 7/27/92  Masahiro Iwadare    Bug fix, FFT modification for Layer 3 * * 7/27/92  Masahiro Iwadare    Bug fix, "new", "old", and "oldest"   * *                              updates                               * * 8/07/92  Mike Coleman        Bug fix, read_absthr()                * * 95/3/21  Jon Rowlands        Removed extra debug statements        * **********************************************************************/#include "common.h"#include "encoder.h"FILE		*fpo;	/* file pointer */void psycho_anal(buffer,savebuf,chn,lay,snr32,sfreq)short int *buffer;short int savebuf[1056];int   chn, lay;FLOAT snr32[32];double sfreq;        /* to match prototype : float args are always double */{ unsigned int   i, j, k; FLOAT          r_prime, phi_prime; FLOAT          freq_mult, bval_lo, minthres, sum_energy; double         tb, temp1, temp2, temp3;/* 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 int     new = 0, old = 1, oldest = 0; static int     init = 0, flush, sync_flush, syncsize, sfreq_idx;/* 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; 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"); 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");     fpo = fopen("out.dat", "wb");	if(fpo==NULL) {		puts("\t The attempt to open the output file failed.\n");		exit(-1);}     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");     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);     if(lay==1){        flush = 384;        syncsize = 1024;        sync_flush = 576;     }     else {        flush = 384*3.0/2.0;        syncsize = 1056;        sync_flush = syncsize - flush;     }/* 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));/* 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];