ldelnr.c

语音CELP压缩解压源代码（C语音)

#define NFRAC 6
#define TRUE 1
#define FALSE 0
#define M1 -10
#define M2   9

/* five fractional delays plus zero delay = 6				*/

static float frac[NFRAC] = {0.25, 0.33333334, 0.5, 0.66666667, 0.75, 0.0};
static int twelfths[NFRAC] = {3, 4, 6, 8, 9, 0};

/**************************************************************************
*
* NAME
*		ldelay_nr
*
* FUNCTION
*		Time delay a bandlimited signal
*		using blockwise recursion (aka "nonrecursive")
*		with a long interpolation filter (see delay_ne.c)
*
* SYNOPSIS
*		subroutine ldelay_nr(x, start, n, d, m, y)
*
*   formal 
*                       data    I/O
*       name            type    type    function
*       -------------------------------------------------------------------
*	x[n]		float	i	signal input
*	start		int	i	beginning of output sequence
*	n		int	i	length of input sequence
*	d		float	i	fractional pitch
*	m		int	i	integer pitch
*	y[n]		float	o	delayed input signal
*
***************************************************************************
*
* DESCRIPTION
*
*	Subroutine to time delay a bandlimited signal by resampling the
*	reconstructed data (aka sinc interpolation).  The well known
*	reconstruction formula is:
*
*              |    M2      sin[pi(t-nT)/T]    M2
*   y(n) = X(t)| = SUM x(n) --------------- = SUM x(n) sinc[(t-nT)/T]
*              |   n=M1         pi(t-nT)/T    n=M1
*              |t=n+d
*
*	The interpolating (sinc) function is Hamming windowed to bandlimit
*	the signal to reduce aliasing.
*
*	Multiple simultaneous time delays may be efficiently calculated
*	by polyphase filtering.  Polyphase filters decimate the unused
*	filter coefficients.  See Chapter 6 in C&R for details. 
*
***************************************************************************
*
* CALLED BY
*
*	pitchvq	  psearch
*
* CALLS
*
*	ham
*
***************************************************************************
*
* REFERENCES
*
*	Crochiere & Rabiner, Multirate Digital Signal Processing,
*	P-H 1983, Chapters 2 and 6.
*
*	Kroon & Atal, "Pitch Predictors with High Temporal Resolution,"
*	ICASSP '90, S12.6
*
*	CELP Documentation, Version 3.1, Figure 9, 15 Dec '89.
*
**************************************************************************/

