home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 6 / AACD06.ISO / AACD / Sound / LAME / Source / psy.c < prev    next >
C/C++ Source or Header  |  1999-05-31  |  10KB  |  295 lines
  1.  
  2. /**********************************************************************
  3.  * ISO MPEG Audio Subgroup Software Simulation Group (1996)
  4.  * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension
  5.  *
  6.  * $Id: psy.c,v 1.2 1998/10/05 17:06:48 larsi Exp $
  7.  *
  8.  * $Log: psy.c,v $
  9.  * Revision 1.2  1998/10/05 17:06:48  larsi
  10.  * *** empty log message ***
  11.  *
  12.  * Revision 1.1.1.1  1998/10/05 14:47:18  larsi
  13.  *
  14.  * Revision 1.1  1996/02/14 04:04:23  rowlands
  15.  * Initial revision
  16.  *
  17.  * Received from Mike Coleman
  18.  **********************************************************************/
  19. #include "common.h"
  20. #include "encoder.h"
  21.  
  22. #if 0
  23. FILE        *fpo;    /* file pointer */
  24. void psycho_anal(buffer,savebuf,chn,lay,snr32,sfreq)
  25. short int *buffer;
  26. short int savebuf[1056];
  27. int   chn, lay;
  28. FLOAT snr32[32];
  29. double sfreq;        /* to match prototype : float args are always double */
  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.      fpo = fopen("out.dat", "wb");
  106.     if(fpo==NULL) {
  107.         puts("\t The attempt to open the output file failed.\n");
  108.         exit(-1);}
  109.      cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");
  110.      rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");
  111.      window = (FLOAT *) mem_alloc(sizeof(FBLK), "window");
  112.      absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr");
  113.      tmn = (double *) mem_alloc(sizeof(DCB), "tmn");
  114.      s = (FCB *) mem_alloc(sizeof(FCBCB), "s");
  115.      lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");
  116.      r = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "r");
  117.      phi_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "phi_sav");
  118.  
  119.      i = sfreq + 0.5;
  120.      switch(i){
  121.         case 32000: sfreq_idx = 0; break;
  122.         case 44100: sfreq_idx = 1; break;
  123.         case 48000: sfreq_idx = 2; break;
  124.         default:    printf("error, invalid sampling frequency: %d Hz\n",i);
  125.         exit(-1);
  126.      }
  127.      printf("absthr[][] sampling frequency index: %d\n",sfreq_idx);
  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.   case 3:
  232.      printf("MFCX1 - should never get here - layer 3 is not currently supported\n");
  233.      break;
  234.   default:
  235.      printf("error, invalid MPEG/audio coding layer: %d\n",lay);
  236.  }
  237.  
  238. /* These mem_free() calls must correspond with the mem_alloc() calls     */
  239. /* used at the beginning of this function to simulate "automatic"        */
  240. /* variables placed on the stack.                                        */
  241.  
  242.  mem_free((void **) &grouped_c);
  243.  mem_free((void **) &grouped_e);
  244.  mem_free((void **) &nb);
  245.  mem_free((void **) &cb);
  246.  mem_free((void **) &ecb);
  247.  mem_free((void **) &bc);
  248.  mem_free((void **) &wsamp_r);
  249.  mem_free((void **) &wsamp_i);
  250.  mem_free((void **) &phi);
  251.  mem_free((void **) &energy);
  252.  mem_free((void **) &c);
  253.  mem_free((void **) &fthr);
  254.  mem_free((void **) &snrtmp);
  255. }
  256. #endif
  257.  
  258.  
  259. /******************************************************************************
  260. routine to read in absthr table from a file.
  261. ******************************************************************************/
  262.  
  263. extern double absthr_0[];
  264. extern double absthr_1[];
  265. extern double absthr_2[];
  266.  
  267. void read_absthr(absthr, table)
  268. FLOAT *absthr;
  269. int table;
  270. {
  271.   double *a= NULL;
  272.   int j;
  273.   
  274.   switch(table){
  275.   case 0:
  276.     a = absthr_0;
  277.     break;
  278.   case 1:
  279.     a = absthr_1;
  280.     break;
  281.   case 2:
  282.     a = absthr_2;
  283.     break;
  284.   default:
  285.     fprintf(stderr, "absthr table: Not valid table number\n");
  286.     exit( EXIT_FAILURE );
  287.   } 
  288.   for (j=0; j<HBLKSIZE; j++){
  289.     absthr[j] =  a[j];
  290.   }
  291.  
  292. }
  293.  
  294.  
  295.