tonal.c

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

/**********************************************************************
 * ISO MPEG Audio Subgroup Software Simulation Group (1996)
 * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension
 *
 * $Id: tonal.c,v 1.1 1996/02/14 04:04:23 rowlands Exp $
 *
 * $Log: tonal.c,v $
 * Revision 1.1  1996/02/14 04:04:23  rowlands
 * Initial revision
 *
 * Received from Mike Coleman
 **********************************************************************/
/**********************************************************************
 *   date   programmers         comment                               *
 * 2/25/91  Douglas Wong        start of version 1.1 records          *
 * 3/06/91  Douglas Wong        rename: setup.h to endef.h            *
 *                              updated I_psycho_one and II_psycho_one*
 * 3/11/91  W. J. Carter        Added Douglas Wong's updates dated    *
 *                              3/9/91 for I_Psycho_One() and for     *
 *                              II_Psycho_One().                      *
 * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
 *                              Located and fixed numerous software   *
 *                              bugs and table data errors.           *
 * 6/11/91  Davis Pan           corrected several bugs                *
 *                              based on comments from H. Fuchs       *
 * 01jul91  dpwe (Aware Inc.)   Made pow() args float                 *
 *                              Removed logical bug in I_tonal_label: *
 *                              Sometimes *tone returned == STOP      *
 * 7/10/91  Earle Jennings      no change necessary in port to MsDos  *
 * 11sep91  [e-mail]      Subtracted 90.3dB from II_f_f_t peaks *
 * 10/1/91  Peter W. Farrett    Updated II_Psycho_One(),I_Psycho_One()*
 *                              to include comments.                  *
 *11/29/91  Masahiro Iwadare    Bug fix regarding POWERNORM           *
 *                              fixed several other miscellaneous bugs*
 * 2/11/92  W. Joseph Carter    Ported new code to Macintosh.  Most   *
 *                              important fixes involved changing     *
 *                              16-bit ints to long or unsigned in    *
 *                              bit alloc routines for quant of 65535 *
 *                              and passing proper function args.     *
 *                              Removed "Other Joint Stereo" option   *
 *                              and made bitrate be total channel     *
 *                              bitrate, irrespective of the mode.    *
 *                              Fixed many small bugs & reorganized.  *
 * 2/12/92  Masahiro Iwadare    Fixed some potential bugs in          *
 *          Davis Pan           subsampling()                         *
 * 2/25/92  Masahiro Iwadare    Fixed some more potential bugs        *
 * 6/24/92  Tan Ah Peng         Modified window for FFT               * 
 *                              (denominator N-1 to N)                *
 *                              Updated all critical band rate &      *
 *                              absolute threshold tables and critical*
 *                              boundaries for use with Layer I & II  *  
 *                              Corrected boundary limits for tonal   *
 *                              component computation                 *
 *                              Placement of non-tonal component at   *
 *                              geometric mean of critical band       *
 *                              (previous placement method commented  *
 *                               out - can be used if desired)        *
 * 3/01/93  Mike Li             Infinite looping fix in noise_label() *
 * 3/19/93  Jens Spille         fixed integer overflow problem in     *
 *                              psychoacoutic model 1                 *
 * 3/19/93  Giorgio Dimino      modifications to better account for   *
 *                              tonal and non-tonal components        *
 * 5/28/93 Sriram Jayasimha     "London" mod. to psychoacoustic model1*
 * 8/05/93 Masahiro Iwadare     noise_label modification "option"     *
 * 1/21/94 Seymore Shlien       fixed another infinite looping problem*
 * 7/12/95 Soeren H. Nielsen    Changes for LSF, new tables           *
 **********************************************************************/

#include "common.h"
#include "encoder.h"
#define LONDON                  /* enable "LONDON" modification */
#define MAKE_SENSE              /* enable "MAKE_SENSE" modification */
#define MI_OPTION               /* enable "MI_OPTION" modification */
/**********************************************************************
*
*        This module implements the psychoacoustic model I for the
* MPEG encoder layer II. It uses simplified tonal and noise masking
* threshold analysis to generate SMR for the encoder bit allocation
* routine.
*
**********************************************************************/

int crit_band;
int FAR *cbound;
int sub_size;

