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