home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga 5 / MA_Cover_5.iso / ppc / wos_mpeg_encode / mpegaudio / tonal.c < prev   
Encoding:
C/C++ Source or Header  |  1998-05-04  |  35.8 KB  |  988 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. tonal.c
  4. **********************************************************************/
  5. /**********************************************************************
  6.  * MPEG/audio coding/decoding software, work in progress              *
  7.  *   NOT for public distribution until verified and approved by the   *
  8.  *   MPEG/audio committee.  For further information, please contact   *
  9.  *   Davis Pan, 508-493-2241, e-mail: pan@3d.enet.dec.com             *
  10.  *                                                                    *
  11.  * VERSION 3.9t                                                       *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Douglas Wong        start of version 1.1 records          *
  15.  * 3/06/91  Douglas Wong        rename: setup.h to endef.h            *
  16.  *                              updated I_psycho_one and II_psycho_one*
  17.  * 3/11/91  W. J. Carter        Added Douglas Wong's updates dated    *
  18.  *                              3/9/91 for I_Psycho_One() and for     *
  19.  *                              II_Psycho_One().                      *
  20.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  21.  *                              Located and fixed numerous software   *
  22.  *                              bugs and table data errors.           *
  23.  * 6/11/91  Davis Pan           corrected several bugs                *
  24.  *                              based on comments from H. Fuchs       *
  25.  * 01jul91  dpwe (Aware Inc.)   Made pow() args float                 *
  26.  *                              Removed logical bug in I_tonal_label: *
  27.  *                              Sometimes *tone returned == STOP      *
  28.  * 7/10/91  Earle Jennings      no change necessary in port to MsDos  *
  29.  * 11sep91  dpwe@aware.com      Subtracted 90.3dB from II_f_f_t peaks *
  30.  * 10/1/91  Peter W. Farrett    Updated II_Psycho_One(),I_Psycho_One()*
  31.  *                              to include comments.                  *
  32.  *11/29/91  Masahiro Iwadare    Bug fix regarding POWERNORM           *
  33.  *                              fixed several other miscellaneous bugs*
  34.  * 2/11/92  W. Joseph Carter    Ported new code to Macintosh.  Most   *
  35.  *                              important fixes involved changing     *
  36.  *                              16-bit ints to long or unsigned in    *
  37.  *                              bit alloc routines for quant of 65535 *
  38.  *                              and passing proper function args.     *
  39.  *                              Removed "Other Joint Stereo" option   *
  40.  *                              and made bitrate be total channel     *
  41.  *                              bitrate, irrespective of the mode.    *
  42.  *                              Fixed many small bugs & reorganized.  *
  43.  * 2/12/92  Masahiro Iwadare    Fixed some potential bugs in          *
  44.  *          Davis Pan           subsampling()                         *
  45.  * 2/25/92  Masahiro Iwadare    Fixed some more potential bugs        *
  46.  * 6/24/92  Tan Ah Peng         Modified window for FFT               * 
  47.  *                              (denominator N-1 to N)                *
  48.  *                              Updated all critical band rate &      *
  49.  *                              absolute threshold tables and critical*
  50.  *                              boundaries for use with Layer I & II  *  
  51.  *                              Corrected boundary limits for tonal   *
  52.  *                              component computation                 *
  53.  *                              Placement of non-tonal component at   *
  54.  *                              geometric mean of critical band       *
  55.  *                              (previous placement method commented  *
  56.  *                               out - can be used if desired)        *
  57.  * 3/01/93  Mike Li             Infinite looping fix in noise_label() *
  58.  * 3/19/93  Jens Spille         fixed integer overflow problem in     *
  59.  *                              psychoacoutic model 1                 *
  60.  * 3/19/93  Giorgio Dimino      modifications to better account for   *
  61.  *                              tonal and non-tonal components        *
  62.  * 5/28/93 Sriram Jayasimha     "London" mod. to psychoacoustic model1*
  63.  * 8/05/93 Masahiro Iwadare     noise_label modification "option"     *
  64.  **********************************************************************/
  65.  
  66. #include "common.h"
  67. #include "encoder.h"
  68. #define LONDON                  /* enable "LONDON" modification */
  69. #define MAKE_SENSE              /* enable "MAKE_SENSE" modification */
  70. #define MI_OPTION               /* enable "MI_OPTION" modification */
  71. /**********************************************************************/
  72. /*
  73. /*        This module implements the psychoacoustic model I for the
  74. /* MPEG encoder layer II. It uses simplified tonal and noise masking
  75. /* threshold analysis to generate SMR for the encoder bit allocation
  76. /* routine.
  77. /*
  78. /**********************************************************************/
  79.  
  80. int crit_band;
  81. int FAR *cbound;
  82. int sub_size;
  83.  
  84. void read_cbound(int lay,int freq)  /* this function reads in critical */
  85. {
  86.  int i,j,k;
  87.  FILE *fp;
  88.  char r[16], t[80];
  89.  
  90.  strcpy(r, "2cb1");
  91.  r[0] = (char) lay + '0';
  92.  r[3] = (char) freq + '0';
  93.  if( !(fp = OpenTableFile(r)) ){       /* check boundary values */
  94.     printf("Please check %s boundary table\n",r);
  95.     exit(1);
  96.  }
  97.  fgets(t,80,fp);               /* read input for critical bands */
  98.  sscanf(t,"%d\n",&crit_band);
  99.  cbound = (int FAR *) mem_alloc(sizeof(int) * crit_band, "cbound");
  100.  for(i=0;i<crit_band;i++){   /* continue to read input for */
  101.     fgets(t,80,fp);            /* critical band boundaries   */
  102.     sscanf(t,"%d %d\n",&j, &k);
  103.     if(i==j) cbound[j] = k;
  104.     else {                     /* error */
  105.        printf("Please check index %d in cbound table %s\n",i,r);
  106.        exit(1);
  107.     }
  108.  }
  109.  fclose(fp);
  110. }        
  111.  
  112. void read_freq_band(g_ptr FAR *ltg,int lay,int freq)  /* this function reads in   */
  113. {
  114.  int i,j, k;
  115.  double b,c;
  116.  FILE *fp;
  117.  char r[16], t[80];
  118.  
  119.  strcpy(r, "2th1");
  120.  r[0] = (char) lay + '0';
  121.  r[3] = (char) freq + '0';
  122.  if( !(fp = OpenTableFile(r)) ){   /* check freq. values  */
  123.     printf("Please check frequency and cband table %s\n",r);
  124.     exit(1);
  125.  }
  126.  fgets(t,80,fp);              /* read input for freq. subbands */
  127.  sscanf(t,"%d\n",&sub_size);
  128.  *ltg = (g_ptr FAR ) mem_alloc(sizeof(g_thres) * sub_size, "ltg");
  129.  (*ltg)[0].line = 0;          /* initialize global masking threshold */
  130.  (*ltg)[0].bark = 0;
  131.  (*ltg)[0].hear = 0;
  132.  for(i=1;i<sub_size;i++){    /* continue to read freq. subband */
  133.     fgets(t,80,fp);          /* and assign                     */
  134.     sscanf(t,"%d %d %lf %lf\n",&j, &k, &b, &c);
  135.     if(i == j){
  136.        (*ltg)[j].line = k;
  137.        (*ltg)[j].bark = b;
  138.        (*ltg)[j].hear = c;
  139.     }
  140.     else {                   /* error */
  141.        printf("Please check index %d in freq-cb table %s\n",i,r);
  142.        exit(1);
  143.     }
  144.  }
  145.  fclose(fp);
  146. }
  147.  
  148. void make_map(mask FAR power[HAN_SIZE], g_thres FAR *ltg)       /* this function calculates the */
  149. {
  150.  int i,j;
  151.  
  152.  for(i=1;i<sub_size;i++) for(j=ltg[i-1].line;j<=ltg[i].line;j++)
  153.     power[j].map = i;
  154. }
  155.  
  156. double add_db(double a,double b)
  157. {
  158.  a = pow(10.0,a/10.0);
  159.  b = pow(10.0,b/10.0);
  160.  return 10 * log10(a+b);
  161. }
  162.  
  163. /****************************************************************/
  164. /*
  165. /*        Fast Fourier transform of the input samples.
  166. /*
  167. /****************************************************************/
  168.  
  169. void II_f_f_t(double FAR sample[FFT_SIZE], mask FAR power[HAN_SIZE])      /* this function calculates an */
  170. {
  171.  int i,j,k,L,l=0;
  172.  int ip, le, le1;
  173.  double t_r, t_i, u_r, u_i;
  174.  static int M, MM1, init = 0, N;
  175.  double *x_r, *x_i, *energy;
  176.  static int *rev;
  177.  static double *w_r, *w_i;
  178.  
  179.  x_r = (double *) mem_alloc(sizeof(DFFT), "x_r");
  180.  x_i = (double *) mem_alloc(sizeof(DFFT), "x_i");
  181.  energy = (double *) mem_alloc(sizeof(DFFT), "energy");
  182.  for(i=0;i<FFT_SIZE;i++) x_r[i] = x_i[i] = energy[i] = 0;
  183.  if(!init){
  184.     rev = (int *) mem_alloc(sizeof(IFFT), "rev");
  185.     w_r = (double *) mem_alloc(sizeof(D10), "w_r");
  186.     w_i = (double *) mem_alloc(sizeof(D10), "w_i");
  187.     M = 10;
  188.     MM1 = 9;
  189.     N = FFT_SIZE;
  190.     for(L=0;L<M;L++){
  191.        le = 1 << (M-L);
  192.        le1 = le >> 1;
  193.        w_r[L] = cos(PI/le1);
  194.        w_i[L] = -sin(PI/le1);
  195.     }
  196.     for(i=0;i<FFT_SIZE;rev[i] = l,i++) for(j=0,l=0;j<10;j++){
  197.        k=(i>>j) & 1;
  198.        l |= (k<<9-j);                
  199.     }
  200.     init = 1;
  201.  }
  202.  memcpy( (char *) x_r, (char *) sample, sizeof(double) * FFT_SIZE);
  203.  for(L=0;L<MM1;L++){
  204.     le = 1 << (M-L);
  205.     le1 = le >> 1;
  206.     u_r = 1;
  207.     u_i = 0;
  208.     for(j=0;j<le1;j++){
  209.        for(i=j;i<N;i+=le){
  210.           ip = i + le1;
  211.           t_r = x_r[i] + x_r[ip];
  212.           t_i = x_i[i] + x_i[ip];
  213.           x_r[ip] = x_r[i] - x_r[ip];
  214.           x_i[ip] = x_i[i] - x_i[ip];
  215.           x_r[i] = t_r;
  216.           x_i[i] = t_i;
  217.           t_r = x_r[ip];
  218.           x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
  219.           x_i[ip] = x_i[ip] * u_r + t_r * u_i;
  220.        }
  221.        t_r = u_r;
  222.        u_r = u_r * w_r[L] - u_i * w_i[L];
  223.        u_i = u_i * w_r[L] + t_r * w_i[L];
  224.     }
  225.  }
  226.  for(i=0;i<N;i+=2){
  227.     ip = i + 1;
  228.     t_r = x_r[i] + x_r[ip];
  229.     t_i = x_i[i] + x_i[ip];
  230.     x_r[ip] = x_r[i] - x_r[ip];
  231.     x_i[ip] = x_i[i] - x_i[ip];
  232.     x_r[i] = t_r;
  233.     x_i[i] = t_i;
  234.     energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
  235.  }
  236.  for(i=0;i<FFT_SIZE;i++) if(i<rev[i]){
  237.     t_r = energy[i];
  238.     energy[i] = energy[rev[i]];
  239.     energy[rev[i]] = t_r;
  240.  }
  241.  for(i=0;i<HAN_SIZE;i++){    /* calculate power density spectrum */
  242.     if (energy[i] < 1E-20) energy[i] = 1E-20;
  243.     power[i].x = 10 * log10(energy[i]) + POWERNORM;
  244.     power[i].next = STOP;
  245.     power[i].type = FALSE;
  246.  }
  247.  mem_free((void **) &x_r);
  248.  mem_free((void **) &x_i);
  249.  mem_free((void **) &energy);
  250. }
  251.  
  252. /****************************************************************/
  253. /*
  254. /*         Window the incoming audio signal.
  255. /*
  256. /****************************************************************/
  257.  
  258. void II_hann_win(double FAR sample[FFT_SIZE])          /* this function calculates a  */
  259. {                                 /* samples for a 1024-pt. FFT  */
  260.  register int i;
  261.  register double sqrt_8_over_3;
  262.  static int init = 0;
  263.  static double FAR *window;
  264.  
  265.  if(!init){  /* calculate window function for the Fourier transform */
  266.     window = (double FAR *) mem_alloc(sizeof(DFFT), "window");
  267.     sqrt_8_over_3 = pow(8.0/3.0, 0.5);
  268.     for(i=0;i<FFT_SIZE;i++){
  269.        /* Hann window formula */
  270.        window[i]=sqrt_8_over_3*0.5*(1-cos(2.0*PI*i/(FFT_SIZE)))/FFT_SIZE;
  271.     }
  272.     init = 1;
  273.  }
  274.  for(i=0;i<FFT_SIZE;i++) sample[i] *= window[i];
  275. }
  276.  
  277. /*******************************************************************/
  278. /*
  279. /*        This function finds the maximum spectral component in each
  280. /* subband and return them to the encoder for time-domain threshold
  281. /* determination.
  282. /*
  283. /*******************************************************************/
  284. #ifndef LONDON
  285. void II_pick_max(mask FAR power[HAN_SIZE], double FAR spike[SBLIMIT])
  286. {
  287.  double max;
  288.  int i,j;
  289.  
  290.  for(i=0;i<HAN_SIZE;spike[i>>4] = max, i+=16)      /* calculate the      */
  291.  for(j=0, max = DBMIN;j<16;j++)                    /* maximum spectral   */
  292.     max = (max>power[i+j].x) ? max : power[i+j].x; /* component in each  */
  293. }                                                  /* subband from bound */
  294.                                                    /* 4-16               */
  295. #else
  296. void II_pick_max(mask FAR power[HAN_SIZE], double FAR spike[SBLIMIT])
  297. {
  298.  double sum;
  299.  int i,j;
  300.  
  301.  for(i=0;i<HAN_SIZE;spike[i>>4] = 10.0*log10(sum), i+=16)
  302.                                                    /* calculate the      */
  303.  for(j=0, sum = pow(10.0,0.1*DBMIN);j<16;j++)      /* sum of spectral   */
  304.    sum += pow(10.0,0.1*power[i+j].x);              /* component in each  */
  305. }                                                  /* subband from bound */
  306.                                                    /* 4-16               */
  307. #endif
  308.  
  309. /****************************************************************/
  310. /*
  311. /*        This function labels the tonal component in the power
  312. /* spectrum.
  313. /*
  314. /****************************************************************/
  315.  
  316. void II_tonal_label(mask FAR power[HAN_SIZE], int *tone)  /* this function extracts (tonal) */
  317. {
  318.  int i,j, last = LAST, first, run, last_but_one = LAST; /* dpwe */
  319.  double max;
  320.  
  321.  *tone = LAST;
  322.  for(i=2;i<HAN_SIZE-12;i++){
  323.     if(power[i].x>power[i-1].x && power[i].x>=power[i+1].x){
  324.        power[i].type = TONE;
  325.        power[i].next = LAST;
  326.        if(last != LAST) power[last].next = i;
  327.        else first = *tone = i;
  328.        last = i;
  329.     }
  330.  }
  331.  last = LAST;
  332.  first = *tone;
  333.  *tone = LAST;
  334.  while(first != LAST){               /* the conditions for the tonal          */
  335.     if(first<3 || first>500) run = 0;/* otherwise k+/-j will be out of bounds */
  336.     else if(first<63) run = 2;       /* components in layer II, which         */
  337.     else if(first<127) run = 3;      /* are the boundaries for calc.          */
  338.     else if(first<255) run = 6;      /* the tonal components                  */
  339.     else run = 12;
  340.     max = power[first].x - 7;        /* after calculation of tonal   */
  341.     for(j=2;j<=run;j++)              /* components, set to local max */
  342.        if(max < power[first-j].x || max < power[first+j].x){
  343.           power[first].type = FALSE;
  344.           break;
  345.        }
  346.     if(power[first].type == TONE){   /* extract tonal components */
  347.        int help=first;
  348.        if(*tone==LAST) *tone = first;
  349.        while((power[help].next!=LAST)&&(power[help].next-first)<=run)
  350.           help=power[help].next;
  351.        help=power[help].next;
  352.        power[first].next=help;
  353.        if((first-last)<=run){
  354.           if(last_but_one != LAST) power[last_but_one].next=first;
  355.        }
  356.        if(first>1 && first<500){     /* calculate the sum of the */
  357.           double tmp;                /* powers of the components */
  358.           tmp = add_db(power[first-1].x, power[first+1].x);
  359.           power[first].x = add_db(power[first].x, tmp);
  360.        }
  361.        for(j=1;j<=run;j++){
  362.           power[first-j].x = power[first+j].x = DBMIN;
  363.           power[first-j].next = power[first+j].next = STOP;
  364.           power[first-j].type = power[first+j].type = FALSE;
  365.        }
  366.        last_but_one=last;
  367.        last = first;
  368.        first = power[first].next;
  369.     }
  370.     else {
  371.        int ll;
  372.        if(last == LAST); /* *tone = power[first].next; dpwe */
  373.        else power[last].next = power[first].next;
  374.        ll = first;
  375.        first = power[first].next;
  376.        power[ll].next = STOP;
  377.     }
  378.  }
  379. }
  380.  
  381. /****************************************************************/
  382. /*
  383. /*        This function groups all the remaining non-tonal
  384. /* spectral lines into critical band where they are replaced by
  385. /* one single line.
  386. /*
  387. /****************************************************************/
  388.         
  389. void noise_label(mask FAR *power, int *noise, g_thres FAR *ltg)
  390. {
  391.  int i,j, centre, last = LAST;
  392.  double index, weight, sum;
  393.                               /* calculate the remaining spectral */
  394.  for(i=0;i<crit_band-1;i++){  /* lines for non-tonal components   */
  395.      for(j=cbound[i],weight = 0.0,sum = DBMIN;j<cbound[i+1];j++){
  396.         if(power[j].type != TONE){
  397.            if(power[j].x != DBMIN){
  398.               sum = add_db(power[j].x,sum);
  399. /* the line below and others under the "MAKE_SENSE" condition are an alternate
  400.    interpretation of "geometric mean". This approach may make more sense but
  401.    it has not been tested with hardware. */
  402. #ifdef MAKE_SENSE
  403.               weight += pow(10.0, power[j].x/10.0) * (ltg[power[j].map].bark-i);
  404. #endif
  405.               power[j].x = DBMIN;
  406.            }
  407.         }   /*  check to see if the spectral line is low dB, and if  */
  408.      }      /* so replace the center of the critical band, which is */
  409.             /* the center freq. of the noise component              */
  410.  
  411. #ifdef MAKE_SENSE
  412.      if(sum <= DBMIN)  centre = (cbound[i+1]+cbound[i]) /2;
  413.      else {
  414.         index = weight/pow(10.0,sum/10.0);
  415.         centre = cbound[i] + (int) (index * (double) (cbound[i+1]-cbound[i]) );
  416.      } 
  417. #else
  418.      index = (double)( ((double)cbound[i]) * ((double)(cbound[i+1]-1)) );
  419.      centre = (int)(pow(index,0.5)+0.5);
  420. #endif
  421.  
  422.     /* locate next non-tonal component until finished; */
  423.     /* add to list of non-tonal components             */
  424. #ifdef MI_OPTION
  425.      /* Masahiro Iwadare's fix for infinite looping problem? */
  426.      if(power[centre].type == TONE) 
  427.        if (power[centre+1].type == TONE) centre++; else centre--;
  428. #else
  429.      /* Mike Li's fix for infinite looping problem */
  430.      if(power[centre].type == FALSE) centre++;
  431.  
  432.      if(power[centre].type == NOISE){
  433.        if(power[centre].x >= ltg[power[i].map].hear){
  434.          if(sum >= ltg[power[i].map].hear) sum = add_db(power[j].x,sum);
  435.          else
  436.          sum = power[centre].x;
  437.        }
  438.      }
  439. #endif
  440.      if(last == LAST) *noise = centre;
  441.      else {
  442.         power[centre].next = LAST;
  443.         power[last].next = centre;
  444.      }
  445.      power[centre].x = sum;
  446.      power[centre].type = NOISE;        
  447.      last = centre;
  448.  }        
  449. }
  450.  
  451. /****************************************************************/
  452. /*
  453. /*        This function reduces the number of noise and tonal
  454. /* component for further threshold analysis.
  455. /*
  456. /****************************************************************/
  457.  
  458. void subsampling(mask FAR power[HAN_SIZE], g_thres FAR *ltg, int *tone, int *noise)
  459. {
  460.  int i, old;
  461.  
  462.  i = *tone; old = STOP;    /* calculate tonal components for */
  463.  while(i!=LAST){           /* reduction of spectral lines    */
  464.     if(power[i].x < ltg[power[i].map].hear){
  465.        power[i].type = FALSE;
  466.        power[i].x = DBMIN;
  467.        if(old == STOP) *tone = power[i].next;
  468.        else power[old].next = power[i].next;
  469.     }
  470.     else old = i;
  471.     i = power[i].next;
  472.  }
  473.  i = *noise; old = STOP;    /* calculate non-tonal components for */
  474.  while(i!=LAST){            /* reduction of spectral lines        */
  475.     if(power[i].x < ltg[power[i].map].hear){
  476.        power[i].type = FALSE;
  477.        power[i].x = DBMIN;
  478.        if(old == STOP) *noise = power[i].next;
  479.        else power[old].next = power[i].next;
  480.     }
  481.     else old = i;
  482.     i = power[i].next;
  483.  }
  484.  i = *tone; old = STOP;
  485.  while(i != LAST){                              /* if more than one */
  486.     if(power[i].next == LAST)break;             /* tonal component  */
  487.     if(ltg[power[power[i].next].map].bark -     /* is less than .5  */
  488.        ltg[power[i].map].bark < 0.5) {          /* bark, take the   */
  489.        if(power[power[i].next].x > power[i].x ){/* maximum          */
  490.           if(old == STOP) *tone = power[i].next;
  491.           else power[old].next = power[i].next;
  492.           power[i].type = FALSE;
  493.           power[i].x = DBMIN;
  494.           i = power[i].next;
  495.        }
  496.        else {
  497.           power[power[i].next].type = FALSE;
  498.           power[power[i].next].x = DBMIN;
  499.           power[i].next = power[power[i].next].next;
  500.           old = i;
  501.        }
  502.     }
  503.     else {
  504.       old = i;
  505.       i = power[i].next;
  506.     }
  507.  }
  508. }
  509.  
  510. /****************************************************************/
  511. /*
  512. /*        This function calculates the individual threshold and
  513. /* sum with the quiet threshold to find the global threshold.
  514. /*
  515. /****************************************************************/
  516.  
  517. void threshold(mask FAR power[HAN_SIZE], g_thres FAR *ltg, int *tone, int *noise, int bit_rate)
  518. {
  519.  int k, t;
  520.  double dz, tmps, vf;
  521.  
  522.  for(k=1;k<sub_size;k++){
  523.     ltg[k].x = DBMIN;
  524.     t = *tone;          /* calculate individual masking threshold for */
  525.     while(t != LAST){   /* components in order to find the global     */
  526.        if(ltg[k].bark-ltg[power[t].map].bark >= -3.0 && /*threshold (LTG)*/
  527.           ltg[k].bark-ltg[power[t].map].bark <8.0){
  528.           dz = ltg[k].bark-ltg[power[t].map].bark; /* distance of bark value*/
  529.           tmps = -1.525-0.275*ltg[power[t].map].bark - 4.5 + power[t].x;
  530.              /* masking function for lower & upper slopes */
  531.           if(-3<=dz && dz<-1) vf = 17*(dz+1)-(0.4*power[t].x +6);
  532.           else if(-1<=dz && dz<0) vf = (0.4 *power[t].x + 6) * dz;
  533.           else if(0<=dz && dz<1) vf = (-17*dz);
  534.           else if(1<=dz && dz<8) vf = -(dz-1) * (17-0.15 *power[t].x) - 17;
  535.           tmps += vf;        
  536.           ltg[k].x = add_db(ltg[k].x, tmps);
  537.        }
  538.        t = power[t].next;
  539.     }
  540.  
  541.     t = *noise;        /* calculate individual masking threshold  */
  542.     while(t != LAST){  /* for non-tonal components to find LTG    */
  543.        if(ltg[k].bark-ltg[power[t].map].bark >= -3.0 &&
  544.           ltg[k].bark-ltg[power[t].map].bark <8.0){
  545.           dz = ltg[k].bark-ltg[power[t].map].bark; /* distance of bark value */
  546.           tmps = -1.525-0.175*ltg[power[t].map].bark -0.5 + power[t].x;
  547.              /* masking function for lower & upper slopes */
  548.           if(-3<=dz && dz<-1) vf = 17*(dz+1)-(0.4*power[t].x +6);
  549.           else if(-1<=dz && dz<0) vf = (0.4 *power[t].x + 6) * dz;
  550.           else if(0<=dz && dz<1) vf = (-17*dz);
  551.           else if(1<=dz && dz<8) vf = -(dz-1) * (17-0.15 *power[t].x) - 17;
  552.           tmps += vf;
  553.           ltg[k].x = add_db(ltg[k].x, tmps);
  554.        }
  555.        t = power[t].next;
  556.     }
  557.     if(bit_rate<96)ltg[k].x = add_db(ltg[k].hear, ltg[k].x);
  558.     else ltg[k].x = add_db(ltg[k].hear-12.0, ltg[k].x);
  559.  }
  560. }
  561.  
  562. /****************************************************************/
  563. /*
  564. /*        This function finds the minimum masking threshold and
  565. /* return the value to the encoder.
  566. /*
  567. /****************************************************************/
  568.  
  569. void II_minimum_mask(g_thres FAR *ltg,double FAR ltmin[SBLIMIT],int sblimit)
  570. {
  571.  double min;
  572.  int i,j;
  573.  
  574.  j=1;
  575.  for(i=0;i<sblimit;i++)
  576.     if(j>=sub_size-1)                   /* check subband limit, and       */
  577.        ltmin[i] = ltg[sub_size-1].hear; /* calculate the minimum masking  */
  578.     else {                              /* level of LTMIN for each subband*/
  579.        min = ltg[j].x;
  580.        while(ltg[j].line>>4 == i && j < sub_size){
  581.        if(min>ltg[j].x)  min = ltg[j].x;
  582.        j++;
  583.     }
  584.     ltmin[i] = min;
  585.  }
  586. }
  587.  
  588. /*****************************************************************/
  589. /*
  590. /*        This procedure is called in musicin to pick out the
  591. /* smaller of the scalefactor or threshold.
  592. /*
  593. /*****************************************************************/
  594.  
  595. void II_smr(double FAR ltmin[SBLIMIT], double FAR spike[SBLIMIT], double FAR scale[SBLIMIT], int sblimit)
  596. {
  597.  int i;
  598.  double max;
  599.                 
  600.  for(i=0;i<sblimit;i++){                     /* determine the signal   */
  601.     max = 20 * log10(scale[i] * 32768) - 10; /* level for each subband */
  602.     if(spike[i]>max) max = spike[i];         /* for the maximum scale  */
  603.     max -= ltmin[i];                         /* factors                */
  604.     ltmin[i] = max;
  605.  }
  606. }
  607.         
  608. /****************************************************************/
  609. /*
  610. /*        This procedure calls all the necessary functions to
  611. /* complete the psychoacoustic analysis.
  612. /*
  613. /****************************************************************/
  614.  
  615. void II_Psycho_One(short FAR buffer[2][1152], double FAR scale[2][SBLIMIT], double FAR ltmin[2][SBLIMIT], frame_params *fr_ps)
  616. {
  617.  layer *info = fr_ps->header;
  618.  int   stereo = fr_ps->stereo;
  619.  int   sblimit = fr_ps->sblimit;
  620.  int k,i, tone=0, noise=0;
  621.  static char init = 0;
  622.  static int off[2] = {256,256};
  623.  double *sample;
  624.  DSBL *spike;
  625.  static D1408 *fft_buf;
  626.  static mask_ptr FAR power;
  627.  static g_ptr FAR ltg;
  628.  
  629.  sample = (double *) mem_alloc(sizeof(DFFT), "sample");
  630.  spike = (DSBL *) mem_alloc(sizeof(D2SBL), "spike");
  631.      /* call functions for critical boundaries, freq. */
  632.  if(!init){  /* bands, bark values, and mapping */
  633.     fft_buf = (D1408 *) mem_alloc((long) sizeof(D1408) * 2, "fft_buf");
  634.     power = (mask_ptr FAR ) mem_alloc(sizeof(mask) * HAN_SIZE, "power");
  635.     read_cbound(info->lay,info->sampling_frequency);
  636.     read_freq_band(<g,info->lay,info->sampling_frequency);
  637.     make_map(power,ltg);
  638.     for (i=0;i<1408;i++) fft_buf[0][i] = fft_buf[1][i] = 0;
  639.     init = 1;
  640.  }
  641.  for(k=0;k<stereo;k++){  /* check pcm input for 3 blocks of 384 samples */
  642.     for(i=0;i<1152;i++) fft_buf[k][(i+off[k])%1408]= (double)buffer[k][i]/SCALE;
  643.     for(i=0;i<FFT_SIZE;i++) sample[i] = fft_buf[k][(i+1216+off[k])%1408];
  644.     off[k] += 1152;
  645.     off[k] %= 1408;
  646.                             /* call functions for windowing PCM samples,*/
  647.     II_hann_win(sample);    /* location of spectral components in each  */
  648.     for(i=0;i<HAN_SIZE;i++) power[i].x = DBMIN;  /*subband with labeling*/
  649.     II_f_f_t(sample, power);                     /*locate remaining non-*/
  650.     II_pick_max(power, &spike[k][0]);            /*tonal sinusoidals,   */
  651.     II_tonal_label(power, &tone);                /*reduce noise & tonal */
  652.     noise_label(power, &noise, ltg);             /*components, find     */
  653.     subsampling(power, ltg, &tone, &noise);      /*global & minimal     */
  654.     threshold(power, ltg, &tone, &noise,         /*threshold, and sgnl- */
  655.       bitrate[info->lay-1][info->bitrate_index]/stereo); /*to-mask ratio*/
  656.     II_minimum_mask(ltg, <min[k][0], sblimit);
  657.     II_smr(<min[k][0], &spike[k][0], &scale[k][0], sblimit);        
  658.  }
  659.  mem_free((void **) &sample);
  660.  mem_free((void **) &spike);
  661. }
  662.  
  663. /**********************************************************************/
  664. /*
  665. /*        This module implements the psychoacoustic model I for the
  666. /* MPEG encoder layer I. It uses simplified tonal and noise masking
  667. /* threshold analysis to generate SMR for the encoder bit allocation
  668. /* routine.
  669. /*
  670. /**********************************************************************/
  671.  
  672. /****************************************************************/
  673. /*
  674. /*        Fast Fourier transform of the input samples.
  675. /*
  676. /****************************************************************/
  677.  
  678. void I_f_f_t(double FAR sample[FFT_SIZE/2], mask FAR power[HAN_SIZE/2])         /* this function calculates */
  679. {
  680.  int i,j,k,L,l=0;
  681.  int ip, le, le1;
  682.  double t_r, t_i, u_r, u_i;
  683.  static int M, MM1, init = 0, N;
  684.  double *x_r, *x_i, *energy;
  685.  static int *rev;
  686.  static double *w_r, *w_i;
  687.  
  688.  x_r = (double *) mem_alloc(sizeof(DFFT2), "x_r");
  689.  x_i = (double *) mem_alloc(sizeof(DFFT2), "x_i");
  690.  energy = (double *) mem_alloc(sizeof(DFFT2), "energy");
  691.  for(i=0;i<FFT_SIZE/2;i++) x_r[i] = x_i[i] = energy[i] = 0;
  692.  if(!init){
  693.     rev = (int *) mem_alloc(sizeof(IFFT2), "rev");
  694.     w_r = (double *) mem_alloc(sizeof(D9), "w_r");
  695.     w_i = (double *) mem_alloc(sizeof(D9), "w_i");
  696.     M = 9;
  697.     MM1 = 8;
  698.     N = FFT_SIZE/2;
  699.     for(L=0;L<M;L++){
  700.        le = 1 << (M-L);
  701.        le1 = le >> 1;
  702.        w_r[L] = cos(PI/le1);
  703.        w_i[L] = -sin(PI/le1);
  704.     }
  705.     for(i=0;i<FFT_SIZE/2;rev[i] = l,i++) for(j=0,l=0;j<9;j++){
  706.        k=(i>>j) & 1;
  707.        l |= (k<<8-j);                
  708.     }
  709.     init = 1;
  710.  }
  711.  memcpy( (char *) x_r, (char *) sample, sizeof(double) * FFT_SIZE/2);
  712.  for(L=0;L<MM1;L++){
  713.     le = 1 << (M-L);
  714.     le1 = le >> 1;
  715.     u_r = 1;
  716.     u_i = 0;
  717.     for(j=0;j<le1;j++){
  718.        for(i=j;i<N;i+=le){
  719.           ip = i + le1;
  720.           t_r = x_r[i] + x_r[ip];
  721.           t_i = x_i[i] + x_i[ip];
  722.           x_r[ip] = x_r[i] - x_r[ip];
  723.           x_i[ip] = x_i[i] - x_i[ip];
  724.           x_r[i] = t_r;
  725.           x_i[i] = t_i;
  726.           t_r = x_r[ip];
  727.           x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
  728.           x_i[ip] = x_i[ip] * u_r + t_r * u_i;
  729.        }
  730.        t_r = u_r;
  731.        u_r = u_r * w_r[L] - u_i * w_i[L];
  732.        u_i = u_i * w_r[L] + t_r * w_i[L];
  733.     }
  734.  }
  735.  for(i=0;i<N;i+=2){
  736.     ip = i + 1;
  737.     t_r = x_r[i] + x_r[ip];
  738.     t_i = x_i[i] + x_i[ip];
  739.     x_r[ip] = x_r[i] - x_r[ip];
  740.     x_i[ip] = x_i[i] - x_i[ip];
  741.     x_r[i] = t_r;
  742.     x_i[i] = t_i;
  743.     energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
  744.  }
  745.  for(i=0;i<FFT_SIZE/2;i++) if(i<rev[i]){
  746.     t_r = energy[i];
  747.     energy[i] = energy[rev[i]];
  748.     energy[rev[i]] = t_r;
  749.  }
  750.  for(i=0;i<HAN_SIZE/2;i++){                     /* calculate power  */
  751.     if(energy[i] < 1E-20) energy[i] = 1E-20;    /* density spectrum */
  752.        power[i].x = 10 * log10(energy[i]) + POWERNORM;
  753.        power[i].next = STOP;
  754.        power[i].type = FALSE;
  755.  }
  756.  mem_free((void **) &x_r);
  757.  mem_free((void **) &x_i);
  758.  mem_free((void **) &energy);
  759. }
  760.  
  761. /****************************************************************/
  762. /*
  763. /*         Window the incoming audio signal.
  764. /*
  765. /****************************************************************/
  766.  
  767. void I_hann_win(double FAR sample[FFT_SIZE/2])             /* this function calculates a  */
  768. {                                   /* samples for a 512-pt. FFT   */
  769.  register int i;
  770.  register double sqrt_8_over_3;
  771.  static int init = 0;
  772.  static double FAR *window;
  773.  
  774.  if(!init){  /* calculate window function for the Fourier transform */
  775.     window = (double FAR *) mem_alloc(sizeof(DFFT2), "window");
  776.     sqrt_8_over_3 = pow(8.0/3.0, 0.5);
  777.     for(i=0;i<FFT_SIZE/2;i++){
  778.       /* Hann window formula */
  779.       window[i]=sqrt_8_over_3*0.5*(1-cos(2.0*PI*i/(FFT_SIZE/2)))/(FFT_SIZE/2);
  780.     }
  781.     init = 1;
  782.  }
  783.  for(i=0;i<FFT_SIZE/2;i++) sample[i] *= window[i];
  784. }
  785.  
  786. /*******************************************************************/
  787. /*
  788. /*        This function finds the maximum spectral component in each
  789. /* subband and return them to the encoder for time-domain threshold
  790. /* determination.
  791. /*
  792. /*******************************************************************/
  793. #ifndef LONDON
  794. void I_pick_max(mask FAR power[HAN_SIZE/2], double FAR spike[SBLIMIT])
  795. {
  796.  double max;
  797.  int i,j;
  798.  
  799.  /* calculate the spectral component in each subband */
  800.  for(i=0;i<HAN_SIZE/2;spike[i>>3] = max, i+=8)
  801.     for(j=0, max = DBMIN;j<8;j++) max = (max>power[i+j].x) ? max : power[i+j].x;
  802. }
  803. #else
  804. void I_pick_max(mask FAR power[HAN_SIZE], double FAR spike[SBLIMIT])
  805. {
  806.  double sum;
  807.  int i,j;
  808.  
  809.  for(i=0;i<HAN_SIZE/2;spike[i>>3] = 10.0*log10(sum), i+=8)
  810.                                                    /* calculate the      */
  811.  for(j=0, sum = pow(10.0,0.1*DBMIN);j<8;j++)       /* sum of spectral   */
  812.    sum += pow(10.0,0.1*power[i+j].x);              /* component in each  */
  813. }                                                  /* subband from bound */
  814. #endif
  815. /****************************************************************/
  816. /*
  817. /*        This function labels the tonal component in the power
  818. /* spectrum.
  819. /*
  820. /****************************************************************/
  821.  
  822. void I_tonal_label(mask FAR power[HAN_SIZE/2], int *tone)     /* this function extracts   */
  823. {
  824.  int i,j, last = LAST, first, run;
  825.  double max;
  826.  int last_but_one= LAST;
  827.  
  828.  *tone = LAST;
  829.  for(i=2;i<HAN_SIZE/2-6;i++){
  830.     if(power[i].x>power[i-1].x && power[i].x>=power[i+1].x){
  831.        power[i].type = TONE;
  832.        power[i].next = LAST;
  833.        if(last != LAST) power[last].next = i;
  834.        else first = *tone = i;
  835.        last = i;
  836.     }
  837.  }
  838.  last = LAST;
  839.  first = *tone;
  840.  *tone = LAST;
  841.  while(first != LAST){                /* conditions for the tonal     */
  842.     if(first<3 || first>250) run = 0; /* otherwise k+/-j will be out of bounds*/
  843.     else if(first<63) run = 2;        /* components in layer I, which */
  844.     else if(first<127) run = 3;       /* are the boundaries for calc.   */
  845.     else run = 6;                     /* the tonal components          */
  846.     max = power[first].x - 7;
  847.     for(j=2;j<=run;j++)  /* after calc. of tonal components, set to loc.*/
  848.        if(max < power[first-j].x || max < power[first+j].x){   /* max   */
  849.           power[first].type = FALSE;
  850.           break;
  851.        }
  852.     if(power[first].type == TONE){    /* extract tonal components */
  853.        int help=first;
  854.        if(*tone == LAST) *tone = first;
  855.        while((power[help].next!=LAST)&&(power[help].next-first)<=run)
  856.           help=power[help].next;
  857.        help=power[help].next;
  858.        power[first].next=help;
  859.        if((first-last)<=run){
  860.           if(last_but_one != LAST) power[last_but_one].next=first;
  861.        }
  862.        if(first>1 && first<255){     /* calculate the sum of the */
  863.           double tmp;                /* powers of the components */
  864.           tmp = add_db(power[first-1].x, power[first+1].x);
  865.           power[first].x = add_db(power[first].x, tmp);
  866.        }
  867.        for(j=1;j<=run;j++){
  868.           power[first-j].x = power[first+j].x = DBMIN;
  869.           power[first-j].next = power[first+j].next = STOP; /*dpwe: 2nd was .x*/
  870.           power[first-j].type = power[first+j].type = FALSE;
  871.        }
  872.        last_but_one=last;
  873.        last = first;
  874.        first = power[first].next;
  875.     }
  876.     else {
  877.        int ll;
  878.        if(last == LAST) ; /* *tone = power[first].next; dpwe */
  879.        else power[last].next = power[first].next;
  880.        ll = first;
  881.        first = power[first].next;
  882.        power[ll].next = STOP;
  883.     }
  884.  }
  885. }                        
  886.                                 
  887. /****************************************************************/
  888. /*
  889. /*        This function finds the minimum masking threshold and
  890. /* return the value to the encoder.
  891. /*
  892. /****************************************************************/
  893.  
  894. void I_minimum_mask(g_thres FAR *ltg,double FAR ltmin[SBLIMIT])
  895. {
  896.  double min;
  897.  int i,j;
  898.  
  899.  j=1;
  900.  for(i=0;i<SBLIMIT;i++)
  901.     if(j>=sub_size-1)                   /* check subband limit, and       */
  902.        ltmin[i] = ltg[sub_size-1].hear; /* calculate the minimum masking  */
  903.     else {                              /* level of LTMIN for each subband*/
  904.        min = ltg[j].x;
  905.        while(ltg[j].line>>3 == i && j < sub_size){
  906.           if (min>ltg[j].x)  min = ltg[j].x;
  907.           j++;
  908.        }
  909.        ltmin[i] = min;
  910.     }
  911. }
  912.  
  913. /*****************************************************************/
  914. /*
  915. /*        This procedure is called in musicin to pick out the
  916. /* smaller of the scalefactor or threshold.
  917. /*
  918. /*****************************************************************/
  919.  
  920. void I_smr(double FAR ltmin[SBLIMIT], double FAR spike[SBLIMIT], double FAR scale[SBLIMIT])
  921. {
  922.  int i;
  923.  double max;
  924.                 
  925.  for(i=0;i<SBLIMIT;i++){                      /* determine the signal   */
  926.     max = 20 * log10(scale[i] * 32768) - 10;  /* level for each subband */
  927.     if(spike[i]>max) max = spike[i];          /* for the scalefactor    */
  928.     max -= ltmin[i];
  929.     ltmin[i] = max;
  930.  }
  931. }
  932.         
  933. /****************************************************************/
  934. /*
  935. /*        This procedure calls all the necessary functions to
  936. /* complete the psychoacoustic analysis.
  937. /*
  938. /****************************************************************/
  939.  
  940. void I_Psycho_One(short FAR buffer[2][1152], double FAR scale[2][SBLIMIT], double FAR ltmin[2][SBLIMIT], frame_params *fr_ps)
  941. {
  942.  int stereo = fr_ps->stereo;
  943.  the_layer info = fr_ps->header;
  944.  int k,i, tone=0, noise=0;
  945.  static char init = 0;
  946.  static int off[2] = {256,256};
  947.  double *sample;
  948.  DSBL *spike;
  949.  static D640 *fft_buf;
  950.  static mask_ptr FAR power;
  951.  static g_ptr FAR ltg;
  952.  
  953.  sample = (double *) mem_alloc(sizeof(DFFT2), "sample");
  954.  spike = (DSBL *) mem_alloc(sizeof(D2SBL), "spike");
  955.             /* call functions for critical boundaries, freq. */
  956.  if(!init){ /* bands, bark values, and mapping              */
  957.     fft_buf = (D640 *) mem_alloc(sizeof(D640) * 2, "fft_buf");
  958.     power = (mask_ptr FAR ) mem_alloc(sizeof(mask) * HAN_SIZE/2, "power");
  959.     read_cbound(info->lay,info->sampling_frequency);
  960.     read_freq_band(<g,info->lay,info->sampling_frequency);
  961.     make_map(power,ltg);
  962.     for(i=0;i<640;i++) fft_buf[0][i] = fft_buf[1][i] = 0;
  963.     init = 1;
  964.  }
  965.  for(k=0;k<stereo;k++){    /* check PCM input for a block of */
  966.     for(i=0;i<384;i++)     /* 384 samples for a 512-pt. FFT  */
  967.        fft_buf[k][(i+off[k])%640]= (double) buffer[k][i]/SCALE;
  968.     for(i=0;i<FFT_SIZE/2;i++)
  969.        sample[i] = fft_buf[k][(i+448+off[k])%640];
  970.     off[k] += 384;
  971.     off[k] %= 640;
  972.                         /* call functions for windowing PCM samples,   */
  973.     I_hann_win(sample); /* location of spectral components in each     */
  974.     for(i=0;i<HAN_SIZE/2;i++) power[i].x = DBMIN;   /* subband with    */
  975.     I_f_f_t(sample, power);              /* labeling, locate remaining */
  976.     I_pick_max(power, &spike[k][0]);     /* non-tonal sinusoidals,     */
  977.     I_tonal_label(power, &tone);         /* reduce noise & tonal com., */
  978.     noise_label(power, &noise, ltg);     /* find global & minimal      */
  979.     subsampling(power, ltg, &tone, &noise);  /* threshold, and sgnl-   */
  980.     threshold(power, ltg, &tone, &noise,     /* to-mask ratio          */
  981.       bitrate[info->lay-1][info->bitrate_index]/stereo);
  982.     I_minimum_mask(ltg, <min[k][0]);
  983.     I_smr(<min[k][0], &spike[k][0], &scale[k][0]);        
  984.  }
  985.  mem_free((void **) &sample);
  986.  mem_free((void **) &spike);
  987. }
  988.