void read_cbound(lay,freq)  /* this function reads in critical */
int lay, freq;              /* band boundaries                 */
{
 int i,j,k;
 FILE *fp;
 char r[16], t[80];

 strcpy(r, "2cb1");
 r[0] = (char) lay + '0';
 r[3] = (char) freq + '0';
 if( !(fp = OpenTableFile(r)) ){       /* check boundary values */
    printf("Please check %s boundary table\n",r);
    exit(1);
 }
 fgets(t,80,fp);               /* read input for critical bands */
 sscanf(t,"%d\n",&crit_band);
 cbound = (int FAR *) mem_alloc(sizeof(int) * crit_band, "cbound");
 for(i=0;i<crit_band;i++){   /* continue to read input for */
    fgets(t,80,fp);            /* critical band boundaries   */
    sscanf(t,"%d %d\n",&j, &k);
    if(i==j) cbound[j] = k;
    else {                     /* error */
       printf("Please check index %d in cbound table %s\n",i,r);
       exit(1);
    }
 }
 fclose(fp);
}        

void read_freq_band(ltg,lay,freq)  /* this function reads in   */
int lay, freq;                     /* frequency bands and bark */
g_ptr FAR *ltg;                /* values                   */
{
 int i,j, k;
 double b,c;
 FILE *fp;
 char r[16], t[80];

 strcpy(r, "2th1");
 r[0] = (char) lay + '0';
 r[3] = (char) freq + '0';
 if( !(fp = OpenTableFile(r)) ){   /* check freq. values  */
    printf("Please check frequency and cband table %s\n",r);
    exit(1);
 }
 fgets(t,80,fp);              /* read input for freq. subbands */
 sscanf(t,"%d\n",&sub_size);
 *ltg = (g_ptr FAR ) mem_alloc(sizeof(g_thres) * sub_size, "ltg");
 (*ltg)[0].line = 0;          /* initialize global masking threshold */
 (*ltg)[0].bark = 0;
 (*ltg)[0].hear = 0;
 for(i=1;i<sub_size;i++){    /* continue to read freq. subband */
    fgets(t,80,fp);          /* and assign                     */
    sscanf(t,"%d %d %lf %lf\n",&j, &k, &b, &c);
    if(i == j){
       (*ltg)[j].line = k;
       (*ltg)[j].bark = b;
       (*ltg)[j].hear = c;
    }
    else {                   /* error */
       printf("Please check index %d in freq-cb table %s\n",i,r);
       exit(1);
    }
 }
 fclose(fp);
}

void make_map(power, ltg)       /* this function calculates the */
mask FAR power[HAN_SIZE];   /* global masking threshold     */
g_thres FAR *ltg;
{
 int i,j;

 for(i=1;i<sub_size;i++) for(j=ltg[i-1].line;j<=ltg[i].line;j++)
    power[j].map = i;
}

double add_db(a,b)
double a,b;
{
 a = pow(10.0,a/10.0);
 b = pow(10.0,b/10.0);
 return 10 * log10(a+b);
}

/****************************************************************
*
*        Fast Fourier transform of the input samples.
*
****************************************************************/

void II_f_f_t(sample, power)      /* this function calculates an */
double FAR sample[FFT_SIZE];  /* FFT analysis for the freq.  */
mask FAR power[HAN_SIZE];     /* domain                      */
{
 int i,j,k,L,l=0;
 int ip, le, le1;
 double t_r, t_i, u_r, u_i;
 static int M, MM1, init = 0, N;
 double *x_r, *x_i, *energy;
 static int *rev;
 static double *w_r, *w_i;

 x_r = (double *) mem_alloc(sizeof(DFFT), "x_r");
 x_i = (double *) mem_alloc(sizeof(DFFT), "x_i");
 energy = (double *) mem_alloc(sizeof(DFFT), "energy");
 for(i=0;i<FFT_SIZE;i++) x_r[i] = x_i[i] = energy[i] = 0;
 if(!init){
    rev = (int *) mem_alloc(sizeof(IFFT), "rev");
    w_r = (double *) mem_alloc(sizeof(D10), "w_r");
    w_i = (double *) mem_alloc(sizeof(D10), "w_i");
    M = 10;
    MM1 = 9;
    N = FFT_SIZE;
    for(L=0;L<M;L++){
       le = 1 << (M-L);
       le1 = le >> 1;
       w_r[L] = cos(PI/le1);
       w_i[L] = -sin(PI/le1);
    }
    for(i=0;i<FFT_SIZE;rev[i] = l,i++) for(j=0,l=0;j<10;j++){
       k=(i>>j) & 1;
       l |= (k<<(9-j));                
    }
    init = 1;
 }
 memcpy( (char *) x_r, (char *) sample, sizeof(double) * FFT_SIZE);
 for(L=0;L<MM1;L++){
    le = 1 << (M-L);
    le1 = le >> 1;
    u_r = 1;
    u_i = 0;
    for(j=0;j<le1;j++){
       for(i=j;i<N;i+=le){
          ip = i + le1;
          t_r = x_r[i] + x_r[ip];
          t_i = x_i[i] + x_i[ip];
          x_r[ip] = x_r[i] - x_r[ip];
          x_i[ip] = x_i[i] - x_i[ip];
          x_r[i] = t_r;
          x_i[i] = t_i;
          t_r = x_r[ip];
          x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
          x_i[ip] = x_i[ip] * u_r + t_r * u_i;
       }
       t_r = u_r;
       u_r = u_r * w_r[L] - u_i * w_i[L];
       u_i = u_i * w_r[L] + t_r * w_i[L];
    }
 }
 for(i=0;i<N;i+=2){
    ip = i + 1;
    t_r = x_r[i] + x_r[ip];
    t_i = x_i[i] + x_i[ip];
    x_r[ip] = x_r[i] - x_r[ip];
    x_i[ip] = x_i[i] - x_i[ip];
    x_r[i] = t_r;
    x_i[i] = t_i;
    energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
 }
 for(i=0;i<FFT_SIZE;i++) if(i<rev[i]){
    t_r = energy[i];
    energy[i] = energy[rev[i]];
    energy[rev[i]] = t_r;
 }
 for(i=0;i<HAN_SIZE;i++){    /* calculate power density spectrum */
    if (energy[i] < 1E-20) energy[i] = 1E-20;
    power[i].x = 10 * log10(energy[i]) + POWERNORM;
    power[i].next = STOP;
    power[i].type = FALSE;
 }
 mem_free((void **) &x_r);
 mem_free((void **) &x_i);
 mem_free((void **) &energy);
}

