语音压缩

开发平台：
Unix_Linux

delay_nr.f：源码内容
							C==========================================================================
C
C ROUTINE
C		delay_nr
C
C FUNCTION
C		Time delay a bandlimited signal
C		using blockwise recursion (aka "nonrecursive")
C
C SYNOPSIS
C		subroutine delay_nr(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
C	bufr		int	i	length of input sequence x(n)
C	start		int	i	beginning of output sequence
C	n		int	i	length of input 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**************************************************************************
C
	subroutine delay_nr(x, bufr, start, n, d, m, y)
	implicit undefined(a-z)
	include 'ccsub.h'
convex	#include "ccsub.h"
	integer n, bufr, start
	real x(bufr), y(0:n-1), d, pitch, pitch2, pitch3
	integer nfrac, i, j, k, k2, k3, m, mtwo, qd_nr, M1, M2
c	* 5 fractional delays plus zero delay = 6
	parameter (nfrac = 6, 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,    0/
	data frac/    0.25, 0.33333333, 0.5, 0.66666667, 0.75, 0.0/
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 allign 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		wsinc(.,6)	frac = 0   zero delay for programming ease
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_nr(d)
c
c	Resample:
c
c	If delay is greater than n (=60), then lay down n values at 
c	the appropriate delay.
	if (m .ge. n) then
	   do 80 i = 0, n-1+M1
	      y(i) = 0.0
	      do 70 j = M1, M2
	         y(i) = y(i) + x(i+j+start-m)*wsinc(j,k)
70	      continue
80	   continue
	   do 81 i=0,M2
	      x(start+i) = y(i)
81	   continue
      	   do 82 i = n+M1, n-1
	         y(i) = 0.0
	      do 83 j = M1, M2
	         y(i) = y(i) + x(i+j+start-m)*wsinc(j,k)
83	      continue
82	   continue
	else
c	If delay is less than n/2 (=30), then lay down two pitch 
c	periods with the appropriate delay, and part of third
c	pitch period with the appropriate delay.
	   pitch = m + d
	   if (m .lt. n/2) then
      	      do 90 i = 0, m-1+M1
	         y(i) = 0.0
	         do 91 j = M1, M2
	            y(i) = y(i) + x(i+j+start-m)*wsinc(j,k)
91	         continue
90	      continue
	      do 990 i=0,M2
	         x(start+i) = y(i)
990	      continue
      	      do 901 i = m+M1, m-1
	         y(i) = 0.0
	         do 911 j = M1, M2
	            y(i) = y(i) + x(i+j+start-m)*wsinc(j,k)
911	         continue
901	      continue
	      pitch2 = 2*pitch
	      k2 = int(pitch2 - 2*m)
	      mtwo = int(pitch2)
	      d = pitch2 - (2*m + k2)
	      if (abs(d) .gt. 0.1 .and. abs(d) .lt. 0.9) then
		 k = qd_nr(d)
	      else
		 k = 6
	      end if
      	      do 92 i = m, mtwo-1
	         y(i) = 0.0
	         do 93 j = M1, M2
	            y(i) = y(i) + x(i+j+start-2*m-k2)*wsinc(j,k)
93	         continue
92	      continue
	      pitch3 = 3*pitch
	      k3 = int(pitch3 - (3*m + k2))
	      d = pitch3 - (3*m + k2 + k3)
	      if (abs(d) .gt. 0.1 .and. abs(d) .lt. 0.9) then
		 k = qd_nr(d)
	      else
		 k = 6
	      end if
      	      do 94 i = mtwo, n-1
	         y(i) = 0.0
	         do 95 j = M1, M2
	            y(i) = y(i) + x(i+j+start-m-mtwo-k3-k2)*wsinc(j,k)
95	         continue
94	      continue
	   else
c	If delay is greater than n/2 (=30), then lay down one pitch
c	period with the appropriate delay, and part of second 
c	pitch period with appropriate delay.
      	      do 96 i = 0, m-1+M1
	         y(i) = 0.0
	         do 97 j = M1, M2
	            y(i) = y(i) + x(i+j+start-m)*wsinc(j,k)
97	         continue
96	      continue
	      do 991 i=0,M2
		 x(start+i)=y(i)
991	      continue
      	      do 960 i = m+M1, m-1
	         y(i) = 0.0
	         do 970 j = M1, M2
	            y(i) = y(i) + x(i+j+start-m)*wsinc(j,k)
970	         continue
960	      continue
	      pitch2 = 2*pitch
	      k2 = int(pitch2 - 2*m)
	      d = pitch2 - (2*m + k2)
	      if (abs(d) .gt. 0.1 .and. abs(d) .lt. 0.9) then
		 k = qd_nr(d)
	      else
		 k = 6
	      end if
      	      do 98 i = m, n-1
	         y(i) = 0.0
	         do 99 j = M1, M2
	            y(i) = y(i) + x(i+j+start-2*m-k2)*wsinc(j,k)
99	         continue
98	      continue
	   end if
	end if 
	do 992 i=0,M2
	   x(start+i)=0
992	continue
	return
	end
c
c	Sinc function
c
copt	function sinc(arg)
copt	implicit undefined(a-z)
copt	real sinc, arg, pi
copt	pi = 4.0*atan(1.0)
copt	if (arg .eq. 0.0) then
copt	   sinc = 1.0
copt	else
copt	   sinc = sin(pi*arg)/(pi*arg)
copt	end if
copt	return
copt	end
c
c	Quantize d function 
c
	function qd_nr(d)
	implicit undefined(a-z)
	real d
	integer qd_nr, i, k, nfrac
	parameter (nfrac = 5)
	real frac(nfrac)
	logical ok
C
C	* The second and fourth elements of the data statement below
C	  when multiplied by 3 must be greater than 1 or 2 respectively.
	data frac/0.25, 0.33334, 0.5, 0.66667, 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_nr=k
	return
	end
c

						