开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > mpeg/mp3 > 查看源码
TONAL.C
资源名称：mpeg.zip　[点击查看]
上传用户：njqiyou
上传日期：2007-01-08
资源大小：574k
文件大小：47k
源码类别：

mpeg/mp3

开发平台：

C/C++

TONAL.C：源码内容
  1. /**********************************************************************
  2.  * ISO MPEG Audio Subgroup Software Simulation Group (1996)
  3.  * ISO 13818-3 MPEG-2 Audio Multichannel Encoder
  4.  *
  5.  * $Id: tonal.c 1.8 1996/02/12 07:13:35 rowlands Exp $
  6.  *
  7.  * $Log: tonal.c $
  8.  * Revision 1.8  1996/02/12 07:13:35  rowlands
  9.  * Release following Munich meeting
  10.  *
  11.  * Revision 1.5.2.1  1995/11/06  04:19:12  rowlands
  12.  * Received from Uwe Felderhoff (IRT)
  13.  *
  14.  * Revision 1.7  1995/08/14  08:06:37  tenkate
  15.  * ML-LSF added Warner ten Kate 7/8/95 (Philips)
  16.  * II_psycho_one() split into II_psycho_one() for audio data and
  17.  * II_psycho_one_ml() for MultiLingual data.
  18.  * Variables crit_band, cbound and sub_size has been made local to
  19.  * these functions.
  20.  * Tables "cb" and "th" are copied from LSF-directory.
  21.  * ltg has found its ml counterpart.
  22.  *
  23.  * Revision 1.6  1995/07/31  07:47:50  tenkate
  24.  * bugs correction (if center==1 in Psycho_One), 25/07/95 WtK
  25.  *
  26.  * Revision 1.4.2.1  1995/06/16  03:46:42  rowlands
  27.  * Input from Susanne Ritscher (IRT)
  28.  *
  29.  **********************************************************************/
  31. /**********************************************************************
  32.  *   date   programmers         comment                               *
  33.  * 2/25/91  Douglas Wong        start of version 1.1 records          *
  34.  * 3/06/91  Douglas Wong        rename: setup.h to endef.h            *
  35.  *                              updated I_psycho_one and II_psycho_one*
  36.  * 3/11/91  W. J. Carter        Added Douglas Wong's updates dated    *
  37.  *                              3/9/91 for I_Psycho_One() and for     *
  38.  *                              II_Psycho_One().                      *
  39.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  40.  *                              Located and fixed numerous software   *
  41.  *                              bugs and table data errors.           *
  42.  * 6/11/91  Davis Pan           corrected several bugs                *
  43.  *                              based on comments from H. Fuchs       *
  44.  * 01jul91  dpwe (Aware Inc.)   Made pow() args float                 *
  45.  *                              Removed logical bug in I_tonal_label: *
  46.  *                              Sometimes *tone returned == STOP      *
  47.  * 7/10/91  Earle Jennings      no change necessary in port to MsDos  *
  48.  * 11sep91  dpwe@aware.com      Subtracted 90.3dB from II_f_f_t peaks *
  49.  * 10/1/91  Peter W. Farrett    Updated II_Psycho_One(),I_Psycho_One()*
  50.  * to include comments.       *
  51.  *11/29/91  Masahiro Iwadare    Bug fix regarding POWERNORM           *
  52.  *                              fixed several other miscellaneous bugs*
  53.  * 2/11/92  W. Joseph Carter    Ported new code to Macintosh.  Most   *
  54.  *                              important fixes involved changing     *
  55.  *                              16-bit ints to long or unsigned in    *
  56.  *                              bit alloc routines for quant of 65535 *
  57.  *                              and passing proper function args.     *
  58.  *                              Removed "Other Joint Stereo" option   *
  59.  *                              and made bitrate be total channel     *
  60.  *                              bitrate, irrespective of the mode.    *
  61.  *                              Fixed many small bugs & reorganized.  *
  62.  * 2/12/92  Masahiro Iwadare    Fixed some potential bugs in          *
  63.  *          Davis Pan           subsampling()                         *
  64.  * 2/25/92  Masahiro Iwadare    Fixed some more potential bugs        *
  65.  * 92-11-06 Soren H. Nielsen Corrected power calculation in I_ and *
  66.  * II_f_f_t.               *
  67.  **********************************************************************
  68.  *                                                                    *
  69.  *                                                                    *
  70.  *  MPEG/audio Phase 2 coding/decoding multichannel                   *
  71.  *                                                                    *
  72.  *  7/27/93        Susanne Ritscher,  IRT Munich                      *
  73.  *  8/27/93        Susanne Ritscher, IRT Munich                       *
  74.  *                 Channel-Switching is working                       *
  75.  *  9/1/93         Susanne Ritscher,  IRT Munich                      *
  76.  *                 all channels normalized                            *
  77.  *                                                                    *
  78.  *  9/20/93        channel-switching is only performed at a           *
  79.  *                 certain limit of TC_ALLOC dB, which is included    *
  80.  *                 in encoder.h                                       *
  81.  *                                                                    *
  82.  * 10/18/93        seperated smr and ltmin                            *
  83.  *                                                                    *
  84.  *  Version 1.0                                                       *
  85.  *                                                                    *
  86.  *  07/12/94       Susanne Ritscher,  IRT Munich                      *
  87.  *                 Tel: +49 89 32399 458                              *
  88.  *                 Fax: +49 89 32399 415                              *
  89.  *                                                                    *
  90.  *  Version 1.1                                                       *
  91.  *                                                                    *
  92.  *  02/23/95    Susanne Ritscher,  IRT Munich                      *
  93.  *                 corrected some bugs                                *
  94.  *                 extension bitstream is working                     *
  95.  *                                                                    *
  96.  **********************************************************************/
  97. #define TONAL_WIDTH 0.1 /* When more than one tonal component is within
  98.                            this width in Bark, the weaker one(s) are
  99.                            eliminated */
  102. #include "common.h"
  103. #include "encoder.h"
  105. /**********************************************************************/
  106. /*
  107. /*        This module implements the psychoacoustic model I for the
  108. /* MPEG encoder layer II. It uses simplified tonal and noise masking
  109. /* threshold analysis to generate SMR for the encoder bit allocation
  110. /* routine.
  111. /*
  112. /**********************************************************************/
  114. int read_crit_band(int lay, int freq) 
  115. {
  116.  int crit_band;
  117.  FILE *fp;
  118.  char r[16], t[80];
  120.  strcpy(r, "2cb1");
  121.  r[0] = (char) lay + '0';
  122.  r[3] = (char) freq + '0';
  123.  if( !(fp = OpenTableFile(r)) ){       /* check boundary values */
  124.     printf("Please check %s boundary tablen",r);
  125.     exit(0);
  126.  }
  127.  fgets(t,80,fp);
  128.  sscanf(t,"%dn",&crit_band);
  129.  fclose(fp);
  130.  return(crit_band);
  131. }        
  133. void read_cbound(int lay, int freq, int crit_band, int *cbound) 
  134.  /* this function reads in critical band boundaries */
  135. {
  136.  int i,j,k;
  137.  FILE *fp;
  138.  char r[16], t[80];
  140.  strcpy(r, "2cb1");
  141.  r[0] = (char) lay + '0';
  142.  r[3] = (char) freq + '0';
  143.  if( !(fp = OpenTableFile(r)) ){       /* check boundary values */
  144.     printf("Please check %s boundary tablen",r);
  145.     exit(0);
  146.  }
  147.  fgets(t,80,fp);               /* skip input for critical bands */
  148.  sscanf(t,"%dn",&k);
  149.  for(i=0;i<crit_band;i++){ 
  150.     fgets(t,80,fp);
  151.     sscanf(t,"%d %dn",&j, &k);
  152.     if(i==j) cbound[j] = k;
  153.     else {                     /* error */
  154.        printf("Please check index %d in cbound table %sn",i,r);
  155.        exit(0);
  156.     }
  157.  }
  158.  fclose(fp);
  159. }        
  161. void read_freq_band(int *sub_size, g_ptr *ltg, int lay, int freq) 
  162.   /* this function reads in frequency bands and bark values  */
  163. {
  164.  int i,j, k;
  165.  double a,b,c;
  166.  FILE *fp;
  167.  char r[16], t[80];
  169.  strcpy(r, "2th1");
  170.  r[0] = (char) lay + '0';
  171.  r[3] = (char) freq + '0';
  172.  if( !(fp = OpenTableFile(r)) ){   /* check freq. values  */
  173.     printf("Please check frequency and cband table %sn",r);
  174.     exit(0);
  175.  }
  176.  fgets(t,80,fp);              /* read input for freq. subbands */
  177.  sscanf(t,"%dn",sub_size);
  178.  *ltg = (g_ptr) mem_alloc (sizeof(g_thres) * (*sub_size), "ltg");
  179.  (*ltg)[0].line = 0;          /* initialize global masking threshold */
  180.  (*ltg)[0].bark = 0;
  181.  (*ltg)[0].hear = 0;
  182.  for(i=1;i<(*sub_size);i++){    /* continue to read freq. subband */
  183.     fgets(t,80,fp);          /* and assign                     */
  184.     sscanf(t,"%d %d %lf %lfn",&j, &k, &b, &c);
  185.     if(i == j){
  186.        (*ltg)[j].line = k;
  187.        (*ltg)[j].bark = b;
  188.        (*ltg)[j].hear = c;
  189.     }
  190.     else {                   /* error */
  191.        printf("Please check index %d in freq-cb table %sn",i,r);
  192.        exit(0);
  193.     }
  194.  }
  195.  fclose(fp);
  196. }
  198. void make_map (int sub_size, mask *power, g_thres *ltg)
  199. /* this function calculates the global masking threshold */
  200. {
  201.     int i, j;
  202.    
  203.     for (i = 1; i < sub_size; i++)
  204. for (j = ltg[i-1].line; j <= ltg[i].line; j++)
  205.     power[j].map = i;
  206. }
  208. double add_db(double a, double b)
  209. {
  210.  a = pow(10.0,a/10.0);
  211.  b = pow(10.0,b/10.0);
  212.  return 10 * log10(a+b);
  213. }
  215. /****************************************************************/
  216. /*
  217. /*        Fast Fourier transform of the input samples.
  218. /*
  219. /****************************************************************/
  221. void II_f_f_t (double *sample, mask *power) /* this function calculates an */
  222.     /* FFT analysis for the freq.  */
  223.     /* domain                      */
  224. {
  225.     int i,j,k,ll,l=0;
  226.     int ip, le, le1;
  227.     double t_r, t_i, u_r, u_i;
  228.     static int M, MM1, init = 0, N, NV2, NM1;
  229.     double *x_r, *x_i, *energy;
  230.     static int *rev;
  231.     static double *w_r, *w_i;
  233.     x_r = (double *) mem_alloc (sizeof (DFFT), "x_r");
  234.     x_i = (double *) mem_alloc (sizeof (DFFT), "x_i");
  235.     energy = (double *) mem_alloc (sizeof (DFFT), "energy");
  237.     if (!init)
  238.     {
  239. rev = (int *) mem_alloc (sizeof (IFFT), "rev");
  240. w_r = (double *) mem_alloc (sizeof (D10), "w_r");
  241. w_i = (double *) mem_alloc (sizeof (D10), "w_i");
  242. M = 10;
  243. MM1 = 9;
  244. N = FFT_SIZE;
  245. NV2 = FFT_SIZE >> 1;
  246. NM1 = FFT_SIZE - 1;
  247. for (ll = 0; ll < M; ll++)
  248. {
  249.     le = 1 << (M-ll);
  250.     le1 = le >> 1;
  251.     w_r[ll] = cos (PI / le1);
  252.     w_i[ll] = -sin (PI / le1);
  253. }
  254. for (i = 0; i < FFT_SIZE; rev[i] = l, i++)
  255.     for (j = 0, l = 0; j < 10; j++)
  256.     {
  257. k = (i >> j) & 1;
  258. l |= (k << 9 - j);                
  259.     }
  260. init = 1;
  261.     }
  263.     for (i = 0; i < FFT_SIZE; i++)
  264.     {
  265.        x_r[i] = sample[i];
  266.        x_i[i] = energy[i] = 0;
  267.     }
  268. /*
  269.     memcpy ((char *) x_r, (char *) sample, sizeof (double) * FFT_SIZE);
  270. */
  271.     for (ll = 0; ll < MM1; ll++)
  272.     {
  273. le = 1 << (M - ll);
  274. le1 = le >> 1;
  275. u_r = 1;
  276. u_i = 0;
  277. for (j = 0; j < le1; j++)
  278. {
  279.     for (i = j; i < N; i += le)
  280.     {
  281. ip = i + le1;
  282. t_r = x_r[i] + x_r[ip];
  283. t_i = x_i[i] + x_i[ip];
  284. x_r[ip] = x_r[i] - x_r[ip];
  285. x_i[ip] = x_i[i] - x_i[ip];
  286. x_r[i] = t_r;
  287. x_i[i] = t_i;
  288. t_r = x_r[ip];
  289. x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
  290. x_i[ip] = x_i[ip] * u_r + t_r * u_i;
  291.     }
  292.     t_r = u_r;
  293.     u_r = u_r * w_r[ll] - u_i * w_i[ll];
  294.     u_i = u_i * w_r[ll] + t_r * w_i[ll];
  295. }
  296.     }
  297.     for (i = 0; i < N; i += 2)
  298.     {
  299. ip = i + 1;
  300. t_r = x_r[i] + x_r[ip];
  301. t_i = x_i[i] + x_i[ip];
  302. x_r[ip] = x_r[i] - x_r[ip];
  303. x_i[ip] = x_i[i] - x_i[ip];
  304. x_r[i] = t_r;
  305. x_i[i] = t_i;
  306. energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
  307.     }
  308.     for (i = 0; i < FFT_SIZE; i++)
  309. if (i < rev[i])
  310. {
  311.     t_r = energy[i];
  312.     energy[i] = energy[rev[i]];
  313.     energy[rev[i]] = t_r;
  314. }
  315.     for (i = 0; i < HAN_SIZE; i++)
  316.     {
  317. /* calculate power density spectrum */
  318. if (energy[i] < 1E-20)
  319.     energy[i] = 1E-20;
  320. /* power calculation corrected with a factor 4, both positive
  321.    and negative frequencies exist, 1992-11-06 shn */
  322. power[i].x = 10 * log10 (energy[i] * 4.0) + POWERNORM;
  323. power[i].next = STOP;
  324. power[i].type = FALSE;
  325.     }
  327.     mem_free ((void **) &x_r);
  328.     mem_free ((void **) &x_i);
  329.     mem_free ((void **) &energy);
  331. }
  333. /****************************************************************/
  334. /*
  335. /*         Window the incoming audio signal.
  336. /*
  337. /****************************************************************/
  339. void II_hann_win (double *sample)   /* this function calculates a  */
  340.     /* Hann window for PCM (input) */
  341. {     /* samples for a 1024-pt. FFT  */
  342.     register int i;
  343.     register double sqrt_8_over_3;
  344.     static int init = 0;
  345.     static double *window;
  346.     
  347.     if (!init)
  348.     {
  349. /* calculate window function for the Fourier transform */
  350. window = (double *) mem_alloc (sizeof (DFFT), "window");
  351. sqrt_8_over_3 = pow (8.0 / 3.0, 0.5);
  352. for (i = 0; i < FFT_SIZE; i++)
  353. {
  354.     /* Hann window formula */
  355.     window[i] = sqrt_8_over_3 * 0.5 * (1 - cos (2.0*PI*i / (FFT_SIZE-1))) / FFT_SIZE;
  356. }
  357. init = 1;
  358.     }
  360.     for (i = 0; i < FFT_SIZE; i++)
  361. sample[i] *= window[i];
  362. }
  364. /*******************************************************************/
  365. /*
  366. /*        This function finds the maximum spectral component in each
  367. /* subband and return them to the encoder for time-domain threshold
  368. /* determination.
  369. /*
  370. /*******************************************************************/
  372. void II_pick_max (mask *power, double *spike)
  373. {
  374.  double max;
  375.  int i,j;
  377.  for(i=0;i<HAN_SIZE;spike[i>>4] = max, i+=16)      /* calculate the      */
  378.  for(j=0, max = DBMIN;j<16;j++)                    /* maximum spectral   */
  379.     max = (max>power[i+j].x) ? max : power[i+j].x; /* component in each  */
  380. }                                                  /* subband from bound */
  381.                                                    /* 4-16               */
  383. /****************************************************************/
  384. /*
  385. /*        This function labels the tonal component in the power
  386. /* spectrum.
  387. /*
  388. /****************************************************************/
  390. void II_tonal_label(mask *power, int *tone) /* this function extracts (tonal) */
  391.     /* sinusoidals from the spectrum  */
  392. {
  393.  int i,j, last = LAST, first, run, last_but_one = LAST; /* dpwe */
  394.  double max;
  396.  *tone = LAST;
  397.  for(i=2;i<HAN_SIZE-12;i++){
  398.     if(power[i].x>power[i-1].x && power[i].x>=power[i+1].x){
  399.        power[i].type = TONE;
  400.        power[i].next = LAST;
  401.        if(last != LAST) power[last].next = i;
  402.        else first = *tone = i;
  403.        last = i;
  404.     }
  405.  }
  406.  last = LAST;
  407.  first = *tone;
  408.  *tone = LAST;
  409.  while(first != LAST){               /* the conditions for the tonal          */
  410.     if(first<2 || first>499) run = 0;/* otherwise k+/-j will be out of bounds */
  411.     else if(first<62) run = 2;       /* components in layer II, which         */
  412.     else if(first<126) run = 3;      /* are the boundaries for calc.          */
  413.     else if(first<254) run = 6;      /* the tonal components                  */
  414.     else run = 12;
  415.     max = power[first].x - 7;        /* after calculation of tonal   */
  416.     for(j=2;j<=run;j++)              /* components, set to local max */
  417.        if(max < power[first-j].x || max < power[first+j].x){
  418.           power[first].type = FALSE;
  419.           break;
  420.        }
  421.     if(power[first].type == TONE){   /* extract tonal components */
  422.        int help=first;
  423.        if(*tone==LAST) *tone = first;
  424.        while((power[help].next!=LAST)&&(power[help].next-first)<=run)
  425.           help=power[help].next;
  426.        help=power[help].next;
  427.        power[first].next=help;
  428.        if((first-last)<=run){
  429.           if(last_but_one != LAST) power[last_but_one].next=first;
  430.        }
  431.        if(first>1 && first<255){     /* calculate the sum of the */
  432.           double tmp;                /* powers of the components */
  433.           tmp = add_db(power[first-1].x, power[first+1].x);
  434.           power[first].x = add_db(power[first].x, tmp);
  435.        }
  436.        for(j=1;j<=run;j++){
  437.           power[first-j].x = power[first+j].x = DBMIN;
  438.           power[first-j].next = power[first+j].next = STOP;
  439.           power[first-j].type = power[first+j].type = FALSE;
  440.        }
  441.        last_but_one=last;
  442.        last = first;
  443.        first = power[first].next;
  444.     }
  445.     else {
  446.        int ll;
  447.        if(last == LAST); /* *tone = power[first].next; dpwe */
  448.        else power[last].next = power[first].next;
  449.        ll = first;
  450.        first = power[first].next;
  451.        power[ll].next = STOP;
  452.     }
  453.  }
  454. }
  456. /****************************************************************/
  457. /*
  458. /*        This function groups all the remaining non-tonal
  459. /* spectral lines into critical band where they are replaced by
  460. /* one single line.
  461. /*
  462. /****************************************************************/
  463.         
  464. void noise_label(int crit_band, int *cbound, mask *power, int *noise, g_thres *ltg)
  465. {
  466.  int i,j, centre, last = LAST;
  467.  double index, weight, sum;
  468.                               /* calculate the remaining spectral */
  469.  for(i=0;i<crit_band-1;i++){  /* lines for non-tonal components   */
  470.      for(j=cbound[i],weight = 0.0,sum = DBMIN;j<cbound[i+1];j++){
  471.         if(power[j].type != TONE){
  472.            if(power[j].x != DBMIN){
  473.               sum = add_db(power[j].x,sum);
  474.               weight += pow(10.0, power[j].x/10.0) * (ltg[power[j].map].bark-i);
  475.               power[j].x = DBMIN;
  476.            }
  477.         }   /* check to see if the spectral line is low dB, and if  */
  478.      }      /* so replace the center of the critical band, which is */
  479.             /* the center freq. of the noise component              */
  480.      if(sum <= DBMIN)  centre = (cbound[i+1]+cbound[i]) /2;
  481.      else {
  482.         index = weight/pow(10.0,sum/10.0);
  483.         centre = cbound[i] + (int) (index * (double) (cbound[i+1]-cbound[i]) );
  484.      }     /* locate next non-tonal component until finished; */
  485.            /* add to list of non-tonal components             */
  486.      if(power[centre].type == TONE) centre++;
  487.      if(last == LAST) *noise = centre;
  488.      else {
  489.         power[centre].next = LAST;
  490.         power[last].next = centre;
  491.      }
  492.      power[centre].x = sum;
  493.      power[centre].type = NOISE;        
  494.      last = centre;
  495.  }        
  496. }
  498. /****************************************************************/
  499. /*
  500. /*        This function reduces the number of noise and tonal
  501. /* component for further threshold analysis.
  502. /*
  503. /****************************************************************/
  505. void subsampling (mask *power, g_thres *ltg, int *tone, int *noise)
  506. {
  507.     int i, old;
  508.    
  509.     i = *tone;
  510.     old = STOP;
  511.     /* calculate tonal components for */
  512.     while (i != LAST)
  513.     {
  514. /* reduction of spectral lines    */
  515. if (power[i].x < ltg[power[i].map].hear)
  516. {
  517.     power[i].type = FALSE;
  518.     power[i].x = DBMIN;
  519.     if (old == STOP)
  520. *tone = power[i].next;
  521.     else
  522. power[old].next = power[i].next;
  523. }
  524. else
  525.     old = i;
  526. i = power[i].next;
  527.     }
  528.     i = *noise;
  529.     old = STOP;
  530.     /* calculate non-tonal components for */
  531.     while (i != LAST)
  532.     {
  533. /* reduction of spectral lines        */
  534.         if (power[i].x < ltg[power[i].map].hear)
  535. {
  536.     power[i].type = FALSE;
  537.     power[i].x = DBMIN;
  538.     if (old == STOP)
  539. *noise = power[i].next;
  540.     else
  541. power[old].next = power[i].next;
  542. }
  543. else
  544.     old = i;
  545. i = power[i].next;
  546.     }
  547.     i = *tone;
  548.     old = STOP;
  549.     while (i != LAST)
  550.     {
  551. /* if more than one */
  552. if (power[i].next == LAST)
  553.     break;             /* tonal component  */
  554. if (ltg[power[power[i].next].map].bark -     /* is less than .5  */
  555.     ltg[power[i].map].bark < TONAL_WIDTH)
  556. {
  557.     /* bark, take the   */
  558.     if (power[power[i].next].x > power[i].x)
  559.     {
  560. /* maximum          */
  561. if (old == STOP)
  562.     *tone = power[i].next;
  563. else
  564.     power[old].next = power[i].next;
  565. power[i].type = FALSE;
  566. power[i].x = DBMIN;
  567.     }
  568.     else
  569.     {
  570. power[power[i].next].type = FALSE;
  571. power[power[i].next].x = DBMIN;
  572. power[i].next = power[power[i].next].next;
  573.     }
  574. }
  575. else
  576.     old = i;
  577. i = power[i].next;
  578.     }
  579. }
  582. /* ----------------------------------------------------------------------------
  583. The masking function parameters are here set according to the parameters in
  584. the IRT real time implementation. The constant definitions are for convenience.
  585. 1993-07-23 shn
  586. ---------------------------------------------------------------------------- */
  588. #define AV_TONAL_K    -9.0 /* Masking index, tonal, constant part  [dB] */
  589. #define AV_NOISE_K    -5.0 /* Masking index, noisy, constant part  [dB] */
  590. #define AV_TONAL_DZ   -0.3 /* Masking index, tonal, CBR dependence [dB/Bark] */
  591. #define AV_NOISE_DZ   -0.3 /* Masking index, noisy, CBR dependence [dB/Bark] */
  593. #define LOW_LIM_1     -1.0 /* 1st lower slope from 0         to LOW_LIM_1 [Bark] */
  594. #define LOW_LIM_2     -3.0 /* 2nd lower slope from LOW_LIM_1 to LOW_LIM_2 [Bark] */
  596. #define LOW_DZ_K_1     6.0 /* 1st lower slope, constant part [dB/Bark] */
  597. #define LOW_DZ_SPL_1   0.4 /* 1st lower slope, SPL dependence [dB/(Bark*dB)] */
  598. #define LOW_DZ_MIN_1  17.0 /* 1st lower slope, minimum value [dB/Bark] */
  599. #define LOW_DZ_2      17.0 /* 2nd lower slope [dB/Bark] */
  601. #define UP_LIM_1       1.0 /* 1st upper slope from 0        to UP_LIM_1 [Bark] */
  602. #define UP_LIM_2       8.0 /* 2nd upper slope from UP_LIM_1 to UP_LIM_2 [Bark] */
  604. #define UP_DZ_1      -18.0 /* 1st upper slope, constant part [dB/Bark] */
  605. #define UP_SPL_1       0.0 /* 1st upper slope, SPL dependence [dB/(Bark*dB)] */
  606. #define UP_DZ_2      -17.0 /* 2nd upper slope, constant part [dB/Bark] */
  607. #define UP_SPL_2      -0.1 /* 2nd upper slope, SPL dependence [dB/(Bark*dB)] */
  609. #define H_THR_OFFSET -12.0 /* Hearing threshold offset [dB] */
  610. #define H_THR_OS_BR    0 /*96*/   /* Minimum datarate for offset, [kbit/s per channel] */
  612. #define MASK_ADD       2.0 /* Addition of maskers                 [dB] */
  613. #define QUIET_ADD      3.0 /* Addition of masker and threshold in quiet [dB] */
  615. /****************************************************************
  616. *
  617. *        This function calculates the individual threshold and
  618. * sum with the quiet threshold to find the global threshold.
  619. *
  620. ****************************************************************/
  622. void threshold(int sub_size, mask *power, g_thres *ltg, int *tone, int *noise, int bit_rate)
  623. {
  624.   int k, t;
  625.   double z, dz, spl, vf, tmps;
  627.   for (k=1; k<sub_size; k++) /* Target frequencies */
  628.   {
  629.     ltg[k].x = DBMIN;
  630.     t = *tone;          /* calculate individual masking threshold  */
  631.     while(t != LAST)    /* for tonal components, t,  to find LTG   */
  632.     {
  633.       z   = ltg[power[t].map].bark; /* critical band rate of masker */
  634.       dz  = ltg[k].bark - z ;       /* distance of bark value*/
  635.       spl = power[t].x;             /* sound pressure level of masker */
  637.       if (dz >= LOW_LIM_2 && dz <  UP_LIM_2)
  638.       {
  639.         tmps = spl + AV_TONAL_K + AV_TONAL_DZ * z;
  641.         /* masking function for lower & upper slopes */
  642.         if (LOW_LIM_2<=dz && dz<LOW_LIM_1)
  643.           if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
  644.             vf = LOW_DZ_2 * (dz - LOW_LIM_1) + 
  645.                 (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * LOW_LIM_1;
  646.           else
  647.             vf = LOW_DZ_2 * (dz - LOW_LIM_1) + LOW_DZ_MIN_1 * LOW_LIM_1;
  648.         else if (LOW_LIM_1<=dz && dz<0)
  649.           if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
  650.             vf = (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * dz;
  651.           else
  652.             vf = LOW_DZ_MIN_1 * dz;
  653.         else if (0<=dz && dz<UP_LIM_1)
  654.           vf = (UP_DZ_1 * dz);
  655.         else if (UP_LIM_1<=dz && dz<UP_LIM_2)
  656.           vf = (dz - UP_LIM_1) * (UP_DZ_2 - UP_SPL_2 * spl) +
  657.                UP_DZ_1 * UP_LIM_1;
  658.         tmps += vf;        
  659.         ltg[k].x = non_lin_add(ltg[k].x, tmps, MASK_ADD);
  660.       }
  661.       t = power[t].next;
  662.     } /* while */
  664.     t = *noise;        /* calculate individual masking threshold   */
  665.     while(t != LAST)   /* for non-tonal components, t, to find LTG */
  666.     {
  667.       z   = ltg[power[t].map].bark; /* critical band rate of masker */
  668.       dz  = ltg[k].bark - z ;       /* distance of bark value*/
  669.       spl = power[t].x;             /* sound pressure level of masker */
  671.       if (dz >= LOW_LIM_2 && dz <  UP_LIM_2)
  672.       {
  673.         tmps = spl + AV_NOISE_K + AV_NOISE_DZ * z;
  675.         /* masking function for lower & upper slopes */
  676.         if (LOW_LIM_2<=dz && dz<LOW_LIM_1)
  677.           if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
  678.             vf = LOW_DZ_2 * (dz - LOW_LIM_1) + 
  679.                 (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * LOW_LIM_1;
  680.           else
  681.             vf = LOW_DZ_2 * (dz - LOW_LIM_1) + LOW_DZ_MIN_1 * LOW_LIM_1;
  682.         else if (LOW_LIM_1<=dz && dz<0)
  683.           if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
  684.             vf = (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * dz;
  685.           else
  686.             vf = LOW_DZ_MIN_1 * dz;
  687.         else if (0<=dz && dz<UP_LIM_1)
  688.           vf = (UP_DZ_1 * dz);
  689.         else if (UP_LIM_1<=dz && dz<UP_LIM_2)
  690.           vf = (dz - UP_LIM_1) * (UP_DZ_2 - UP_SPL_2 * spl) +
  691.                UP_DZ_1 * UP_LIM_1;
  692.        tmps += vf;        
  693.        ltg[k].x = non_lin_add(ltg[k].x, tmps, MASK_ADD);
  694.      }
  695.       t = power[t].next;
  696.     }
  698.     if (bit_rate < H_THR_OS_BR)
  699.       ltg[k].x = non_lin_add(ltg[k].hear, ltg[k].x, QUIET_ADD);
  700.     else
  701.       ltg[k].x = non_lin_add(ltg[k].hear + H_THR_OFFSET, ltg[k].x, QUIET_ADD);
  702.  
  703.   } /* for */
  704.   fflush(stderr);
  705. }
  709. /* --------------------------------------------------------------
  710. non_lin_add
  711. A flexible addition function for levels.
  712. Input: a,b: the levels to be added.
  713.  c: the number of dB increase when a and b are equal.
  714. Common values for c are 3.01 (power addition)
  715.     and 6.02 (voltage addition).
  716. 10.0/(10*log10(2)) = 3.3219
  717. Function added 1993-04-14 Soren H. Nielsen
  718. -------------------------------------------------------------- */
  719. double non_lin_add(double a, double b, double c)
  720. {
  721.   c *= 3.3219;
  722.   a = pow(10.0, a/c); b = pow(10.0, b/c);
  723.   return(c*log10(a+b));
  724. }
  727. /****************************************************************/
  728. /*
  729. /*        This function finds the minimum masking threshold and
  730. /* return the value to the encoder.
  731. /*
  732. /****************************************************************/
  734. void II_minimum_mask(int sub_size, g_thres *ltg, double *ltmin, int sblimit)
  735. {
  736.  double min;
  737.  int i,j;
  738.  
  740.  j=1;
  741.  for(i=0;i<sblimit;i++)
  742.     if(j>=sub_size-1)                   /* check subband limit, and       */
  743.        ltmin[i] = ltg[sub_size-1].hear; /* calculate the minimum masking  */
  744.     else {                              /* level of LTMIN for each subband*/
  745.        min = ltg[j].x;
  746.        while(ltg[j].line>>4 == i && j < sub_size){
  747.        if(min>ltg[j].x)  min = ltg[j].x;
  748.        j++;
  749.     }
  750.     ltmin[i] = min;
  751.  }
  752. }
  754. /*****************************************************************/
  755. /*
  756. /*        This procedure is called in musicin to pick out the
  757. /* smaller of the scalefactor or threshold.
  758. /*
  759. /*****************************************************************/
  761. void II_smr (double *ltmin, double *smr, double *spike, double *scale, int sblimit, int l, int m)
  762. {
  763.     int i,j;
  764.     double max;
  765.    
  766.     for (i = l; i < m; i++)
  767.     {
  768. /* determine the signal   */
  769. max = 20 * log10 (scale[i] * 32768) - 10;   /* level for each subband */
  770. if (spike[i] > max)
  771.     max = spike[i];     /* for the maximum scale  */
  772. max -= ltmin[i];     /* factors                */
  773. smr[i] = max;
  774.     }
  775. }
  776.         
  777. /****************************************************************/
  778. /*
  779. /*        This procedure calls all the necessary functions to
  780. /* complete the psychoacoustic analysis.
  781. /*
  782. /****************************************************************/
  784. void II_Psycho_One (double (*buffer)[1152],
  785.     double (*scale)[32],
  786.     double (*ltmin)[32],
  787.     frame_params *fr_ps,
  788.     double (*smr)[32],
  789.     double (*spiki)[32],
  790.     int aiff)
  791. {
  792.     layer *info = fr_ps->header;
  793.     int   stereo = fr_ps->stereo;
  794.     int   stereomc = fr_ps->stereomc;
  795.     int   stereoaug = fr_ps->stereoaug;
  796.     int   sblimit = fr_ps->sblimit;
  797.     int k,i, tone=0, noise=0;
  798.     static char init = 0;
  799.     static int  off[12] = {256, 256, 256, 256, 256, 256, 256, 256, 256, 256, 256, 256}; /* max 7 MC channels + 5 compatible channels */
  800.     double *sample;
  801.     DSBL *spike;
  802.     static D1408 *fft_buf;
  803.     static mask_ptr power;
  804.     static g_ptr ltg;
  805.     int j, l, z, q;
  806.     static int crit_band;
  807.     static int *cbound;
  808.     static int sub_size;
  809.    
  810.     sample = (double *) mem_alloc (sizeof (DFFT), "sample");
  811.     spike = (DSBL *) mem_alloc (sizeof (D12SBL), "spike");
  812.    
  813.     if (!init)
  814.     {
  815. /* bands, bark values, and mapping */
  816. /* changed 5 to 7 for matricing 8/10/93,SR */
  817. /* changed 7 to 12 for aug matricing 8/11/96,FdB */
  818. fft_buf = (D1408 *) mem_alloc ((long) sizeof (D1408) * 12, "fft_buf");
  819.      power = (mask_ptr) mem_alloc (sizeof (mask) * HAN_SIZE, "power");
  820. /* call functions for critical boundaries, freq. */
  821. crit_band = read_crit_band (info->lay, info->sampling_frequency);
  822. cbound = (int *) mem_alloc (sizeof (int) * crit_band, "cbound");
  823. read_cbound (info->lay, info->sampling_frequency, crit_band, cbound);
  824. read_freq_band (&sub_size, &ltg, info->lay, info->sampling_frequency);
  825. make_map (sub_size, power, ltg);
  826. for (i = 0; i < 1408; i++)
  827.     fft_buf[0][i] = fft_buf[1][i] = fft_buf[2][i] =
  828.     fft_buf[3][i] = fft_buf[4][i] = fft_buf[5][i] =
  829.     fft_buf[6][i] = fft_buf[7][i] = fft_buf[8][i] =
  830.     fft_buf[9][i] = fft_buf[10][i] = fft_buf[11][i] = 0;
  831. init = 1;
  832.     }
  833.  
  834.     if (aiff != 1)
  835.     {
  836. j = 0;
  837. l = 2;
  838.     }
  839.     else
  840.     {    
  841. j = 0;
  842. if (stereoaug == 2) l = 12; /*fr_ps->stereo + fr_ps->stereomc + fr_ps->stereoaug + 5 compatible */
  843. else                l =  7; /*fr_ps->stereo + fr_ps->stereomc + 2 compatible */
  844.     } 
  846.     for (k = j; k < l; k++)
  847.     {
  848.         for (i = 0; i < 1152; i++)
  849.     fft_buf[k][(i + off[k]) % 1408] = (double) buffer[k][i] / SCALE;
  850. for (i = 0; i < FFT_SIZE; i++)
  851.     sample[i] = fft_buf[k][(i + 1216 + off[k]) % 1408];
  852.    
  853. off[k] += 1152;
  854. off[k] %= 1408;
  855. /* call functions for windowing PCM samples,*/
  856. II_hann_win (sample); /* location of spectral components in each  */
  857. for (i = 0; i < HAN_SIZE; i++)
  858.     power[i].x = DBMIN; /* subband with labeling */
  859. II_f_f_t (sample, power); /* locate remaining non- */
  860.        
  861. if (fr_ps->header->center == 3 && k == 2)
  862. {
  863.     /* set center to 0, 9/2/93,SR*/
  864.     /* add to Left and Right ? WtK */
  865.     for (z = 184; z < HAN_SIZE; z++)
  866. power[z].x = -103.670;    /* DBMIN + 96.330; */
  867. }
  868. II_pick_max (power, spike[k]); /* tonal sinusoidals,   */
  870. #ifdef PRINTOUT
  871. if (verbosity >= 3)
  872. {
  873.     printf ("nChannel %d", k);
  874.     printf ("nSignal value per subband, from the FFT:n");
  875.     for (i = 0; i < sblimit; i++)
  876. printf ("%5.1f dB  ", spike[k][i]);
  877.     printf ("nMax. signal peak per subband, SCF SPL:n");
  878.     for (i = 0; i < sblimit; i++)   /* from [II_smr] determine the SCF SPL */
  879. printf ("%5.1f dB  ", 20 * log10(scale[k][i] * 32768)); 
  880.     fflush (stdout);
  881. }
  882. #endif
  884. II_tonal_label (power, &tone);                /* reduce noise & tonal components , find */
  885. noise_label (crit_band, cbound, power, &noise, ltg); 
  887. #ifdef PRINTOUT
  888. if (verbosity >= 3)
  889. {
  890.     printf ("nMaskers before sorting, FFT based levels:n");
  891.     for (i = 0; i < 511; i++)
  892.     {
  893. if ((power[i].type == NOISE) && (power[i].x > -200))
  894.     printf ("N:%3u %5.1f dB  ", i, power[i].x);
  895. if ((power[i].type == TONE) && (power[i].x > -200))
  896.     printf ("T:%3u %5.1f dB  ", i, power[i].x);
  897.     }
  898.     printf ("tone = %d noise = %d n", tone, noise);
  899.     fflush (stdout);
  900. }
  901. #endif
  903. subsampling (power, ltg, &tone, &noise);      /* global & minimal     */
  904.     
  905. #ifdef PRINTOUT
  906. if (verbosity >= 3)
  907. {
  908.     printf ("nMaskers after sorting:n");
  909.     for (i = 0; i < 511; i++)
  910.     {
  911. if ((power[i].type == NOISE) && (power[i].x > -200))
  912.     printf ("N:%3u %5.1f dB  ", i, power[i].x);
  913. if ((power[i].type == TONE) && (power[i].x > -200))
  914.     printf ("T:%3u %5.1f dB  ", i, power[i].x);
  915.     }
  916.     fflush (stdout);
  917. }
  918. #endif
  920. threshold (sub_size, power, ltg, &tone, &noise,
  921.    bitrate[info->lay-1][info->bitrate_index] / (stereo+stereomc+stereoaug)); /*to-mask ratio*//* 21/03/1995 JMZ BUG ???!!!???*/
  922. II_minimum_mask (sub_size, ltg, &ltmin[k][0], sblimit);
  923.     
  924. #ifdef PRINTOUT
  925. if (verbosity >= 3)
  926. {
  927.     printf ("nMinimum masking threshold:n");
  928.     for (i = 0; i < sblimit; i++)
  929. printf ("%5.1f dB  ", ltmin[k][i]);
  930.     fflush (stdout);
  931. }
  932. #endif
  934. for (i = 0; i < SBLIMIT; i++)
  935.     spiki[k][i] = spike[k][i];
  936.      
  937. i = 0;
  938. q = sblimit;
  939. II_smr (ltmin[k], smr[k], spike[k], scale[k], sblimit, i, q);
  940.     }
  941.     
  942.     mem_free ((void **) &sample);
  943.     mem_free ((void **) &spike); 
  944. } /*II_Psycho_One*/
  946. void II_Psycho_One_ml (
  947. double (*buffer)[1152],
  948. double (*scale)[32],
  949. double (*ltmin)[32],
  950. frame_params *fr_ps,
  951. double (*smr)[32],
  952. double (*spiki)[32]
  953. )
  954. {
  955.     layer *info = fr_ps->header;
  956.     int n_ml_ch   = info->multiling_ch;
  957.     int sblimit_ml = fr_ps->sblimit_ml;
  958.     int k,i, tone=0, noise=0;
  959.     static char init = 0;
  960.     static int off[7] = {256, 256, 256, 256, 256, 256, 256};   /* max 7 ML channels */
  961.     double *sample;
  962.     DSBL *spike;
  963.     static D1408 *fft_buf;
  964.     static mask_ptr power;
  965.     static g_ptr ltg_ml;
  966.     static int crit_band_ml;
  967.     static int *cbound_ml;
  968.     static int sub_size_ml;
  969.     int j, l, z, q;
  972.     sample = (double *) mem_alloc (sizeof (DFFT), "sample");
  973.     spike = (DSBL *) mem_alloc (sizeof (D7SBL), "spike");
  974.    
  975.     if (!init)
  976.     {
  977. /* bands, bark values, and mapping */
  978. fft_buf = (D1408 *) mem_alloc ((long) sizeof (D1408) * 7, "fft_buf");
  979. power = (mask_ptr) mem_alloc (sizeof (mask) * HAN_SIZE, "power");
  980. /* call functions for critical boundaries, freq. */
  981. if (info->multiling_fs == 0)
  982. {
  983.     crit_band_ml = read_crit_band (info->lay, info->sampling_frequency);
  984.     cbound_ml = (int *) mem_alloc (sizeof (int) * crit_band_ml, "cbound_ml");
  985.     read_cbound (info->lay, info->sampling_frequency, crit_band_ml, cbound_ml);
  986.     read_freq_band (&sub_size_ml, &ltg_ml, info->lay, info->sampling_frequency);
  987.     /* values are equal to those for the mc audio data */
  988. }
  989. else
  990. {
  991.     crit_band_ml = read_crit_band (info->lay, info->sampling_frequency + 4);
  992.     cbound_ml = (int *) mem_alloc (sizeof (int) * crit_band_ml, "cbound_ml");
  993.     read_cbound (info->lay, info->sampling_frequency + 4, crit_band_ml, cbound_ml);
  994.     read_freq_band (&sub_size_ml, &ltg_ml, info->lay, info->sampling_frequency + 4);
  995. }
  996. make_map (sub_size_ml, power, ltg_ml);
  997. for (i = 0; i < 1408; i++)
  998.     fft_buf[0][i] = fft_buf[1][i] = fft_buf[2][i] =     
  999.     fft_buf[3][i] = fft_buf[4][i] = fft_buf[5][i] = fft_buf[6][i] = 0;
  1000. init = 1;
  1001.     }
  1002.  
  1003.     if (n_ml_ch > 0) 
  1004.     {
  1005. for (k = 0; k < n_ml_ch; k++)
  1006. { 
  1007.     for (i = 0; i < 1152; i++)
  1008. fft_buf[k][(i+off[k]) % 1408]= (double) buffer[7+k][i] / SCALE;
  1009.     for (i = 0; i < FFT_SIZE; i++)
  1010. sample[i] = fft_buf[k][(i+1216+off[k]) % 1408];
  1011.     
  1012.     off[k] += 1152;
  1013.     off[k] %= 1408;
  1015.     /* call functions for windowing PCM samples,*/
  1016.     II_hann_win (sample); /* location of spectral components in each  */
  1017.     for (i = 0; i < HAN_SIZE; i++)
  1018. power[i].x = DBMIN; /* subband with labeling */
  1019.     II_f_f_t (sample, power); /* locate remaining non- */
  1020.     
  1021.     II_pick_max (power, spike[k]);
  1022. /* tonal sinusoidals,   */
  1024.     II_tonal_label (power, &tone); /* reduce noise & tonal components , find */
  1025.     noise_label (crit_band_ml, cbound_ml, power, &noise, ltg_ml); 
  1027.     subsampling (power, ltg_ml, &tone, &noise); /* global & minimal     */
  1028.     
  1029.     threshold (sub_size_ml, power, ltg_ml, &tone, &noise,
  1030.        bitrate[info->lay-1][info->bitrate_index] / 2); /* to-mask ratio */
  1031. /* threshold, and sgnl- */
  1032.     /* fprintf(stderr,  "sblimit_ml : %dn",  sblimit_ml); fflush(stderr); */
  1033.     II_minimum_mask (sub_size_ml, ltg_ml, &ltmin[7+k][0], sblimit_ml);
  1034.     
  1035.     for (i = 0; i < SBLIMIT; i++)
  1036. spiki[7+k][i] = spike[k][i];
  1037.  
  1038.     i = 0;
  1039.     q = sblimit_ml;
  1040.     II_smr (ltmin[7+k], smr[7+k], spike[k], scale[7+k], sblimit_ml, i, q);
  1042. } /* k-loop*/
  1043.     } /*n_ml_ch>0*/
  1045.     mem_free ((void **) &sample);
  1046.     mem_free ((void **) &spike); 
  1047. } /* II_psycho_One_ml */
  1049. /**********************************************************************/
  1050. /*
  1051. /*        This module implements the psychoacoustic model I for the
  1052. /* MPEG encoder layer I. It uses simplified tonal and noise masking
  1053. /* threshold analysis to generate SMR for the encoder bit allocation
  1054. /* routine.
  1055. /*
  1056. /**********************************************************************/
  1058. /****************************************************************/
  1059. /*
  1060. /*        Fast Fourier transform of the input samples.
  1061. /*
  1062. /****************************************************************/
  1064. void I_f_f_t (double *sample, mask *power)  /* this function calculates */
  1065.     /* an FFT analysis for the  */
  1066.     /* freq. domain             */
  1067. {
  1068.  int i,j,k,ll,l=0;
  1069.  int ip, le, le1;
  1070.  double t_r, t_i, u_r, u_i;
  1071.  static int M, MM1, init = 0, N, NV2, NM1;
  1072.  double *x_r, *x_i, *energy;
  1073.  static int *rev;
  1074.  static double *w_r, *w_i;
  1076.  x_r = (double *) mem_alloc(sizeof(DFFT2), "x_r");
  1077.  x_i = (double *) mem_alloc(sizeof(DFFT2), "x_i");
  1078.  energy = (double *) mem_alloc(sizeof(DFFT2), "energy");
  1079.  for(i=0;i<FFT_SIZE/2;i++) x_r[i] = x_i[i] = energy[i] = 0;
  1080.  if(!init){
  1081.     rev = (int *) mem_alloc(sizeof(IFFT2), "rev");
  1082.     w_r = (double *) mem_alloc(sizeof(D9), "w_r");
  1083.     w_i = (double *) mem_alloc(sizeof(D9), "w_i");
  1084.     M = 9;
  1085.     MM1 = 8;
  1086.     N = FFT_SIZE/2;
  1087.     NV2 = FFT_SIZE/2 >> 1;
  1088.     NM1 = FFT_SIZE/2 - 1;
  1089.     for(ll=0;ll<M;ll++){
  1090.        le = 1 << (M-ll);
  1091.        le1 = le >> 1;
  1092.        w_r[ll] = cos(PI/le1);
  1093.        w_i[ll] = -sin(PI/le1);
  1094.     }
  1095.     for(i=0;i<FFT_SIZE/2;rev[i] = l,i++) for(j=0,l=0;j<9;j++){
  1096.        k=(i>>j) & 1;
  1097.        l |= (k<<8-j);                
  1098.     }
  1099.     init = 1;
  1100.  }
  1101.  memcpy( (char *) x_r, (char *) sample, sizeof(double) * FFT_SIZE/2);
  1102.  for(ll=0;ll<MM1;ll++){
  1103.     le = 1 << (M-ll);
  1104.     le1 = le >> 1;
  1105.     u_r = 1;
  1106.     u_i = 0;
  1107.     for(j=0;j<le1;j++){
  1108.        for(i=j;i<N;i+=le){
  1109.           ip = i + le1;
  1110.           t_r = x_r[i] + x_r[ip];
  1111.           t_i = x_i[i] + x_i[ip];
  1112.           x_r[ip] = x_r[i] - x_r[ip];
  1113.           x_i[ip] = x_i[i] - x_i[ip];
  1114.           x_r[i] = t_r;
  1115.           x_i[i] = t_i;
  1116.           t_r = x_r[ip];
  1117.           x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
  1118.           x_i[ip] = x_i[ip] * u_r + t_r * u_i;
  1119.        }
  1120.        t_r = u_r;
  1121.        u_r = u_r * w_r[ll] - u_i * w_i[ll];
  1122.        u_i = u_i * w_r[ll] + t_r * w_i[ll];
  1123.     }
  1124.  }
  1125.  for(i=0;i<N;i+=2){
  1126.     ip = i + 1;
  1127.     t_r = x_r[i] + x_r[ip];
  1128.     t_i = x_i[i] + x_i[ip];
  1129.     x_r[ip] = x_r[i] - x_r[ip];
  1130.     x_i[ip] = x_i[i] - x_i[ip];
  1131.     x_r[i] = t_r;
  1132.     x_i[i] = t_i;
  1133.     energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
  1134.  }
  1135.  for(i=0;i<FFT_SIZE/2;i++) if(i<rev[i]){
  1136.     t_r = energy[i];
  1137.     energy[i] = energy[rev[i]];
  1138.     energy[rev[i]] = t_r;
  1139.  }
  1140.  for(i=0;i<HAN_SIZE/2;i++){                     /* calculate power  */
  1141.     if(energy[i] < 1E-20) energy[i] = 1E-20;    /* density spectrum */
  1142.        /* power calculation corrected with a factor 4, both positive
  1143.   and negative frequencies exist, 1992-11-06 shn */
  1144.        power[i].x = 10 * log10(energy[i]*4) + POWERNORM;
  1145.        power[i].next = STOP;
  1146.        power[i].type = FALSE;
  1147.  }
  1148.     mem_free ((void **) &x_r);
  1149.     mem_free ((void **) &x_i);
  1150.     mem_free ((void **) &energy);
  1151. }
  1153. /****************************************************************/
  1154. /*
  1155. /*         Window the incoming audio signal.
  1156. /*
  1157. /****************************************************************/
  1159. void I_hann_win (double *sample)    /* this function calculates a  */
  1160.     /* Hann window for PCM (input) */
  1161. {     /* samples for a 512-pt. FFT   */
  1162.  register int i;
  1163.  register double sqrt_8_over_3;
  1164.  static int init = 0;
  1165.  static double *window;
  1167.  if(!init){  /* calculate window function for the Fourier transform */
  1168.     window = (double *) mem_alloc (sizeof (DFFT2), "window");
  1169.     sqrt_8_over_3 = pow(8.0/3.0, 0.5);
  1170.     for(i=0;i<FFT_SIZE/2;i++){
  1171.       /* Hann window formula */
  1172.       window[i]=sqrt_8_over_3*0.5*(1-cos(2.0*PI*i/(FFT_SIZE/2-1)))/(FFT_SIZE/2);
  1173.     }
  1174.     init = 1;
  1175.  }
  1176.  for(i=0;i<FFT_SIZE/2;i++) sample[i] *= window[i];
  1177. }
  1179. /*******************************************************************/
  1180. /*
  1181. /*        This function finds the maximum spectral component in each
  1182. /* subband and return them to the encoder for time-domain threshold
  1183. /* determination.
  1184. /*
  1185. /*******************************************************************/
  1187. void I_pick_max (mask *power, double *spike)
  1188. {
  1189.  double max;
  1190.  int i,j;
  1192.  /* calculate the spectral component in each subband */
  1193.  for(i=0;i<HAN_SIZE/2;spike[i>>3] = max, i+=8)
  1194.     for(j=0, max = DBMIN;j<8;j++) max = (max>power[i+j].x) ? max : power[i+j].x;
  1195. }
  1197. /****************************************************************/
  1198. /*
  1199. /*        This function labels the tonal component in the power
  1200. /* spectrum.
  1201. /*
  1202. /****************************************************************/
  1204. void I_tonal_label (mask *power, int *tone) /* this function extracts   */
  1205.     /* (tonal) sinusoidals from */
  1206.     /* the spectrum             */
  1207. {
  1208.  int i,j, last = LAST, first, run;
  1209.  double max;
  1210.  int last_but_one=LAST;
  1212.  *tone = LAST;
  1213.  for(i=2;i<HAN_SIZE/2-6;i++){
  1214.     if(power[i].x>power[i-1].x && power[i].x>=power[i+1].x){
  1215.        power[i].type = TONE;
  1216.        power[i].next = LAST;
  1217.        if(last != LAST) power[last].next = i;
  1218.        else first = *tone = i;
  1219.        last = i;
  1220.     }
  1221.  }
  1222.  last = LAST;
  1223.  first = *tone;
  1224.  *tone = LAST;
  1225.  while(first != LAST){                /* conditions for the tonal     */
  1226.     if(first<2 || first>249) run = 0; /* otherwise k+/-j will be out of bounds*/
  1227.     else if(first<62) run = 2;        /* components in layer II, which */
  1228.     else if(first<126) run = 3;       /* are the boundaries for calc.   */
  1229.     else run = 6;                     /* the tonal components          */
  1230.     max = power[first].x - 7;
  1231.     for(j=2;j<=run;j++)  /* after calc. of tonal components, set to loc.*/
  1232.        if(max < power[first-j].x || max < power[first+j].x){   /* max   */
  1233.           power[first].type = FALSE;
  1234.           break;
  1235.        }
  1236.     if(power[first].type == TONE){    /* extract tonal components */
  1237.        int help=first;
  1238.        if(*tone == LAST) *tone = first;
  1239.        while((power[help].next!=LAST)&&(power[help].next-first)<=run)
  1240.           help=power[help].next;
  1241.        help=power[help].next;
  1242.        power[first].next=help;
  1243.        if((first-last)<=run){
  1244.           if(last_but_one != LAST) power[last_but_one].next=first;
  1245.        }
  1246.        if(first>1 && first<255){     /* calculate the sum of the */
  1247.           double tmp;                /* powers of the components */
  1248.           tmp = add_db(power[first-1].x, power[first+1].x);
  1249.           power[first].x = add_db(power[first].x, tmp);
  1250.        }
  1251.        for(j=1;j<=run;j++){
  1252.           power[first-j].x = power[first+j].x = DBMIN;
  1253.           power[first-j].next = power[first+j].next = STOP; /*dpwe: 2nd was .x*/
  1254.           power[first-j].type = power[first+j].type = FALSE;
  1255.        }
  1256.        last_but_one=last;
  1257.        last = first;
  1258.        first = power[first].next;
  1259.     }
  1260.     else {
  1261.        int ll;
  1262.        if(last == LAST) ; /* *tone = power[first].next; dpwe */
  1263.        else power[last].next = power[first].next;
  1264.        ll = first;
  1265.        first = power[first].next;
  1266.        power[ll].next = STOP;
  1267.     }
  1268.  }
  1269. }                        
  1270.                                 
  1271. /****************************************************************/
  1272. /*
  1273. /*        This function finds the minimum masking threshold and
  1274. /* return the value to the encoder.
  1275. /*
  1276. /****************************************************************/
  1278. void I_minimum_mask(int sub_size, g_thres *ltg, double *ltmin)
  1279. {
  1280.  double min;
  1281.  int i,j;
  1283.  j=1;
  1284.  for(i=0;i<SBLIMIT;i++)
  1285.     if(j>=sub_size-1)                   /* check subband limit, and       */
  1286.        ltmin[i] = ltg[sub_size-1].hear; /* calculate the minimum masking  */
  1287.     else {                              /* level of LTMIN for each subband*/
  1288.        min = ltg[j].x;
  1289.        while(ltg[j].line>>3 == i && j < sub_size){
  1290.           if (min>ltg[j].x)  min = ltg[j].x;
  1291.           j++;
  1292.        }
  1293.        ltmin[i] = min;
  1294.     }
  1295. }
  1297. /*****************************************************************/
  1298. /*
  1299. /*        This procedure is called in musicin to pick out the
  1300. /* smaller of the scalefactor or threshold.
  1301. /*
  1302. /*****************************************************************/
  1304. void I_smr (double *ltmin, double *spike, double *scale)
  1305. {
  1306.  int i,j;
  1307.  double max;
  1308.                 
  1309.  for(i=0;i<SBLIMIT;i++){                      /* determine the signal   */
  1310.     max = 20 * log10(scale[i] * 32768) - 10;  /* level for each subband */
  1311.     if(spike[i]>max) max = spike[i];          /* for the scalefactor    */
  1312.     max -= ltmin[i];
  1313.     ltmin[i] = max;
  1314.  }
  1315. }
  1316.         
  1317. /****************************************************************/
  1318. /*
  1319. /*        This procedure calls all the necessary functions to
  1320. /* complete the psychoacoustic analysis.
  1321. /*
  1322. /****************************************************************/
  1324. void
  1325. I_Psycho_One(
  1326. double (*buffer)[1152],
  1327. double (*scale)[32],
  1328. double (*ltmin)[32],
  1329. frame_params *fr_ps
  1330. ) {
  1331.  int stereo = fr_ps->stereo;
  1332.  the_layer info = fr_ps->header;
  1333.  int k,i, tone=0, noise=0;
  1334.  static char init = 0;
  1335.  static int off[2] = {256,256};
  1336.  double *sample;
  1337.  DSBL *spike;
  1338.  static D640 *fft_buf;
  1339.  static mask_ptr power;
  1340.  static g_ptr ltg;
  1341.  static int crit_band;
  1342.  static int *cbound;
  1343.  static int sub_size;
  1345.  sample = (double *) mem_alloc(sizeof(DFFT2), "sample");
  1346.  spike = (DSBL *) mem_alloc(sizeof(D2SBL), "spike");
  1347.             /* call functions for critical boundaries, freq. */
  1348.  if(!init){ /* bands, bark values, and mapping              */
  1349.     fft_buf = (D640 *) mem_alloc(sizeof(D640) * 2, "fft_buf");
  1350.     power = (mask_ptr) mem_alloc(sizeof(mask) * HAN_SIZE/2, "power");
  1351.     crit_band = read_crit_band(info->lay,info->sampling_frequency);
  1352.     cbound = (int *) mem_alloc(sizeof(int) * crit_band, "cbound");
  1353.     read_cbound(info->lay,info->sampling_frequency,crit_band,cbound);
  1354.     read_freq_band(&sub_size,&ltg,info->lay,info->sampling_frequency);
  1355.     make_map(sub_size,power,ltg);
  1356.     for(i=0;i<640;i++) fft_buf[0][i] = fft_buf[1][i] = 0;
  1357.     init = 1;
  1358.  }
  1359.  for(k=0;k<stereo;k++){    /* check PCM input for a block of */
  1360.     for(i=0;i<384;i++)     /* 384 samples for a 512-pt. FFT  */
  1361.        fft_buf[k][(i+off[k])%640]= (double) buffer[k][i]/SCALE;
  1362.     for(i=0;i<FFT_SIZE/2;i++)
  1363.        sample[i] = fft_buf[k][(i+448+off[k])%640];
  1364.     off[k] += 384;
  1365.     off[k] %= 640;
  1366.                         /* call functions for windowing PCM samples,   */
  1367.     I_hann_win(sample); /* location of spectral components in each     */
  1368.     for(i=0;i<HAN_SIZE/2;i++) power[i].x = DBMIN;   /* subband with    */
  1369.     I_f_f_t(sample, power);              /* labeling, locate remaining */
  1370.     I_pick_max(power, &spike[k][0]);     /* non-tonal sinusoidals,     */
  1371.     I_tonal_label(power, &tone);         /* reduce noise & tonal com., */
  1372.     noise_label(crit_band,cbound,power, &noise, ltg);     /* find global & minimal      */
  1373.     subsampling (power, ltg, &tone, &noise);  /* threshold, and sgnl-   */
  1374.     threshold(sub_size,power, ltg, &tone, &noise,     /* to-mask ratio          */
  1375.       bitrate[info->lay-1][info->bitrate_index]/stereo);
  1376.     I_minimum_mask(sub_size,ltg, &ltmin[k][0]);
  1377.     I_smr(&ltmin[k][0], &spike[k][0], &scale[k][0]);        
  1378.  }
  1379.     mem_free ((void **) &sample);
  1380.     mem_free ((void **) &spike);
  1381. }