home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / converters / audiompeg_1 / !AudioMPEG / c / psy < prev    next >
Text File  |  1995-05-14  |  19KB  |  443 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. psy.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.9                                                        *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  15.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  16.  * 7/10/91  Earle Jennings      Ported to MsDos.                      *
  17.  *                              replace of floats with FLOAT          *
  18.  * 2/11/92  W. Joseph Carter    Fixed mem_alloc() arg for "absthr".   *
  19.  * 7/24/92  M. Iwadare          HANN window coefficients modified.    *
  20.  * 7/27/92  Masahiro Iwadare    Bug fix, FFT modification for Layer 3 *
  21.  * 7/27/92  Masahiro Iwadare    Bug fix, "new", "old", and "oldest"   *
  22.  *                              updates                               *
  23.  * 8/07/92  Mike Coleman        Bug fix, read_absthr()                *
  24.  **********************************************************************/
  25.  
  26. #include "common.h"
  27. #include "encoder.h"
  28.  
  29. void psycho_anal(buffer,savebuf,chn,lay,snr32,sfreq)
  30. short int *buffer;
  31. short int savebuf[1056];
  32. int   chn, lay;
  33. FLOAT snr32[32];
  34. double sfreq;        /* to match prototype : float args are always double */
  35. {
  36.  unsigned int   i, j, k;
  37.  FLOAT          r_prime, phi_prime;
  38.  FLOAT          freq_mult, bval_lo, minthres, sum_energy;
  39.  double         tb, temp1, temp2, temp3;
  40.  
  41. /* The static variables "r", "phi_sav", "new", "old" and "oldest" have    */
  42. /* to be remembered for the unpredictability measure.  For "r" and        */
  43. /* "phi_sav", the first index from the left is the channel select and     */
  44. /* the second index is the "age" of the data.                             */
  45.  
  46.  static int     new = 0, old = 1, oldest = 0;
  47.  static int     init = 0, flush, sync_flush, syncsize, sfreq_idx;
  48.  
  49. /* The following static variables are constants.                           */
  50.  
  51.  static double  nmt = 5.5;
  52.  
  53.  static FLOAT   crit_band[27] = {0,  100,  200, 300, 400, 510, 630,  770,
  54.                                920, 1080, 1270,1480,1720,2000,2320, 2700,
  55.                               3150, 3700, 4400,5300,6400,7700,9500,12000,
  56.                              15500,25000,30000};
  57.  
  58.  static FLOAT   bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,
  59.                             10.0,  7.0,  4.4,  4.5,  4.5,  4.5,  4.5,
  60.                              4.5,  4.5,  4.5,  4.5,  4.5,  4.5,  4.5,
  61.                              4.5,  4.5,  4.5,  3.5,  3.5,  3.5};
  62.  
  63. /* The following pointer variables point to large areas of memory         */
  64. /* dynamically allocated by the mem_alloc() function.  Dynamic memory     */
  65. /* allocation is used in order to avoid stack frame or data area          */
  66. /* overflow errors that otherwise would have occurred at compile time     */
  67. /* on the Macintosh computer.                                             */
  68.  
  69.  FLOAT          *grouped_c, *grouped_e, *nb, *cb, *ecb, *bc;
  70.  FLOAT          *wsamp_r, *wsamp_i, *phi, *energy;
  71.  FLOAT          *c, *fthr;
  72.  F32            *snrtmp;
  73.  
  74.  static int     *numlines;
  75.  static int     *partition;
  76.  static FLOAT   *cbval, *rnorm;
  77.  static FLOAT   *window;
  78.  static FLOAT   *absthr;
  79.  static double  *tmn;
  80.  static FCB     *s;
  81.  static FHBLK   *lthr;
  82.  static F2HBLK  *r, *phi_sav;
  83.  
  84. /* These dynamic memory allocations simulate "automatic" variables        */
  85. /* placed on the stack.  For each mem_alloc() call here, there must be    */
  86. /* a corresponding mem_free() call at the end of this function.           */
  87.  
  88.  grouped_c = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_c");
  89.  grouped_e = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_e");
  90.  nb = (FLOAT *) mem_alloc(sizeof(FCB), "nb");
  91.  cb = (FLOAT *) mem_alloc(sizeof(FCB), "cb");
  92.  ecb = (FLOAT *) mem_alloc(sizeof(FCB), "ecb");
  93.  bc = (FLOAT *) mem_alloc(sizeof(FCB), "bc");
  94.  wsamp_r = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_r");
  95.  wsamp_i = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_i");
  96.  phi = (FLOAT *) mem_alloc(sizeof(FBLK), "phi");
  97.  energy = (FLOAT *) mem_alloc(sizeof(FBLK), "energy");
  98.  c = (FLOAT *) mem_alloc(sizeof(FHBLK), "c");
  99.  fthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "fthr");
  100.  snrtmp = (F32 *) mem_alloc(sizeof(F2_32), "snrtmp");
  101.  
  102.  if(init==0){
  103.  
  104. /* These dynamic memory allocations simulate "static" variables placed    */
  105. /* in the data space.  Each mem_alloc() call here occurs only once at     */
  106. /* initialization time.  The mem_free() function must not be called.      */
  107.  
  108.      numlines = (int *) mem_alloc(sizeof(ICB), "numlines");
  109.      partition = (int *) mem_alloc(sizeof(IHBLK), "partition");
  110.      cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");
  111.      rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");
  112.      window = (FLOAT *) mem_alloc(sizeof(FBLK), "window");
  113.      absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr");
  114.      tmn = (double *) mem_alloc(sizeof(DCB), "tmn");
  115.      s = (FCB *) mem_alloc(sizeof(FCBCB), "s");
  116.      lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");
  117.      r = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "r");
  118.      phi_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "phi_sav");
  119.  
  120.      i = sfreq + 0.5;
  121.      switch(i){
  122.         case 32000: sfreq_idx = 0; break;
  123.         case 44100: sfreq_idx = 1; break;
  124.         case 48000: sfreq_idx = 2; break;
  125.         default:    printf("error, invalid sampling frequency: %d Hz\n",i);
  126.         exit(-1);
  127.      }
  128.      read_absthr(absthr, sfreq_idx);
  129.      if(lay==1){
  130.         flush = 384;
  131.         syncsize = 1024;
  132.         sync_flush = 576;
  133.      }
  134.      else {
  135.         flush = 384*3.0/2.0;
  136.         syncsize = 1056;
  137.         sync_flush = syncsize - flush;
  138.      }
  139. /* calculate HANN window coefficients */
  140. /*   for(i=0;i<BLKSIZE;i++)window[i]=0.5*(1-cos(2.0*PI*i/(BLKSIZE-1.0))); */
  141.      for(i=0;i<BLKSIZE;i++)window[i]=0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE));
  142. /* reset states used in unpredictability measure */
  143.      for(i=0;i<HBLKSIZE;i++){
  144.         r[0][0][i]=r[1][0][i]=r[0][1][i]=r[1][1][i]=0;
  145.         phi_sav[0][0][i]=phi_sav[1][0][i]=0;
  146.         phi_sav[0][1][i]=phi_sav[1][1][i]=0;
  147.         lthr[0][i] = 60802371420160.0;
  148.         lthr[1][i] = 60802371420160.0;
  149.      }
  150. /*****************************************************************************
  151.  * Initialization: Compute the following constants for use later             *
  152.  *    partition[HBLKSIZE] = the partition number associated with each        *
  153.  *                          frequency line                                   *
  154.  *    cbval[CBANDS]       = the center (average) bark value of each          *
  155.  *                          partition                                        *
  156.  *    numlines[CBANDS]    = the number of frequency lines in each partition  *
  157.  *    tmn[CBANDS]         = tone masking noise                               *
  158.  *****************************************************************************/
  159. /* compute fft frequency multiplicand */
  160.      freq_mult = sfreq/BLKSIZE;
  161.  
  162. /* calculate fft frequency, then bval of each line (use fthr[] as tmp storage)*/
  163.      for(i=0;i<HBLKSIZE;i++){
  164.         temp1 = i*freq_mult;
  165.         j = 1;
  166.         while(temp1>crit_band[j])j++;
  167.         fthr[i]=j-1+(temp1-crit_band[j-1])/(crit_band[j]-crit_band[j-1]);
  168.      }
  169.      partition[0] = 0;
  170. /* temp2 is the counter of the number of frequency lines in each partition */
  171.      temp2 = 1;
  172.      cbval[0]=fthr[0];
  173.      bval_lo=fthr[0];
  174.      for(i=1;i<HBLKSIZE;i++){
  175.         if((fthr[i]-bval_lo)>0.33){
  176.            partition[i]=partition[i-1]+1;
  177.            cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  178.            cbval[partition[i]] = fthr[i];
  179.            bval_lo = fthr[i];
  180.            numlines[partition[i-1]] = temp2;
  181.            temp2 = 1;
  182.         }
  183.         else {
  184.            partition[i]=partition[i-1];
  185.            cbval[partition[i]] += fthr[i];
  186.            temp2++;
  187.         }
  188.      }
  189.      numlines[partition[i-1]] = temp2;
  190.      cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  191.  
  192. /************************************************************************
  193.  * Now compute the spreading function, s[j][i], the value of the spread-*
  194.  * ing function, centered at band j, for band i, store for later use    *
  195.  ************************************************************************/
  196.      for(j=0;j<CBANDS;j++){
  197.         for(i=0;i<CBANDS;i++){
  198.            temp1 = (cbval[i] - cbval[j])*1.05;
  199.            if(temp1>=0.5 && temp1<=2.5){
  200.               temp2 = temp1 - 0.5;
  201.               temp2 = 8.0 * (temp2*temp2 - 2.0 * temp2);
  202.            }
  203.            else temp2 = 0;
  204.            temp1 += 0.474;
  205.            temp3 = 15.811389+7.5*temp1-17.5*sqrt((double) (1.0+temp1*temp1));
  206.            if(temp3 <= -100) s[i][j] = 0;
  207.            else {
  208.               temp3 = (temp2 + temp3)*LN_TO_LOG10;
  209.               s[i][j] = exp(temp3);
  210.            }
  211.         }
  212.      }
  213.  
  214.   /* Calculate Tone Masking Noise values */
  215.      for(j=0;j<CBANDS;j++){
  216.         temp1 = 15.5 + cbval[j];
  217.         tmn[j] = (temp1>24.5) ? temp1 : 24.5;
  218.   /* Calculate normalization factors for the net spreading functions */
  219.         rnorm[j] = 0;
  220.         for(i=0;i<CBANDS;i++){
  221.            rnorm[j] += s[j][i];
  222.         }
  223.      }
  224.      init++;
  225.  }
  226.  
  227. /************************* End of Initialization *****************************/
  228.  switch(lay) {
  229.   case 1:
  230.   case 2:
  231.      for(i=0; i<lay; i++){
  232. /*****************************************************************************
  233.  * Net offset is 480 samples (1056-576) for layer 2; this is because one must*
  234.  * stagger input data by 256 samples to synchronize psychoacoustic model with*
  235.  * filter bank outputs, then stagger so that center of 1024 FFT window lines *
  236.  * up with center of 576 "new" audio samples.                                *
  237.  *                                                                           *
  238.  * For layer 1, the input data still needs to be staggered by 256 samples,   *
  239.  * then it must be staggered again so that the 384 "new" samples are centered*
  240.  * in the 1024 FFT window.  The net offset is then 576 and you need 448 "new"*
  241.  * samples for each iteration to keep the 384 samples of interest centered   *
  242.  *****************************************************************************/
  243.         for(j=0; j<syncsize; j++){
  244.            if(j<(sync_flush))savebuf[j] = savebuf[j+flush];
  245.            else savebuf[j] = *buffer++;
  246.            if(j<BLKSIZE){
  247. /**window data with HANN window***********************************************/
  248.               wsamp_r[j] = window[j]*((FLOAT) savebuf[j]);
  249.               wsamp_i[j] = 0;
  250.            }
  251.         }
  252. /**Compute FFT****************************************************************/
  253.         fft(wsamp_r,wsamp_i,energy,phi,1024);
  254. /*****************************************************************************
  255.  * calculate the unpredictability measure, given energy[f] and phi[f]        *
  256.  *****************************************************************************/
  257. /*only update data "age" pointers after you are done with both channels      */
  258. /*for layer 1 computations, for the layer 2 double computations, the pointers*/
  259. /*are reset automatically on the second pass                                 */
  260.          if(lay==2 || (lay==1 && chn==0) ){
  261.            if(new==0){new = 1; oldest = 1;}
  262.            else {new = 0; oldest = 0;}
  263.            if(old==0)old = 1; else old = 0;
  264.         }
  265.         for(j=0; j<HBLKSIZE; j++){
  266.            r_prime = 2.0 * r[chn][old][j] - r[chn][oldest][j];
  267.            phi_prime = 2.0 * phi_sav[chn][old][j] - phi_sav[chn][oldest][j];
  268.            r[chn][new][j] = sqrt((double) energy[j]);
  269.            phi_sav[chn][new][j] = phi[j];
  270. temp1=r[chn][new][j] * cos((double) phi[j]) - r_prime * cos((double) phi_prime);
  271. temp2=r[chn][new][j] * sin((double) phi[j]) - r_prime * sin((double) phi_prime);
  272.            temp3=r[chn][new][j] + fabs((double)r_prime);
  273.            if(temp3 != 0)c[j]=sqrt(temp1*temp1+temp2*temp2)/temp3;
  274.            else c[j] = 0;
  275.         }
  276. /*****************************************************************************
  277.  * Calculate the grouped, energy-weighted, unpredictability measure,         *
  278.  * grouped_c[], and the grouped energy. grouped_e[]                          *
  279.  *****************************************************************************/
  280.         for(j=1;j<CBANDS;j++){
  281.            grouped_e[j] = 0;
  282.            grouped_c[j] = 0;
  283.         }
  284.         grouped_e[0] = energy[0];
  285.         grouped_c[0] = energy[0]*c[0];
  286.         for(j=1;j<HBLKSIZE;j++){
  287.            grouped_e[partition[j]] += energy[j];
  288.            grouped_c[partition[j]] += energy[j]*c[j];
  289.         }
  290. /*****************************************************************************
  291.  * convolve the grouped energy-weighted unpredictability measure             *
  292.  * and the grouped energy with the spreading function, s[j][k]               *
  293.  *****************************************************************************/
  294.         for(j=0;j<CBANDS;j++){
  295.            ecb[j] = 0;
  296.            cb[j] = 0;
  297.            for(k=0;k<CBANDS;k++){
  298.               if(s[j][k] != 0.0){
  299.                  ecb[j] += s[j][k]*grouped_e[k];
  300.                  cb[j] += s[j][k]*grouped_c[k];
  301.               }
  302.            }
  303.            if(ecb[j] !=0)cb[j] = cb[j]/ecb[j];
  304.            else cb[j] = 0;
  305.         }
  306. /*****************************************************************************
  307.  * Calculate the required SNR for each of the frequency partitions           *
  308.  *         this whole section can be accomplished by a table lookup          *
  309.  *****************************************************************************/
  310.         for(j=0;j<CBANDS;j++){
  311.            if(cb[j]<.05)cb[j]=0.05;
  312.            else if(cb[j]>.5)cb[j]=0.5;
  313.            tb = -0.434294482*log((double) cb[j])-0.301029996;
  314.            bc[j] = tmn[j]*tb + nmt*(1.0-tb);
  315.            k = cbval[j] + 0.5;
  316.            bc[j] = (bc[j] > bmax[k]) ? bc[j] : bmax[k];
  317.            bc[j] = exp((double) -bc[j]*LN_TO_LOG10);
  318.         }
  319. /*****************************************************************************
  320.  * Calculate the permissible noise energy level in each of the frequency     *
  321.  * partitions. Include absolute threshold and pre-echo controls              *
  322.  *         this whole section can be accomplished by a table lookup          *
  323.  *****************************************************************************/
  324.         for(j=0;j<CBANDS;j++)
  325.            if(rnorm[j] && numlines[j])
  326.               nb[j] = ecb[j]*bc[j]/(rnorm[j]*numlines[j]);
  327.            else nb[j] = 0;
  328.         for(j=0;j<HBLKSIZE;j++){
  329. /*temp1 is the preliminary threshold */
  330.            temp1=nb[partition[j]];
  331.            temp1=(temp1>absthr[j])?temp1:absthr[j];
  332. /*do not use pre-echo control for layer 2 because it may do bad things to the*/
  333. /*  MUSICAM bit allocation algorithm                                         */
  334.            if(lay==1){
  335.               fthr[j] = (temp1 < lthr[chn][j]) ? temp1 : lthr[chn][j];
  336.               temp2 = temp1 * 0.00316;
  337.               fthr[j] = (temp2 > fthr[j]) ? temp2 : fthr[j];
  338.            }
  339.            else fthr[j] = temp1;
  340.            lthr[chn][j] = LXMIN*temp1;
  341.         }
  342. /*****************************************************************************
  343.  * Translate the 512 threshold values to the 32 filter bands of the coder    *
  344.  *****************************************************************************/
  345.         for(j=0;j<193;j += 16){
  346.            minthres = 60802371420160.0;
  347.            sum_energy = 0.0;
  348.            for(k=0;k<17;k++){
  349.               if(minthres>fthr[j+k])minthres = fthr[j+k];
  350.               sum_energy += energy[j+k];
  351.            }
  352.            snrtmp[i][j/16] = sum_energy/(minthres * 17.0);
  353.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  354.         }
  355.         for(j=208;j<(HBLKSIZE-1);j += 16){
  356.            minthres = 0.0;
  357.            sum_energy = 0.0;
  358.            for(k=0;k<17;k++){
  359.               minthres += fthr[j+k];
  360.               sum_energy += energy[j+k];
  361.            }
  362.            snrtmp[i][j/16] = sum_energy/minthres;
  363.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  364.         }
  365. /*****************************************************************************
  366.  * End of Psychoacuostic calculation loop                                    *
  367.  *****************************************************************************/
  368.      }
  369.      for(i=0; i<32; i++){
  370.         if(lay==2)
  371.            snr32[i]=(snrtmp[0][i]>snrtmp[1][i])?snrtmp[0][i]:snrtmp[1][i];
  372.         else snr32[i]=snrtmp[0][i];
  373.      }
  374.      break;
  375.   case 3:
  376.      printf("layer 3 is not currently supported\n");
  377.      break;
  378.   default:
  379.      printf("error, invalid MPEG/audio coding layer: %d\n",lay);
  380.  }
  381.  
  382. /* These mem_free() calls must correspond with the mem_alloc() calls     */
  383. /* used at the beginning of this function to simulate "automatic"        */
  384. /* variables placed on the stack.                                        */
  385.  
  386.  mem_free((void **) &grouped_c);
  387.  mem_free((void **) &grouped_e);
  388.  mem_free((void **) &nb);
  389.  mem_free((void **) &cb);
  390.  mem_free((void **) &ecb);
  391.  mem_free((void **) &bc);
  392.  mem_free((void **) &wsamp_r);
  393.  mem_free((void **) &wsamp_i);
  394.  mem_free((void **) &phi);
  395.  mem_free((void **) &energy);
  396.  mem_free((void **) &c);
  397.  mem_free((void **) &fthr);
  398.  mem_free((void **) &snrtmp);
  399. }
  400.  
  401. /******************************************************************************
  402. routine to read in absthr table from a file.
  403. ******************************************************************************/
  404.  
  405. void read_absthr(absthr, table)
  406. FLOAT *absthr;
  407. int table;
  408. {
  409.  FILE *fp;
  410.  long j,index;
  411.  float a;
  412.  char t[80];
  413.  char ta[16];
  414.  
  415.  strcpy( ta, "absthr_0" );
  416.  
  417.  switch(table){
  418.     case 0 : ta[7] = '0';
  419.              break;
  420.     case 1 : ta[7] = '1';
  421.              break;
  422.     case 2 : ta[7] = '2';
  423.              break;
  424.     default : printf("absthr table: Not valid table number\n");
  425.  }
  426.  if(!(fp = OpenTableFile(ta) ) ){
  427.     printf("Please check %s table\n", ta);
  428.     exit(1);
  429.  }
  430.  fgets(t, 150, fp);
  431.  sscanf(t, "table %ld", &index);
  432.  if(index != table){
  433.     printf("error in absthr table %s",ta);
  434.     exit(1);
  435.  }
  436.  for(j=0; j<HBLKSIZE; j++){
  437.     fgets(t,80,fp);
  438.     sscanf(t,"%f", &a);
  439.     absthr[j] =  a;
  440.  }
  441.  fclose(fp);
  442. }
  443.