tonal.c

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

    if (!init)
    {
	/* calculate window function for the Fourier transform */
	window = (double *) 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-1))) / 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.
/*
/*******************************************************************/

void II_pick_max (mask *power, double *spike)
{
 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               */

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

void II_tonal_label(mask *power, int *tone) /* this function extracts (tonal) */
					    /* sinusoidals from the spectrum  */
{
 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;
 while(first != LAST){               /* the conditions for the tonal          */
    if(first<2 || first>499) run = 0;/* otherwise k+/-j will be out of bounds */
    else if(first<62) run = 2;       /* components in layer II, which         */
    else if(first<126) run = 3;      /* are the boundaries for calc.          */
    else if(first<254) run = 6;      /* the tonal components                  */
    else run = 12;
    max = power[first].x - 7;        /* after calculation of tonal   */
    for(j=2;j<=run;j++)              /* components, set to local max */
       if(max < power[first-j].x || max < power[first+j].x){
          power[first].type = FALSE;
          break;
       }
    if(power[first].type == TONE){   /* extract tonal components */
       int help=first;
       if(*tone==LAST) *tone = first;
       while((power[help].next!=LAST)&&(power[help].next-first)<=run)
          help=power[help].next;
       help=power[help].next;
       power[first].next=help;
       if((first-last)<=run){
          if(last_but_one != LAST) power[last_but_one].next=first;
       }
       if(first>1 && first<255){     /* calculate the sum of the */
          double tmp;                /* powers of the components */
          tmp = add_db(power[first-1].x, power[first+1].x);
          power[first].x = add_db(power[first].x, tmp);
       }
       for(j=1;j<=run;j++){
          power[first-j].x = power[first+j].x = DBMIN;
          power[first-j].next = power[first+j].next = STOP;
          power[first-j].type = power[first+j].type = FALSE;
       }
       last_but_one=last;
       last = first;
       first = power[first].next;
    }
    else {
       int ll;
       if(last == LAST); /* *tone = power[first].next; dpwe */
       else power[last].next = power[first].next;
       ll = first;
       first = power[first].next;
       power[ll].next = STOP;
    }
 }
}

/****************************************************************/
/*
/*        This function groups all the remaining non-tonal
/* spectral lines into critical band where they are replaced by
/* one single line.
/*
/****************************************************************/
        
void noise_label(int crit_band, int *cbound, mask *power, int *noise, g_thres *ltg)
{
 int i,j, centre, last = LAST;
 double index, weight, sum;
                              /* calculate the remaining spectral */
 for(i=0;i<crit_band-1;i++){  /* lines for non-tonal components   */
     for(j=cbound[i],weight = 0.0,sum = DBMIN;j<cbound[i+1];j++){
        if(power[j].type != TONE){
           if(power[j].x != DBMIN){
              sum = add_db(power[j].x,sum);
              weight += pow(10.0, power[j].x/10.0) * (ltg[power[j].map].bark-i);
              power[j].x = DBMIN;
           }
        }   /* check to see if the spectral line is low dB, and if  */
     }      /* so replace the center of the critical band, which is */
            /* the center freq. of the noise component              */
     if(sum <= DBMIN)  centre = (cbound[i+1]+cbound[i]) /2;
     else {
        index = weight/pow(10.0,sum/10.0);
        centre = cbound[i] + (int) (index * (double) (cbound[i+1]-cbound[i]) );
     }     /* locate next non-tonal component until finished; */
           /* add to list of non-tonal components             */
     if(power[centre].type == TONE) centre++;
     if(last == LAST) *noise = centre;
     else {
        power[centre].next = LAST;
        power[last].next = centre;
     }
     power[centre].x = sum;
     power[centre].type = NOISE;        
     last = centre;
 }        
}

/****************************************************************/
/*
/*        This function reduces the number of noise and tonal
/* component for further threshold analysis.
/*
/****************************************************************/

void subsampling (mask *power, g_thres *ltg, int *tone, int *noise)
{
    int i, old;
   
    i = *tone;
    old = STOP;
    /* calculate tonal components for */
    while (i != LAST)
    {
	/* reduction of spectral lines    */
	if (power[i].x < ltg[power[i].map].hear)
	{
	    power[i].type = FALSE;
	    power[i].x = DBMIN;
	    if (old == STOP)
		*tone = power[i].next;
	    else
		power[old].next = power[i].next;
	}
	else
	    old = i;
	i = power[i].next;
    }
    i = *noise;
    old = STOP;
    /* calculate non-tonal components for */
    while (i != LAST)
    {
	/* reduction of spectral lines        */
        if (power[i].x < ltg[power[i].map].hear)
	{
	    power[i].type = FALSE;
	    power[i].x = DBMIN;
	    if (old == STOP)
		*noise = power[i].next;
	    else
		power[old].next = power[i].next;
	}
	else
	    old = i;
	i = power[i].next;
    }
    i = *tone;
    old = STOP;
    while (i != LAST)
    {
	/* if more than one */
	if (power[i].next == LAST)
	    break;             /* tonal component  */
	if (ltg[power[power[i].next].map].bark -     /* is less than .5  */
	    ltg[power[i].map].bark < TONAL_WIDTH)
	{
	    /* bark, take the   */
	    if (power[power[i].next].x > power[i].x)
	    {
		/* maximum          */
		if (old == STOP)
		    *tone = power[i].next;
		else
		    power[old].next = power[i].next;
		power[i].type = FALSE;
		power[i].x = DBMIN;
	    }
	    else
	    {
		power[power[i].next].type = FALSE;
		power[power[i].next].x = DBMIN;
		power[i].next = power[power[i].next].next;
	    }
	}
	else
	    old = i;
	i = power[i].next;
    }
}


/* ----------------------------------------------------------------------------
The masking function parameters are here set according to the parameters in
the IRT real time implementation. The constant definitions are for convenience.
1993-07-23 shn
---------------------------------------------------------------------------- */

#define AV_TONAL_K    -9.0 /* Masking index, tonal, constant part  [dB] */
#define AV_NOISE_K    -5.0 /* Masking index, noisy, constant part  [dB] */
#define AV_TONAL_DZ   -0.3 /* Masking index, tonal, CBR dependence [dB/Bark] */
#define AV_NOISE_DZ   -0.3 /* Masking index, noisy, CBR dependence [dB/Bark] */

#define LOW_LIM_1     -1.0 /* 1st lower slope from 0         to LOW_LIM_1 [Bark] */
#define LOW_LIM_2     -3.0 /* 2nd lower slope from LOW_LIM_1 to LOW_LIM_2 [Bark] */

#define LOW_DZ_K_1     6.0 /* 1st lower slope, constant part [dB/Bark] */
#define LOW_DZ_SPL_1   0.4 /* 1st lower slope, SPL dependence [dB/(Bark*dB)] */
#define LOW_DZ_MIN_1  17.0 /* 1st lower slope, minimum value [dB/Bark] */
#define LOW_DZ_2      17.0 /* 2nd lower slope [dB/Bark] */

#define UP_LIM_1       1.0 /* 1st upper slope from 0        to UP_LIM_1 [Bark] */
#define UP_LIM_2       8.0 /* 2nd upper slope from UP_LIM_1 to UP_LIM_2 [Bark] */

#define UP_DZ_1      -18.0 /* 1st upper slope, constant part [dB/Bark] */
#define UP_SPL_1       0.0 /* 1st upper slope, SPL dependence [dB/(Bark*dB)] */
#define UP_DZ_2      -17.0 /* 2nd upper slope, constant part [dB/Bark] */
#define UP_SPL_2      -0.1 /* 2nd upper slope, SPL dependence [dB/(Bark*dB)] */

#define H_THR_OFFSET -12.0 /* Hearing threshold offset [dB] */
#define H_THR_OS_BR    0 /*96*/   /* Minimum datarate for offset, [kbit/s per channel] */

#define MASK_ADD       2.0 /* Addition of maskers	                [dB] */
#define QUIET_ADD      3.0 /* Addition of masker and threshold in quiet [dB] */

/****************************************************************
*
*        This function calculates the individual threshold and
* sum with the quiet threshold to find the global threshold.
*
****************************************************************/

void threshold(int sub_size, mask *power, g_thres *ltg, int *tone, int *noise, int bit_rate)
{
  int k, t;
  double z, dz, spl, vf, tmps;

  for (k=1; k<sub_size; k++) /* Target frequencies */
  {
    ltg[k].x = DBMIN;
    t = *tone;          /* calculate individual masking threshold  */
    while(t != LAST)    /* for tonal components, t,  to find LTG   */
    {
      z   = ltg[power[t].map].bark; /* critical band rate of masker */
      dz  = ltg[k].bark - z ;       /* distance of bark value*/
      spl = power[t].x;             /* sound pressure level of masker */

      if (dz >= LOW_LIM_2 && dz <  UP_LIM_2)
      {
        tmps = spl + AV_TONAL_K + AV_TONAL_DZ * z;

        /* masking function for lower & upper slopes */
        if (LOW_LIM_2<=dz && dz<LOW_LIM_1)
          if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
            vf = LOW_DZ_2 * (dz - LOW_LIM_1) + 
                (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * LOW_LIM_1;
          else
            vf = LOW_DZ_2 * (dz - LOW_LIM_1) + LOW_DZ_MIN_1 * LOW_LIM_1;
        else if (LOW_LIM_1<=dz && dz<0)
          if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
            vf = (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * dz;
          else
            vf = LOW_DZ_MIN_1 * dz;
        else if (0<=dz && dz<UP_LIM_1)
          vf = (UP_DZ_1 * dz);
        else if (UP_LIM_1<=dz && dz<UP_LIM_2)
          vf = (dz - UP_LIM_1) * (UP_DZ_2 - UP_SPL_2 * spl) +
               UP_DZ_1 * UP_LIM_1;
        tmps += vf;        
        ltg[k].x = non_lin_add(ltg[k].x, tmps, MASK_ADD);
      }
      t = power[t].next;
    } /* while */

    t = *noise;        /* calculate individual masking threshold   */
    while(t != LAST)   /* for non-tonal components, t, to find LTG */
    {
      z   = ltg[power[t].map].bark; /* critical band rate of masker */
      dz  = ltg[k].bark - z ;       /* distance of bark value*/
      spl = power[t].x;             /* sound pressure level of masker */

      if (dz >= LOW_LIM_2 && dz <  UP_LIM_2)
      {
        tmps = spl + AV_NOISE_K + AV_NOISE_DZ * z;

        /* masking function for lower & upper slopes */
        if (LOW_LIM_2<=dz && dz<LOW_LIM_1)
          if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
            vf = LOW_DZ_2 * (dz - LOW_LIM_1) + 
                (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * LOW_LIM_1;
          else
            vf = LOW_DZ_2 * (dz - LOW_LIM_1) + LOW_DZ_MIN_1 * LOW_LIM_1;
        else if (LOW_LIM_1<=dz && dz<0)
          if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
            vf = (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * dz;
          else
            vf = LOW_DZ_MIN_1 * dz;
        else if (0<=dz && dz<UP_LIM_1)
          vf = (UP_DZ_1 * dz);
        else if (UP_LIM_1<=dz && dz<UP_LIM_2)
          vf = (dz - UP_LIM_1) * (UP_DZ_2 - UP_SPL_2 * spl) +
               UP_DZ_1 * UP_LIM_1;
       tmps += vf;