语音压缩

开发平台：
Unix_Linux

pgain.f：源码内容
							C==========================================================================
C
C ROUTINE
C		pgain
C
C FUNCTION
C		Find pitch gain and error
C
C SYNOPSIS
C               function pgain(ex, l, first, m, len, match)
C
C   formal 
C
C                       data    I/O
C       name            type    type    function
C       -------------------------------------------------------------------
C       ex(l)		real	i	excitation vector
C       l		int	i	size of ex
C       first		logical	i	first call flag
C	m		int	i	pitch lag
C	len		int	i	length to truncate impulse response
C       match		real	o	negative partial squared error
C       pgain		real	fun	optimal gain for ex
C
C   global 
C                       data    I/O
C       name            type    type    function
C       -------------------------------------------------------------------
C 	/ccsub/		see description include file
C
C==========================================================================
C	
C DESCRIPTION
C
C	For each lag:
C	   a.  Filter first error signal (v0) through truncated
C	       impulse response of perceptual weighting filter
C	       (LPC filter with bandwidth broadening).
C	   b.  Correlate filtered result with actual first error
C	       signal (e0).
C	   c.  Compute first order pitch filter coefficient (pgain)
C	       and error (er) for each lag.
C
C	Note:  Proper selection of the convolution length (len) depends on
C	       the perceptual weighting filter's expansion factor (gamma)
C	       which controls the damping of the impulse response.
C
C	        This is one of CELP's most computationally intensive
C	        routines.  Neglecting overhead, the approximate number of
C		DSP instructions (add, multiply, multiply accumulate, or
C		compare) per second (IPS) is:
C
C	        C:  convolution (recursive truncated end-point correction)
C               C': convolution (recursive truncated end-point correction)
C               R = E:  full correlation & energy
C               R'= E': delta correlation & energy
C               G:      gain quantization
C               G':     delta gain quantization
C
C	        IPS = 2.34 M (for integer delays)
C
c	pitch search complexity for integer delays:
C	implicit undefined(a-z)
c
c	DSP chip instructions/operation:
C	integer MUL, ADD, SUB, MAC, MAD, CMP
C	parameter (MUL=1)	!multiply
C	parameter (ADD=1)	!add
C	parameter (SUB=1)	!subtract
C	parameter (MAC=1)	!multiply & accumulate
C	parameter (MAD=2)	!multiply & add
C	parameter (CMP=1)	!compare
c
c	CELP algorithm parameters:
C	integer L, len, K, Kp, Np, dmin
C	real F
C	parameter (L=60)	!subframe length
C	parameter (len=30)	!length to truncate calculations (<= L)
C	parameter (K=2)		!number of full search subframes
C	parameter (Kp=2)	!number of delta search subframes
C	parameter (F=30.e-3)	!time (seconds)/frame
C	parameter (Np=32)	!number of delta subframe delays
C	parameter (dmin=20)	!minimum delay
c
C	integer j
C	parameter (j=4)
C	integer N(j), i
C	real C, R, E, G, Cp, Rp, Ep, Gp, IPS
C	data N/32, 64, 128, 256/
C	print 1
C1	format(10x,'N',10x,'C',13x,'R',12x,'E',15x,'G',13x,'MIPS')
C	do i = 1, j
C	   C = (len*(len+1)/2 + len*len)*MAC + (N(i)-1)*len*MAD
C     &	     + min(L-dmin,N(i))*(L-dmin)*ADD + (L/2-dmin)*(L-2*dmin)*ADD
C	   Cp= (len*(len+1)/2 + len*len)*MAC + (Np-1)*len*MAD
C     &	     + min(L-dmin,Np)*(L-dmin)*ADD + (L/2-dmin)*(L-2*dmin)*ADD
C	   R = N(i)*L*MAC
C	   Rp= Np*L*MAC
C	   E = R
C	   Ep= Rp
C	   G = N(i)*(1*CMP+1*MUL + 1*MUL)
C	   Gp= Np*(1*CMP+1*MUL + 1*MUL)
C	   IPS = ((C+R+E+G)*K + (Cp+Rp+Ep+Gp)*Kp)/F
C	   print *,N(i),C*K/1.e6/F,R*K/1.e6/F,E*K/1.e6/F,G*K/1.e6/F,IPS/1.e6
C	   print *,Np,Cp*Kp/1.e6/F,Rp*Kp/1.e6/F,Ep*Kp/1.e6/F,Gp*Kp/1.e6/F
C	end do
C	end
C
C  N          C             R            E               G             MIPS
C  32  0.3136667      0.1280000      0.1280000      6.4000003E-03   1.152133    
C  32  0.3136667      0.1280000      0.1280000      6.4000003E-03
C  64  0.4630000      0.2560000      0.2560000      1.2800001E-02   1.563867    
C  32  0.3136667      0.1280000      0.1280000      6.4000003E-03
C 128  0.7190000      0.5120000      0.5120000      2.5600001E-02   2.344667    
C  32  0.3136667      0.1280000      0.1280000      6.4000003E-03
C 256   1.231000       1.024000       1.024000      5.1200002E-02   3.906267    
C  32  0.3136667      0.1280000      0.1280000      6.4000003E-03
C
C==========================================================================
C
C CALLED BY
C
C	psearch
C
C CALLS
C
C
C
C==========================================================================
C
C REFERENCES
C
C	Tremain, Thomas E., Joseph P. Campbell, Jr and Vanoy C. Welch,
C	"A 4.8 kbps Code Excited Linear Predictive Coder," Proceedings
C	of the Mobile Satellite Conference, 3-5 May 1988, pp. 491-496.
C
C	Campbell, Joseph P. Jr., Vanoy C. Welch and Thomas E. Tremain,
C	"An Expandable Error-Protected 4800 bps CELP Coder (U.S. Federal
C	Standard 4800 bps Voice Coder)," Proceedings of ICASSP, 1989.
C	(and Proceedings of Speech Tech, 1989.)
C
C==========================================================================
C*-
	function pgain(ex, l, first, m, len, match)
	implicit undefined(a-z)
	integer l, m, len
	real ex(l), match, pgain
	logical first
	include 'ccsub.com'
