tonal.c

Motion JPEG编解码器源代码

/**********************************************************************
Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
tonal.c
**********************************************************************/
/**********************************************************************
 * MPEG/audio coding/decoding software, work in progress              *
 *   NOT for public distribution until verified and approved by the   *
 *   MPEG/audio committee.  For further information, please contact   *
 *   Davis Pan, 508-493-2241, e-mail: [e-mail]             *
 *                                                                    *
 * VERSION 3.9t                                                       *
 *   changes made since last update:                                  *
 *   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"     *
 * 2004/7/29 Steven Schultz     cleanup and modernize                 *
 **********************************************************************/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#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.
 *
 **********************************************************************/

extern int crit_band;
extern int *cbound;
extern int sub_size;

void make_map(power, ltg)       /* this function calculates the */
mask power[HAN_SIZE];   /* global masking threshold     */
g_thres *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 sample[FFT_SIZE];  /* FFT analysis for the freq.  */
mask 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 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 *window;

 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)))/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 spike[SBLIMIT];
mask 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 spike[SBLIMIT];
mask 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 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;
 while(first != LAST){               /* the conditions for the tonal          */
    if(first<3 || first>500) run = 0;/* otherwise k+/-j will be out of bounds */
    else if(first<63) run = 2;       /* components in layer II, which         */
    else if(first<127) run = 3;      /* are the boundaries for calc.          */
    else if(first<255) 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<500){     /* 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 */