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 Multichannel Encoder
 *
 * $Id: tonal.c 1.8 1996/02/12 07:13:35 rowlands Exp $
 *
 * $Log: tonal.c $
 * Revision 1.8  1996/02/12 07:13:35  rowlands
 * Release following Munich meeting
 *
 * Revision 1.5.2.1  1995/11/06  04:19:12  rowlands
 * Received from Uwe Felderhoff (IRT)
 *
 * Revision 1.7  1995/08/14  08:06:37  tenkate
 * ML-LSF added Warner ten Kate 7/8/95 (Philips)
 * II_psycho_one() split into II_psycho_one() for audio data and
 * II_psycho_one_ml() for MultiLingual data.
 * Variables crit_band, cbound and sub_size has been made local to
 * these functions.
 * Tables "cb" and "th" are copied from LSF-directory.
 * ltg has found its ml counterpart.
 *
 * Revision 1.6  1995/07/31  07:47:50  tenkate
 * bugs correction (if center==1 in Psycho_One), 25/07/95 WtK
 *
 * Revision 1.4.2.1  1995/06/16  03:46:42  rowlands
 * Input from Susanne Ritscher (IRT)
 *
 **********************************************************************/

/**********************************************************************
 *   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        *
 * 92-11-06 Soren H. Nielsen	Corrected power calculation in I_ and *
 *				II_f_f_t.		              *	
 **********************************************************************
 *                                                                    *
 *                                                                    *
 *  MPEG/audio Phase 2 coding/decoding multichannel                   *
 *                                                                    *
 *  7/27/93        Susanne Ritscher,  IRT Munich                      *
 *  8/27/93        Susanne Ritscher, IRT Munich                       *
 *                 Channel-Switching is working                       *
 *  9/1/93         Susanne Ritscher,  IRT Munich                      *
 *                 all channels normalized                            *
 *                                                                    *
 *  9/20/93        channel-switching is only performed at a           *
 *                 certain limit of TC_ALLOC dB, which is included    *
 *                 in encoder.h                                       *
 *                                                                    *
 * 10/18/93        seperated smr and ltmin                            *
 *                                                                    *
 *  Version 1.0                                                       *
 *                                                                    *
 *  07/12/94       Susanne Ritscher,  IRT Munich                      *
 *                 Tel: [phone]                              *
 *                 Fax: [phone]                              *
 *                                                                    *
 *  Version 1.1                                                       *
 *                                                                    *
 *  02/23/95	   Susanne Ritscher,  IRT Munich                      *
 *                 corrected some bugs                                *
 *                 extension bitstream is working                     *
 *                                                                    *
 **********************************************************************/
#define TONAL_WIDTH 0.1 /* When more than one tonal component is within
                           this width in Bark, the weaker one(s) are
                           eliminated */


#include "common.h"
#include "encoder.h"

/**********************************************************************/
/*
/*        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 read_crit_band(int lay, int freq) 
{
 int crit_band;
 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(0);
 }
 fgets(t,80,fp);
 sscanf(t,"%d\n",&crit_band);
 fclose(fp);
 return(crit_band);
}        

void read_cbound(int lay, int freq, int crit_band, int *cbound) 
 /* this function reads in critical 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(0);
 }
 fgets(t,80,fp);               /* skip input for critical bands */
 sscanf(t,"%d\n",&k);
 for(i=0;i<crit_band;i++){ 
    fgets(t,80,fp);
    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(0);
    }
 }
 fclose(fp);
}        

void read_freq_band(int *sub_size, g_ptr *ltg, int lay, int freq) 
  /* this function reads in frequency bands and bark values  */
{
 int i,j, k;
 double a,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(0);
 }
 fgets(t,80,fp);              /* read input for freq. subbands */
 sscanf(t,"%d\n",sub_size);
 *ltg = (g_ptr) 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(0);
    }
 }
 fclose(fp);
}

void make_map (int sub_size, mask *power, g_thres *ltg)
/* this function calculates the global masking threshold */
{
    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(double a, double 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 (double *sample, mask *power) /* this function calculates an */
					    /* FFT analysis for the freq.  */
					    /* domain                      */
{
    int i,j,k,ll,l=0;
    int ip, le, le1;
    double t_r, t_i, u_r, u_i;
    static int M, MM1, init = 0, N, NV2, NM1;
    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");

    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;
	NV2 = FFT_SIZE >> 1;
	NM1 = FFT_SIZE - 1;
	for (ll = 0; ll < M; ll++)
	{
	    le = 1 << (M-ll);
	    le1 = le >> 1;
	    w_r[ll] = cos (PI / le1);
	    w_i[ll] = -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;
    }

    for (i = 0; i < FFT_SIZE; i++)
    {
       x_r[i] = sample[i];
       x_i[i] = energy[i] = 0;
    }
/*
    memcpy ((char *) x_r, (char *) sample, sizeof (double) * FFT_SIZE);
*/
    for (ll = 0; ll < MM1; ll++)
    {
	le = 1 << (M - ll);
	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[ll] - u_i * w_i[ll];
	    u_i = u_i * w_r[ll] + t_r * w_i[ll];
	}
    }
    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 calculation corrected with a factor 4, both positive
	   and negative frequencies exist, 1992-11-06 shn */
	power[i].x = 10 * log10 (energy[i] * 4.0) + 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 (double *sample)   /* this function calculates a  */
				    /* 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 *window;