语音压缩

开发平台：
Unix_Linux

delay.f：源码内容
							C==========================================================================
C
C ROUTINE
C		delay
C
C FUNCTION
C		Time delay a bandlimited signal
C		using point-by-point recursion
C
C SYNOPSIS
C		subroutine delay(x, bufr, start, n, d, m, y)
C
C   formal 
C                       data    I/O
C       name            type    type    function
C       -------------------------------------------------------------------
C	x(n)		real	i	signal input (output in last 60)
C	bufr		int	i	length of input sequence x(n)
C	start		int	i	beginning of output sequence
C	n		int	i	length of output sequence
C	d		real	i	fractional pitch
C	m		int	i	integer pitch
C	y(n)		real	o	delayed input signal
C==========================================================================
C	
C DESCRIPTION
C
C	Subroutine to time delay a bandlimited signal by resampling the
C	reconstructed data (aka sinc interpolation).  The well known
C	reconstruction formula is:
C
C              |    M2      sin[pi(t-nT)/T]    M2
C   y(n) = X(t)| = SUM x(n) --------------- = SUM x(n) sinc[(t-nT)/T]
C              |   n=M1         pi(t-nT)/T    n=M1
C              |t=n+d
C
C	The interpolating (sinc) function is Hamming windowed to bandlimit
C	the signal to reduce aliasing.
C
C	Multiple simultaneous time delays may be efficiently calculated
C	by polyphase filtering.  Polyphase filters decimate the unused
C	filter coefficients.  See Chapter 6 in C&R for details. 
C==========================================================================
C	
C REFERENCES
C
C	Crochiere & Rabiner, Multirate Digital Signal Processing,
C	P-H 1983, Chapters 2 and 6.
C
C	Kroon & Atal, "Pitch Predictors with High Temporal Resolution,"
C	ICASSP '90, S12.6
C
C**************************************************************************
C
	subroutine delay(x, bufr, start, n, d, m, y)
	implicit undefined(a-z)
	integer n, bufr, start
	real x(bufr), y(0:n-1), d
	integer nfrac, i, j, k, m, qd, M1, M2
c	* 5 fractional delays calculated over an 8 point interpolation (-4 to 3)
	parameter (nfrac = 5, M1 = -4, M2 = 3)
	integer twelfths(nfrac)
	real sinc, wsinc(M1:M2, nfrac), hwin(12*M1:12*(M2+1)), frac(nfrac)
	logical first
c
	save first, hwin, wsinc
c
	data first /.true./
	data twelfths/3,    4,          6,   8,          9/
	data frac/    0.25, 0.33333333, 0.5, 0.66666667, 0.75/
c
c	Generate Hamming windowed sinc interpolating function
c	for each allowable fraction.  The interpolating functions
c	are stored in time reverse order (i.e., delay appears as
c	advance) to align with the adaptive code book v0 array.
c	The interp filters are:
c		wsinc(.,1)	frac = 1/4 (3/12)
c		wsinc(.,2)	frac = 1/3 (4/12)
c		.		.
c		wsinc(.,5)	frac = 3/4 (9/12)
c
	if (first) then
	   call ham(hwin, 12*(M2-M1+1)+1)
	   do 20 i = 1, nfrac
	      do 10 j = M1, M2
	         wsinc(j,i) = sinc(frac(i)+j)*hwin(12*j+twelfths(i))
10	      continue
20	   continue
	   first = .false.
	end if
c
	k=qd(d)
c
c	Resample:
	do 80 i = 0, n-1
	   x(start+i) = 0.0
	   do 70 j = M1, M2
	      x(start+i) = x(start+i) + x(start-m+i+j)*wsinc(j,k)
70	   continue
80	continue
c
c	The v0 array in psearch.F/pgain.f must be zero above "start"
c	because of overlap and add convolution techniques used in pgain.
c
	do 100 i = 0, n-1
	   y(i) = x(start + i)
	   x(start + i) = 0
100	continue
	return
	end
c
c	Sinc function
c
	function sinc(arg)
	implicit undefined(a-z)
	real sinc, arg, pi
	pi = 4.0*atan(1.0)
	if (arg .eq. 0.0) then
	   sinc = 1.0
	else
	   sinc = sin(pi*arg)/(pi*arg)
	end if
	return
	end
c
c	Quantize d function 
c
	function qd(d)
	implicit undefined(a-z)
	real d
	integer qd, i, k, nfrac
	parameter (nfrac = 5)
	real frac(nfrac)
	logical ok
	data frac/0.25, 0.33333333, 0.5, 0.66666667, 0.75/
c
	ok = .false.
	do 50 i = 1, nfrac
	   if (abs(d-frac(i)) .lt. 1.e-2) then
	      k = i
	      ok = .true.
	   end if
50	continue
	if (.not. ok) then
	   print *,'delay.f:  Invalid pitch delay = ',d
	   stop
	end if
	qd=k
	return
	end
c

						