#define SIZE (M2 - M1 + 1)
ldelay_nr(x, start, n, d, m, y)
float x[], d, y[];
int m, n, start;
{
  static float wsinc[SIZE][NFRAC], hwin[12*SIZE+1];
  float pitch, pitch2, pitch3;
  static int first = TRUE;
  float sinc();
  int i, j, k, index, k2, k3, mtwo, qd();

  /* Generate Hamming windowed sinc interpolating function for each	*/
  /* allowable fraction.  The interpolating functions are stored in  	*/
  /* time reverse order (i.e., delay appears as advance) to align    	*/
  /* with the adaptive code book v0 array.  The interp filters are:  	*/
  /* 		wsinc[.,0]	frac = 1/4 (3/12)		     	*/
  /* 		wsinc[.,1]	frac = 1/3 (4/12)		     	*/
  /* 		.		.				     	*/
  /* 		wsinc[.,4]	frac = 3/4 (9/12)		     	*/
  /* 		wsinc[.,5]	frac = 0 zero delay for programming ease*/

  if (first)
  {
    ham(hwin, 12*SIZE+1);
    for (i = 0; i < NFRAC; i++)
      for (k = M1, j = 0; j < SIZE; j++)
      {
        wsinc[j][i] = sinc(frac[i] + k) * hwin[12*j+twelfths[i]];
        k++;
      }
    first = FALSE;
  }

  index = qd_nr(d);

  /* *Resample:							 	*/
  /* *If delay is greater than n(=60), then lay down n values at  	*/
  /* *the appropriate delay.				 		*/

  if (m >= n)
  {
    for (i = 0; i < n + M1; i++)
    {
      y[i] = 0.0;
      for (k = M1, j = 0; j < SIZE; j++)
      {
        y[i] += x[i+k+start-m-1] * wsinc[j][index];
        k++;
      }
    }
    for (i = 0; i < M2; i++)
      x[start+i-1] = y[i];
    for (i = n+M1; i < n; i++)
    {
      y[i] = 0.0;
      for (k = M1, j = 0; j < SIZE; j++)
      {
        y[i] += x[i+k+start-m-1] * wsinc[j][index];
        k++;
      }
    }
  }
  else
  /* *If delay is less than n/2 (=30), then lay down two pitch periods  */
  /* *with the appropriate delay, and part of third pitch period with	*/
  /* *the appropriate delay.						*/
  {
    pitch = m + d;
    if (m < n / 2)
    {
      for (i = 0; i < m + M1; i++)
      {
        y[i] = 0.0;
        for (k = M1, j = 0; j < SIZE; j++)
        {
          y[i] += x[i+k+start-m-1] * wsinc[j][index];
          k++;
        }
      }
      for (i = 0; i < M2 + 1; i++)
      {
        x[start+i-1] = y[i];
      }
      for (i = m + M1; i < m; i++)
      {
        y[i] = 0.0;
        for (k = M1, j = 0; j < SIZE; j++)
        {
          y[i] += x[i+k+start-m-1] * wsinc[j][index];
          k++;
        }
       }
      pitch2 = 2 * pitch;
      k2 = (int)(pitch2 - 2 * m);
      mtwo = (int)(pitch2);
      d = pitch2 - (2 * m + k2);
      if (fabs(d) > 0.1 && fabs(d) < 0.9)
        index = qd_nr(d);
      else 
        index = 5;
      for (i = m; i < mtwo; i++)
      {
        y[i] = 0.0;
        for (k = M1, j = 0; j < SIZE; j++)
        {
          y[i] += x[i+k+start-2*m-k2-1] * wsinc[j][index];
          k++;
        }
      }
      pitch3 = 3 * pitch;
      k3 = (int)(pitch3 - (3 * m + k2));
      d = pitch3 - (3 * m + k2 + k3);
      if (fabs(d) > 0.1 && fabs(d) < 0.9)
        index = qd_nr(d);
      else
        index = 5;
      for (i = mtwo; i < n; i++)
      {
        y[i] = 0.0;
         for (k = M1, j = 0; j < SIZE; j++)
         {
           y[i] += x[i+k+start-m-mtwo-k3-k2-1] * wsinc[j][index];
           k++;
         }
      }     
    }
    else
    /* *If delay is greater than n/2 (=30), then lay down one pitch	*/
    /* *period with the appropriate delay, and part of the second pitch	*/
    /* *period with appropriate delay.		  			*/
    {
      for (i = 0; i < m + M1; i++)
      {
        y[i] = 0.0;
        for (k = M1, j = 0; j < SIZE; j++)
        {
          y[i] += x[i+k+start-m-1] * wsinc[j][index];
          k++;
        }
      }
      for (i = 0; i < M2 + 1; i++)
        x[start+i-1] = y[i];
      for (i = m + M1; i < m; i++)
      {
        y[i] = 0.0;
        for (k = M1, j = 0; j < SIZE; j++)
        {
          y[i] += x[i+k+start-m-1] * wsinc[j][index];
          k++;
        }
      }
      pitch2 = 2 * pitch;
      k2 = (int)(pitch2 - 2 * m);
      d = pitch2 - (2 * m + k2);
      if (fabs(d) > 0.1 && fabs(d) < 0.9)
        index = qd_nr(d);
      else
        index = 5;
      for (i = m; i < n; i ++)
      {
        y[i] = 0.0;
        for (k = M1, j = 0; j < SIZE; j++)
        {
          y[i] += x[i+k+start-2*m-k2-1] * wsinc[j][index];
          k++;
        }
      }
    }
  }
  for (i = 0; i < M2 + 1; i++)
    x[start+i-1] = 0;
}