predisto.c

ISO mp3 sources (distribution 10) Layer 1/2/3, C Source, 512 k Sources of the Mpeg 1,2 layer 1,2

			m = 0; n = 1;
			scale_factor_calc(scalar,sblimit, m, n, i, sb_sample);
			trans_pattern(scalar, scfsi, fr_ps,scfsi_dyn);
			pre_quant(sbs_sample, m, n, sb_sample, scalar, bit_alloc, subband, fr_ps);

	 	    for( j = 0; j < 3; ++j)
			for(l = 0; l < 12; l ++)
			    for( k = sb_groups[i-1] + 1; k <= sb_groups[i]; k++)
			    {
				sbs_sample[2][j][l][k] = sbs_sample[0][j][l][k] - sbs_sample[5][j][l][k] - sbs_sample[3][j][l][k];
				sb_sample[1][j][l][k] = sbs_sample[6][j][l][k] + sbs_sample[2][j][l][k] + sb_sample[4][j][l][k];
			    }
			m = 1; n = 2;
			scale_factor_calc(scalar,sblimit, m, n, i, sb_sample);
			
		     break;
	    }
	  }
	}
 }




buffer_sample(unsigned int (*keep_it)[36][32], unsigned int (*sample)[3][32], unsigned int (*bit_alloc)[32], frame_params *fr_ps, int l, int z, int q)
{
   int i,j,k,m;
   int stereo = fr_ps->stereo;
   int sblimit = fr_ps->sblimit;
   int jsbound = fr_ps->jsbound;
   al_table *alloc = fr_ps->alloc;
   unsigned int x;

   for (i=0;i<sblimit;i++) for (j= l;j<((i<jsbound)?z:l+1);j++) 
   {
       if (bit_alloc[j][i]) 
       {
         /* check for grouping in subband */
	 /* if ((*alloc)[i][bit_alloc[j][i]].group==3)*/
           for (m=0;m<3;m++) 
	   {
              sample[j][m][i] = keep_it[j][q*3+m][i];
           }         
       }
       else 
       {      /* for no sample transmitted */
         for (k=0;k<3;k++) sample[j][k][i] = 0;
       }
       if(stereo == 2 && i>= jsbound)  /* joint stereo : copy L to R */
          for (k=0;k<3;k++) sample[1][k][i] = sample[0][k][i];
   }
   for (i=sblimit;i<SBLIMIT;i++) for (j= l;j< z;j++) for (k=0;k<3;k++)
      sample[j][k][i] = 0;
}      



static double a[17] = {
  0.750000000, 0.625000000, 0.875000000, 0.562500000, 0.937500000,
  0.968750000, 0.984375000, 0.992187500, 0.996093750, 0.998046875,
  0.999023438, 0.999511719, 0.999755859, 0.999877930, 0.999938965,
  0.999969482, 0.999984741 };
 
static double b[17] = {
  -0.250000000, -0.375000000, -0.125000000, -0.437500000, -0.062500000,
  -0.031250000, -0.015625000, -0.007812500, -0.003906250, -0.001953125,
  -0.000976563, -0.000488281, -0.000244141, -0.000122070, -0.000061035,
  -0.000030518, -0.000015259 };
 