/****************************************************************
*
*         Window the incoming audio signal.
*
****************************************************************/

void II_hann_win(sample)          /* this function calculates a  */
double FAR sample[FFT_SIZE];  /* Hann window for PCM (input) */
{                                 /* samples for a 1024-pt. FFT  */
 register int i;
 register double sqrt_8_over_3;
 static int init = 0;
 static double FAR *window;

 if(!init){  /* calculate window function for the Fourier transform */
    window = (double FAR *) mem_alloc(sizeof(DFFT), "window");
    sqrt_8_over_3 = pow(8.0/3.0, 0.5);
    for(i=0;i<FFT_SIZE;i++){
       /* Hann window formula */
       window[i]=sqrt_8_over_3*0.5*(1-cos(2.0*PI*i/(FFT_SIZE)))/FFT_SIZE;
    }
    init = 1;
 }
 for(i=0;i<FFT_SIZE;i++) sample[i] *= window[i];
}

/*******************************************************************
*
*        This function finds the maximum spectral component in each
* subband and return them to the encoder for time-domain threshold
* determination.
*
*******************************************************************/
#ifndef LONDON
void II_pick_max(power, spike)
double FAR spike[SBLIMIT];
mask FAR power[HAN_SIZE];
{
 double max;
 int i,j;

 for(i=0;i<HAN_SIZE;spike[i>>4] = max, i+=16)      /* calculate the      */
 for(j=0, max = DBMIN;j<16;j++)                    /* maximum spectral   */
    max = (max>power[i+j].x) ? max : power[i+j].x; /* component in each  */
}                                                  /* subband from bound */
                                                   /* 4-16               */
#else
void II_pick_max(power, spike)
double FAR spike[SBLIMIT];
mask FAR power[HAN_SIZE];
{
 double sum;
 int i,j;

 for(i=0;i<HAN_SIZE;spike[i>>4] = 10.0*log10(sum), i+=16)
                                                   /* calculate the      */
 for(j=0, sum = pow(10.0,0.1*DBMIN);j<16;j++)      /* sum of spectral   */
   sum += pow(10.0,0.1*power[i+j].x);              /* component in each  */
}                                                  /* subband from bound */
                                                   /* 4-16               */
#endif

/****************************************************************
*
*        This function labels the tonal component in the power
* spectrum.
*
****************************************************************/

void II_tonal_label(power, tone)  /* this function extracts (tonal) */
mask FAR power[HAN_SIZE];     /* sinusoidals from the spectrum  */
int *tone;
{
 int i,j, last = LAST, first, run, last_but_one = LAST; /* dpwe */
 double max;

 *tone = LAST;
 for(i=2;i<HAN_SIZE-12;i++){
    if(power[i].x>power[i-1].x && power[i].x>=power[i+1].x){
       power[i].type = TONE;
       power[i].next = LAST;
       if(last != LAST) power[last].next = i;
       else first = *tone = i;
       last = i;
    }
 }
 last = LAST;
 first = *tone;
 *tone = LAST;