psy.c

ISO mp3 sources (distribution 10) Layer 1/2/3, C Source, 512 k Sources of the Mpeg 1,2 layer 1,2

        }
     }
     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++;
           if(j<BLKSIZE){
/**window data with HANN window***********************************************/
              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]        *
 *****************************************************************************/
/*only update data "age" pointers after you are done with both channels      */
/*for layer 1 computations, for the layer 2 double computations, the pointers*/
/*are reset automatically on the second pass                                 */
         if(lay==2 || (lay==1 && chn==0) ){
           if(new==0){new = 1; oldest = 1;}
           else {new = 0; oldest = 0;}
           if(old==0)old = 1; else old = 0;
        }
        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((double) phi_prime);
temp2=r[chn][new][j] * sin((double) phi[j]) - r_prime * sin((double) phi_prime);
           temp3=r[chn][new][j] + fabs((double)r_prime);
           if(temp3 != 0)c[j]=sqrt(temp1*temp1+temp2*temp2)/temp3;
           else c[j] = 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++){
           ecb[j] = 0;
           cb[j] = 0;
           for(k=0;k<CBANDS;k++){
              if(s[j][k] != 0.0){
                 ecb[j] += s[j][k]*grouped_e[k];
                 cb[j] += s[j][k]*grouped_c[k];
              }
           }
           if(ecb[j] !=0)cb[j] = cb[j]/ecb[j];
           else cb[j] = 0;
        }

/*****************************************************************************
 * Calculate the required SNR for each of the frequency partitions           *
 *         this whole section can be accomplished by a table lookup          *
 *****************************************************************************/
        for(j=0;j<CBANDS;j++){
           if(cb[j]<.05)cb[j]=0.05;
           else if(cb[j]>.5)cb[j]=0.5;
           tb = -0.434294482*log((double) cb[j])-0.301029996;
	   cb[j]=tb;
           bc[j] = tmn[j]*tb + nmt*(1.0-tb);
           k = cbval[j] + 0.5;
           bc[j] = (bc[j] > bmax[k]) ? bc[j] : bmax[k];
           bc[j] = exp((double) -bc[j]*LN_TO_LOG10);
        }

/*****************************************************************************
 * Calculate the permissible noise energy level in each of the frequency     *
 * partitions. Include absolute threshold and pre-echo controls              *
 *         this whole section can be accomplished by a table lookup          *
 *****************************************************************************/
        for(j=0;j<CBANDS;j++)
           if(rnorm[j] && numlines[j])
              nb[j] = ecb[j]*bc[j]/(rnorm[j]*numlines[j]);
           else nb[j] = 0;
        for(j=0;j<HBLKSIZE;j++){
/*temp1 is the preliminary threshold */
           temp1=nb[partition[j]];
           temp1=(temp1>absthr[j])?temp1:absthr[j];
/*do not use pre-echo control for layer 2 because it may do bad things to the*/
/*  MUSICAM bit allocation algorithm                                         */
           if(lay==1){
              fthr[j] = (temp1 < lthr[chn][j]) ? temp1 : lthr[chn][j];
              temp2 = temp1 * 0.00316;
              fthr[j] = (temp2 > fthr[j]) ? temp2 : fthr[j];
           }
           else fthr[j] = temp1;
           lthr[chn][j] = LXMIN*temp1;
        }

/*****************************************************************************
 * Translate the 512 threshold values to the 32 filter bands of the coder    *
 *****************************************************************************/
        for(j=0;j<193;j += 16){
           minthres = 60802371420160.0;
           sum_energy = 0.0;
           for(k=0;k<17;k++){
              if(minthres>fthr[j+k])minthres = fthr[j+k];
              sum_energy += energy[j+k];
           }
           snrtmp[i][j/16] = sum_energy/(minthres * 17.0);
           snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
        }
        for(j=208;j<(HBLKSIZE-1);j += 16){
           minthres = 0.0;
           sum_energy = 0.0;
           for(k=0;k<17;k++){
              minthres += fthr[j+k];
              sum_energy += energy[j+k];
           }
           snrtmp[i][j/16] = sum_energy/minthres;
           snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
        }
/*****************************************************************************
 * End of Psychoacuostic calculation loop                                    *
 *****************************************************************************/
     }
     for(i=0; i<32; i++){
        if(lay==2)
           snr32[i]=(snrtmp[0][i]>snrtmp[1][i])?snrtmp[0][i]:snrtmp[1][i];
        else snr32[i]=snrtmp[0][i];
     }
     break;
  case 3:
     printf("layer 3 is not currently supported\n");
     break;
  default:
     printf("error, invalid MPEG/audio coding layer: %d\n",lay);
 }

/* These mem_free() calls must correspond with the mem_alloc() calls     */
/* used at the beginning of this function to simulate "automatic"        */
/* variables placed on the stack.                                        */

 mem_free((void **) &grouped_c);
 mem_free((void **) &grouped_e);
 mem_free((void **) &nb);
 mem_free((void **) &cb);
 mem_free((void **) &ecb);
 mem_free((void **) &bc);
 mem_free((void **) &wsamp_r);
 mem_free((void **) &wsamp_i);
 mem_free((void **) &phi);
 mem_free((void **) &energy);
 mem_free((void **) &c);
 mem_free((void **) &fthr);
 mem_free((void **) &snrtmp);
}

/******************************************************************************
routine to read in absthr table from a file.
******************************************************************************/

void read_absthr(absthr, table)
FLOAT *absthr;
int table;
{
 FILE *fp;
 long j,index;
 float a;
 char t[80];
 char ta[16];

 strcpy( ta, "absthr_0" );
 
 switch(table){
    case 0 : ta[7] = '0';
             break;
    case 1 : ta[7] = '1';
             break;
    case 2 : ta[7] = '2';
             break;
    default : printf("absthr table: Not valid table number\n");
 }
 if(!(fp = OpenTableFile(ta) ) ){
    printf("Please check %s table\n", ta);
    exit(1);
 }
 fgets(t, 150, fp);
 sscanf(t, "table %ld", &index);
 if(index != table){
    printf("error in absthr table %s",ta);
    exit(1);
 }
 for(j=0; j<HBLKSIZE; j++){
    fgets(t,80,fp);
    sscanf(t,"%f", &a);
    absthr[j] =  a;
 }
 fclose(fp);
}