delay.c
上传用户:szhypcb168
上传日期:2007-01-06
资源大小:2187k
文件大小:5k
- #define NFRAC 5
- #define TRUE 1
- #define FALSE 0
- #define M1 -4
- #define M2 3
- /* five fractional delays calculated over an 8 point interpolation */
- /* (-4 to 3) */
- static float frac[NFRAC] = {0.25, 0.33333333, 0.5, 0.66666667, 0.75};
- static int twelfths[NFRAC] = {3, 4, 6, 8, 9};
- /**************************************************************************
- *
- * NAME
- * sinc
- *
- * FUNCTION
- *
- *
- * SYNOPSIS
- * subroutine sinc(arg)
- *
- * formal
- * data I/O
- * name type type function
- * -------------------------------------------------------------------
- * arg float i
- * y(n) float func
- *
- ***************************************************************************
- *
- * DESCRIPTION
- *
- * This interpolating (sinc) function is Hamming windowed to bandlimit
- * the signal to reduce aliasing.
- *
- ***************************************************************************
- *
- * CALLED BY
- *
- * pitchvq psearch
- *
- * CALLS
- *
- *
- ***************************************************************************/
- #include <math.h>
- #include <stdio.h>
- float
- sinc(arg)
- float arg;
- {
- float pi;
- pi = 4.0 * atan(1.0);
- if (arg == 0.0)
- return(1.0);
- else
- return(sin(pi * arg) / (pi * arg));
- }
- /**************************************************************************
- *
- * NAME
- * qd
- *
- * FUNCTION
- *
- *
- * SYNOPSIS
- * subroutine qd(d)
- *
- * formal
- * data I/O
- * name type type function
- * -------------------------------------------------------------------
- * d float i
- * qd(d) int o quantize d function
- *
- ***************************************************************************
- *
- * DESCRIPTION
- *
- * Quantize d function
- *
- ***************************************************************************
- *
- * CALLED BY
- *
- * delay
- *
- * CALLS
- *
- *
- ***************************************************************************/
- int
- qd(d)
- float d;
- {
- int i, index, ok;
- ok = FALSE;
- for (i = 0; i < NFRAC; i++)
- {
- if (fabs(d - frac[i]) < 1.e-2)
- {
- index = i;
- ok = TRUE;
- }
- }
- if (!ok)
- {
- fprintf(stderr, "delay: Invalid pitch delay = %fn", d);
- exit(1);
- }
- return(index);
- }
-
- /**************************************************************************
- *
- * NAME
- * delay
- *
- * FUNCTION
- * Time delay a bandlimited signal
- * using point-by-point recursion
- *
- * SYNOPSIS
- * subroutine delay(x, start, n, d, m, y)
- *
- * formal
- * data I/O
- * name type type function
- * -------------------------------------------------------------------
- * x[n] float i signal input (output in last 60)
- * 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
- *
- **************************************************************************/
- #define SIZE (M2 - M1 + 1)
- delay(x, start, n, d, m, y)
- float x[], d, y[];
- int m, n, start;
- {
- static float wsinc[SIZE][NFRAC], hwin[12*SIZE+1];
- static int first = TRUE;
- int i, j, k, index;
- /* 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) */
- 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(d);
- /* *Resample: */
- for (i = 0; i < n; i++)
- {
- x[start+i-1] = 0.0;
- for (k = M1, j = 0; j < SIZE; j++)
- {
- x[start+i-1] += x[start-m+i+k-1] * wsinc[j][index];
- k++;
- }
- }
- /* *The v0 array in psearch.c/pgain.c must be zero above "start" */
- /* *because of overlap and add convolution techniques used in pgain. */
- for (i = 0; i < n; i++)
- {
- y[i] = x[start+i-1];
- x[start+i-1] = 0.0;
- }
- }