void subband_quantization_pre(unsigned int (*scalar)[3][32], double (*sb_samples)[3][12][32], unsigned int (*bit_alloc)[32], unsigned int (*sbband)[3][12][32], frame_params *fr_ps, int ch1, int ch2)
{
   int i, j, k, s, n, qnt, sig, m, l, ll, z;
   int stereo	= fr_ps->stereo;
   int stereomc = fr_ps->stereomc;
   int sblimit = fr_ps->sblimit;
   unsigned int stps, y;
   double d;
   al_table *alloc = fr_ps->alloc;

 

   for (s=0;s<3;s++)
     for (j=0;j<12;j= j+3)
	for (i=0;i<sblimit;i++)
	  for(k = ch1; k < ch2; k++)
	  {
	    for(z = 0; z < 3; z ++)
	    {
		    if(bit_alloc[k][i])
	               d = sb_samples[k][s][j+z][i] / multiple[scalar[k][s][i]];
		    else
		       d = 0;
	   
	   
	    if (mod(d) >= 1.0) /* > changed to >=, 1992-11-06 shn */
	    {
		printf("In predistortion, not scaled properly, %d %d %d %d\n",k,s,j+z,i);
		printf("Value %1.10f\n",sb_samples[k][s][j+z][i]);
	    }
	    qnt = (*alloc)[i][bit_alloc[k][i]].quant;
            d = d * a[qnt] + b[qnt];
             /* extract MSB N-1 bits from the floating point sample */
             if (d >= 0) sig = 1;
             else { sig = 0; d += 1.0; }
             n = 0;
#ifndef MS_DOS
             stps = (*alloc)[i][bit_alloc[k][i]].steps;
             while ((1L<<n) < stps) n++;
#else
             while  ( ( (unsigned long)(1L<<(long)n) <
	                 ((unsigned long) ((*alloc)[i][bit_alloc[k][i]].steps)
                          & 0xffff
                         )
                      ) && ( n <16)
                    ) n++;
#endif
             n--;
             sbband[k][s][j+z][i] = (unsigned int) (d * (double) (1L<<n));
             /* tag the inverted sign bit to sbband at position N */
             /* The bit inversion is a must for grouping with 3,5,9 steps
                so it is done for all subbands */
             if (sig) sbband[k][s][j+z][i] |= 1<<n;
	     
	     if(!bit_alloc[k][i]) sbband[k][s][j+z][i] = 0;
	    }

/*try try try ! 7/7/94 */

/*             if ((*alloc)[i][bit_alloc[k][i]].group != 3) 
	     {
               y =(*alloc)[i][bit_alloc[k][i]].steps;
               sbband[k][s][j][i] = sbband[k][s][j][i] +
                      		    sbband[k][s][j+1][i] * y +
                                    sbband[k][s][j+2][i] * y * y;
               sbband[k][s][j+1][i] = sbband[k][s][j][i];
               sbband[k][s][j+2][i] = sbband[k][s][j][i];
	      }*/

	     
	     
	 }
           for (s=0;s<3;s++)
             for (j=sblimit;j<SBLIMIT;j++)
		for (i=0;i<12;i++) for (k = 0; k < 7; k++) sbband[k][s][i][j] = 0;
}






void trans_pattern(unsigned int (*scalar)[3][32], unsigned int (*scfsi)[32], frame_params *fr_ps, unsigned int (*scfsi_dyn)[32])
{
   int stereo  = fr_ps->stereo;
   int stereomc = fr_ps->stereomc;
   int sblimit = fr_ps->sblimit;
	int dscf[2];
	int class[2],i,j,k;
static int pattern[5][5] = {0x123, 0x122, 0x122, 0x133, 0x123,
                            0x113, 0x111, 0x111, 0x444, 0x113,
                            0x111, 0x111, 0x111, 0x333, 0x113,
                            0x222, 0x222, 0x222, 0x333, 0x123,
                            0x123, 0x122, 0x122, 0x133, 0x123};
 
  for (k = 0; k < 2; k++)
     for (i=0;i<sblimit;i++) 
     {
       dscf[0] =  (scalar[k][0][i]-scalar[k][1][i]);
       dscf[1] =  (scalar[k][1][i]-scalar[k][2][i]);
       for (j=0;j<2;j++)
       {
         if (dscf[j]<=-3) class[j] = 0;
         else if (dscf[j] > -3 && dscf[j] <0) class[j] = 1;
              else if (dscf[j] == 0) class[j] = 2;
                   else if (dscf[j] > 0 && dscf[j] < 3) class[j] = 3;
                        else class[j] = 4;
       }
       switch (pattern[class[0]][class[1]]) 
       {
         case 0x123 :    scfsi[k][i] = 0;
                         break;
         case 0x122 :    scfsi[k][i] = 3;
                         scalar[k][2][i] = scalar[k][1][i];
                         break;
         case 0x133 :    scfsi[k][i] = 3;
                         scalar[k][1][i] = scalar[k][2][i];
                         break;
         case 0x113 :    scfsi[k][i] = 1;
                         scalar[k][1][i] = scalar[k][0][i];
                         break;
         case 0x111 :    scfsi[k][i] = 2;
                         scalar[k][1][i] = scalar[k][2][i] = scalar[k][0][i];
                         break;
         case 0x222 :    scfsi[k][i] = 2;
                         scalar[k][0][i] = scalar[k][2][i] = scalar[k][1][i];
                         break;
         case 0x333 :    scfsi[k][i] = 2;
                         scalar[k][0][i] = scalar[k][1][i] = scalar[k][2][i];
                         break;
         case 0x444 :    scfsi[k][i] = 2;
                         if (scalar[k][0][i] > scalar[k][2][i])
                              scalar[k][0][i] = scalar[k][2][i];
                         scalar[k][1][i] = scalar[k][2][i] = scalar[k][0][i];
   
	}
    }
}    	   