convex	#include "ccsub.com"
	real cor, eng, y(maxlp), y2(maxlp)
	integer i, j
c
	save y
c
	if (first) then
c
c	   Calculate and save convolution of truncated (to len)
c	   impulse response for first lag of t (=mmin) samples:
c
c	        min(i, len)
c	   y     =  SUM  h * ex       , where i = 1, ..., L points
c	    i, t    j=1   j    i-j+1
c
c	              h |1 2...len x x|
c	   ex |L . . . 2 1|             = y(1)
c	     ex |L . . . 2 1|           = y(2)
c	                   :              :
c	             ex |L  . . .  2 1| = y(L)
c
	   do 40 i = 1, l
	      y(i) = 0.0
	      do 30 j = 1, min(i, len)
	         y(i) = y(i) + h(j)*ex(i-j+1)
30	      continue
40	   continue
	else
c
c	   End correct the convolution sum on subsequent pitch lags:
c
c	   y     =  0
c	    0, t
c	   y     =  y        + ex  * h   where i = 1, ..., L points
c	    i, m     i-1, m-1   -m    i  and   m = t+1, ..., tmax lags
c
	   do 50 i = len, 2, -1
	      y(i-1) = y(i-1) + ex(1)*h(i)
50	   continue
	   do 51 i = l, 2, -1
	      y(i) = y(i-1)
51	   continue
	   y(1) = ex(1)*h(1)
	end if
c
c	For lags (m) shorter than frame size (l), replicate the short
c	adaptive codeword to the full codeword length by
c	overlapping and adding the convolutions:
c
	do 60 i = 1, l
	   y2(i) = y(i)
60	continue
	if (m .lt. l) then
c	                            add in 2nd convolution
	   do 65 i = m+1, l
	      y2(i) = y(i) + y(i-m)
65	   continue
	   if (m .lt. l/2) then
c	                            add in 3rd convolution
	      do 69 i = 2*m+1, l
		 y2(i) = y2(i) + y(i-2*m)
69	      continue
	   end if
	end if
c
c	Calculate correlation and energy:
c	e0 = r(n)   = spectrum prediction residual
c	y2 = r(n-m) = error weighting filtered reconstructed
c	              pitch prediction signal (m = correlation lag)
c
	cor = 0.
	eng = 0.
	do 70 i = 1, l
	   cor = cor + y2(i)*e0(i)
	   eng = eng + y2(i)*y2(i)
70	continue
c
c	Compute gain and error:
c	  NOTE: Actual MSPE = e0.e0 - pgain(2*cor-pgain*eng)
c		since e0.e0 is independent of the code word,
c		minimizing MSPE is equivalent to maximizing:
c		     match = pgain(2*cor-pgain*eng)   (1)
c		If unquantized pgain is used, this simplifies:
c		     match = cor*pgain
c
c	  NOTE: Inferior results were obtained when quantized
c		pgain was used in equation (1)???
c
c	  NOTE: When the delay is less than the frame length, "match"
c	        is only an approximation to the actual error.
c
c	Independent (open-loop) quantization of gain and match (index):
c
	if (eng .le. 0.) eng = 1.0
	pgain = cor/eng
	match = cor*pgain
c
	return
	end

						