void pre_quant(double (*sbs_sample)[3][12][32], int l, int m, double (*sb_sample)[3][12][32], 
		unsigned int (*scalar)[3][32], unsigned int (*bit_alloc)[32], 
		unsigned int (*subband)[3][12][32], frame_params *fr_ps)
{
    double sb_pre_sample[7][3][12][SBLIMIT];  /*predistortion*/
    unsigned int sample[7][3][SBLIMIT];
    unsigned int keep_it[7][36][SBLIMIT];
    int ch, sb, gr, i;
    int hi,  hu;

    subband_quantization_pre(scalar, sb_sample, bit_alloc, subband, fr_ps, l, m);

		
    for (ch = l; ch < m; ch++)
      for (sb = 0; sb < SBLIMIT; sb++)
	for (gr = 0; gr < 3; gr++)
	    for (i = 0; i < 12; i++)
		keep_it[ch][i+gr*12][sb] = subband[ch][gr][i][sb];

		for(i = 0; i < 12; i++) 
		{      
		    buffer_sample(keep_it,sample,bit_alloc,fr_ps, l, m, i);

		    II_dequantize_sample(sample, bit_alloc, sb_pre_sample, fr_ps, l, m, i);
		    II_denormalize_sample(sb_pre_sample, scalar, fr_ps, i>>2, l, m, i);
	        }      	    
		
    for (ch = l; ch < m; ch++)
      for (sb=0;sb<SBLIMIT;sb++)
	for (i = 0; i < 4; i++)
	    for (gr=0;gr<3;gr++)
		sbs_sample[ch][0][i*3 + gr][sb] = sb_pre_sample[ch][gr][i][sb];
		
    for (ch = l; ch < m; ch++)
      for (sb=0;sb<SBLIMIT;sb++)
	for (i = 4; i < 8; i++)
	    for (gr=0;gr<3;gr++)
		sbs_sample[ch][1][(i-4)*3 + gr][sb] = sb_pre_sample[ch][gr][i][sb];

    for (ch = l; ch < m; ch++)
      for (sb=0;sb<SBLIMIT;sb++)
	for (i = 8; i < 12; i++)
	    for (gr=0;gr<3;gr++)
		sbs_sample[ch][2][(i-8)*3 + gr][sb] = sb_pre_sample[ch][gr][i][sb];


}

/**************************************************************
/*
/*   Restore the compressed sample to a factional number.
/*   first complement the MSB of the sample
/*   for Layer II :
/*   Use the formula s = s' * c + d
/*
/*  taken out of decoder,modified for prediction 9/16/93,SR
/*
/**************************************************************/

static double c[17] = { 1.33333333333, 1.60000000000, 1.14285714286,
                        1.77777777777, 1.06666666666, 1.03225806452,
                        1.01587301587, 1.00787401575, 1.00392156863,
                        1.00195694716, 1.00097751711, 1.00048851979,
                        1.00024420024, 1.00012208522, 1.00006103888,
                        1.00003051851, 1.00001525902 };

static double d[17] = { 0.500000000, 0.500000000, 0.250000000, 0.500000000,
                        0.125000000, 0.062500000, 0.031250000, 0.015625000,
                        0.007812500, 0.003906250, 0.001953125, 0.0009765625,
                        0.00048828125, 0.00024414063, 0.00012207031,
                        0.00006103516, 0.00003051758 };

/************************** Layer II stuff ************************/

void II_dequantize_sample(unsigned int (*sample)[3][32], unsigned int (*bit_alloc)[32], double (*fraction)[3][12][32], frame_params *fr_ps, int l, int m, int z)
{
   int i, j, k, x;
   int stereo = fr_ps->stereo;
   int sblimit = fr_ps->sblimit;
   al_table *alloc = fr_ps->alloc;

   for (i=0;i<sblimit;i++)  for (j=0;j<3;j++) for (k = l;k < m ;k++)
     if (bit_alloc[k][i]) {

       /* locate MSB in the sample */
       x = 0;
#ifndef MSDOS
       while ((1L<<x) < (*alloc)[i][bit_alloc[k][i]].steps) x++;
#else
       /* microsoft C thinks an int is a short */
       while (( (unsigned long) (1L<<(long)x) <
                (unsigned long)( (*alloc)[i][bit_alloc[k][i]].steps)
              ) && ( x < 16) ) x++;
#endif

       /* MSB inversion */
       if (((sample[k][j][i] >> x-1) & 1) == 1)
          fraction[k][j][z][i] = 0.0;
       else  fraction[k][j][z][i] = -1.0;

       /* Form a 2's complement sample */
       fraction[k][j][z][i] += (double) (sample[k][j][i] & ((1<<x-1)-1)) /
                            (double) (1L<<x-1);

       /* Dequantize the sample */
       fraction[k][j][z][i] += d[(*alloc)[i][bit_alloc[k][i]].quant];
       fraction[k][j][z][i] *= c[(*alloc)[i][bit_alloc[k][i]].quant];
     }
     else fraction[k][j][z][i] = 0.0;   
   
   for (i=sblimit;i<SBLIMIT;i++) for (j=0;j<3;j++) for(k = l; k < m; k++)
     fraction[k][j][z][i] = 0.0;
}












void II_denormalize_sample(double (*fraction)[3][12][32], unsigned int (*scale_index)[3][32], frame_params *fr_ps, int x, int l, int m, int z)
       /*far*/                             
                                        
                    
      
            
{
   int i,j,k;
   int stereo = fr_ps->stereo;
   int sblimit = fr_ps->sblimit;

   for (i=0;i<sblimit;i++) for (j= l;j < m; j++) {
      fraction[j][0][z][i] *= multiple[scale_index[j][x][i]];
      fraction[j][1][z][i] *= multiple[scale_index[j][x][i]];
      fraction[j][2][z][i] *= multiple[scale_index[j][x][i]];
   }
}



void scale_factor_calc(unsigned int (*scalar)[3][32], int sblimit, int l, int m, int z, double (*sb_sample)[3][12][32])
                                   
            
            
                                          /* (see above) */
{
  int i,j, k,t;
  double s[SBLIMIT];
 
  for (k = l; k< m; k++) for (t=0;t<3;t++) 
  {
   for (i=((z == 0) ? 0:(sb_groups[z-1]+1));i<=sb_groups[z];i++) 
   {
	for (j=1, s[i] = mod(sb_sample[k][t][0][i]);j<SCALE_BLOCK;j++)
		  if (mod(sb_sample[k][t][j][i]) > s[i])
			 s[i] = mod(sb_sample[k][t][j][i]);
   } 
    for (i=((z == 0) ? 0:(sb_groups[z-1]+1));i<=sb_groups[z];i++)
    { 
        for (j=SCALE_RANGE-1,scalar[k][t][i]=0;j>=0;j--)
		if (s[i] < multiple[j]) 
		{                       /* <= changed to <, 1992-11-06 shn*/
                   scalar[k][t][i] = j;
                   break;
                }
     }
   }
}