语音压缩

开发平台：
Unix_Linux

cgain.f：源码内容
							C==========================================================================
C
C ROUTINE
C		cgain
C
C FUNCTION
C		Find codeword gain and error (TERNARY CODE BOOK ASSUMED!)
C
C SYNOPSIS
C               function cgain(ex, l, first, len, match)
C
C   formal 
C
C                       data    I/O
C       name            type    type    function
C       -------------------------------------------------------------------
C       ex(l)		real	i	excitation vector (ternary codeword)
C       l		int	i	size of ex (dimension of codeword)
C       first		logical	i	first call flag
C	len		int	i	length to truncate impulse response
C       match		real	o	negative partial squared error
C       cgain		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 code word find its gain and error:
C	   a.  Filter code words through impulse response
C	       of perceptual weighting filter (LPC filter with
C	       bandwidth broadening).
C	   b.  Correlate filtered result with actual second error
C	       signal (e0).
C	   c.  Compute MSPE gain and error for code book vector.
C
C	Notes:  Codewords may contain many zeros (i.e., ex(1)=0).  The
C		code book could be accessed by a pointer to nonzero samples.
C		Because the code book is static, it`s silly to test its
C		samples as in the code below.
C
C		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	               Code book size   MIPS
C	               --------------   ----
C	                   64           1.1
C	                  128           2.1
C	                  256           4.2
C	                  512 (max)     8.3
C
C	        C:  convolution (recursive truncated end-point correction)
C	        R:  correlation
C	        E:  energy (recursive end-point correction)
C	        G:  gain quantization
C
C	celp code book search complexity (doesn't fully exploit ternary values!):
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, shift, g_bits
C	real p, F
C	parameter (L=60)	!subframe length
C	parameter (len=30)	!length to truncate calculations (<= L)
C	parameter (p=0.77)	!code book sparsity
C	parameter (shift=2)	!shift between code words
C	parameter (K=4)		!number of subframes/frame
C	parameter (F=30.e-3)	!time (seconds)/frame
C	parameter (g_bits=5)	!cbgain bit allocation
c
C	integer j
C	parameter (j=10)
C	integer N(j), i
C	real C, R, E, G, IPS
C	data N/1, 2, 4, 8, 16, 32, 64, 128, 256, 512/
C	print 1
C1	format(10x,'N',10x,'C',13x,'R',12x,'E',15x,'G',13x,'MIPS')
C	do i = 1, j
C	   C = (335)*MAD + (N(i)-1)*shift*(1.-p)*len*ADD
C	   R = N(i)*L*MAC
C	   E = L*MAC + (N(i)-1)*((1.-p*p)*L*MAC + (p*p)*2*MAD)
C	   G = N(i)*(g_bits*(CMP+MUL+ADD) + 3*MUL+1*SUB)
C	   IPS = (C+R+E+G)*K/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	end do
C	end
C
C    N        C             R             E               G           MIPS
C    1  8.9333333E-02  8.0000004E-03  8.0000004E-03  2.5333334E-03  0.1078667   
C    2  9.1173328E-02  1.6000001E-02  1.1573013E-02  5.0666668E-03  0.1238130   
C    4  9.4853334E-02  3.2000002E-02  1.8719040E-02  1.0133334E-02  0.1557057   
C    8  0.1022133      6.4000003E-02  3.3011094E-02  2.0266667E-02  0.2194911   
C   16  0.1169333      0.1280000      6.1595201E-02  4.0533334E-02  0.3470618   
C   32  0.1463733      0.2560000      0.1187634      8.1066668E-02  0.6022034   
C   64  0.2052534      0.5120000      0.2330998      0.1621333       1.112487   
C  128  0.3230133       1.024000      0.4617727      0.3242667       2.133053   
C  256  0.5585334       2.048000      0.9191184      0.6485333       4.174185   
C  512   1.029573       4.096000       1.833810       1.297067       8.256450  
C
C==========================================================================
C
C CALLED BY
C
C	cbsearch
C
C CALLS
C
C	gainencode2
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	"The New 4800 bps Voice Coding Standard," Military and 
C	Government Speech Tech, 1989, p. 64-70.
C
C	Lin, Daniel, "New Approaches to Stochastic Coding of Speech
C	Sources at Very Low Bit Rates," Signal Processing III:  Theories
C	and Applications (Proceedings of EUSIPCO-86), 1986, p.445.
C
C	Xydeas, C.S., M.A. Ireton and D.K. Baghbadrani, "Theory and
C	Real Time Implementation of a CELP Coder at 4.8 and 6.0 kbits/s
C	Using Ternary Code Excitation," Fifth International Conference on
C	Digital Processing of Signals in Communications, 1988, p. 167.
C
C	Ess, Mike, "Simple Convolution on the Cray X-MP,"
C	Supercomputer, March 1988, p. 35
C
C	Supercomputer, July 1988, p. 24
C
C==========================================================================
C*-
	function cgain(ex, l, first, len, match)
	implicit undefined(a-z)
	integer l, len
	real ex(l), match, cgain
	logical first
	include 'ccsub.com'
convex	#include "ccsub.com"
	real cor, eng, y(maxl), y59save, y60save, gainencode2
	integer i, j
c
	save y, y59save, y60save, eng
c
	if (first) then
c
c	   For first code word, calculate and save convolution
c	   of code word with truncated (to len) impulse response:
c
c	   NOTES: A standard convolution of two L point sequences
c	          produces 2L-1 points, however, this convolution
c	          generates only the first L points.
c
c	          A "scalar times vector accumulation" method is used
c	          to exploit (skip) zero samples of the code words:
c
c	          min(L-i+1, len)
c	   y         =  SUM  ex * h   , where i = 1, ..., L points
c	    i+j-1, t    j=1    i   j
c
c	                ex |1 2  .  .  .  L|
c	   h |x x len ... 2 1|               = y(1)
c	     h |x x len ... 2 1|             = y(2)
c	                      :                 :
c	                 h |x x len ... 2 1| = y(L)
c
	   do 10 i = 1, l
	      y(i) = 0.0
10	   continue
	   do 40 i = 1, l
	      if (nint(ex(i)) .ne. 0) then
	         do 30 j = 1, min(l-i+1, len)
	            y(i+j-1) = y(i+j-1) + ex(i)*h(j)
30	         continue
	      end if
40	   continue
	else
c
c	   End correct the convolution sum on subsequent code words:
c	   (Do two end corrections for a shift by 2 code book)
c
c	   y     =  0
c	    0, 0
c	   y     =  y        + ex * h   where i = 1, ..., L points
c	    i, m     i-1, m-1   -m   i  and   m = 1, ..., cbsize-1 code words
c
c	   NOTE:  The data movements in loops "do 59 ..." and "do 69 ..."
c	          are performed many times and can be quite time consuming.
c	          Therefore, special precautions should be taken when
c	          implementing this.  Some implementation suggestions:
c	          1.  Circular buffers with pointers to eliminate data moves.
c	          2.  Fast "block move" operation as offered on some DSPs.
c
c
c	   First shift
c
	   if (nint(ex(2)) .ne. 0) then
c	      *ternary stochastic code book (-1, 0, +1)
	      if (nint(ex(2)) .eq. 1) then
	         do 50 i = len, 2, -1
	            y(i-1) = y(i-1) + h(i)
50	         continue
	      else
	         do 55 i = len, 2, -1
	            y(i-1) = y(i-1) - h(i)
55	         continue
	      end if
	   end if
	   do 59 i = l, 2, -1
	      y(i) = y(i-1)
59	   continue
	   y(1) = ex(2)*h(1)
c
c	   Second shift
c
	   if (nint(ex(1)) .ne. 0) then
c	      *ternary stochastic code book (-1, 0, +1)
	      if (nint(ex(1)) .eq. 1) then
	         do 60 i = len, 2, -1
	            y(i-1) = y(i-1) + h(i)
60	         continue
	      else
	         do 65 i = len, 2, -1
	            y(i-1) = y(i-1) - h(i)
65	         continue
	      end if
	   end if
	   do 69 i = l, 2, -1
	      y(i) = y(i-1)
69	   continue
	   y(1) = ex(1)*h(1)
	end if
c
c	Calculate correlation and energy:
c	e0 = spectrum & pitch prediction residual
c	y  = error weighting filtered code words
c
c	///  CELP's computations are focused in this correlation ///
c		- For a 512 code book this correlation takes 4 MIPS!
c		- Decimation?, Down-sample & decimate?, FEC codes?
	cor = 0.0
	do 70 i = 1, l
	   cor = cor + y(i)*e0(i)
70	continue
c
c			*End correct energy on subsequent code words:
	if (nint(ex(1)).eq.0 .and. nint(ex(2)).eq.0 .and. .not.first) then
	   eng = eng - y59save*y59save-y60save*y60save
	else
	   eng = 0.0
	   do 80 i = 1, l
	      eng = eng + y(i)*y(i)
80	   continue
	end if
	y59save = y(l-1)
	y60save = y(l)
c
c	Compute gain and error:
c	  NOTE: Actual MSPE = e0.e0 - cgain(2*cor-cgain*eng)
c		since e0.e0 is independent of the code word,
c		minimizing MSPE is equivalent to maximizing:
c		     match = cgain(2*cor-cgain*eng)
c		If unquantized cgain is used, this simplifies:
c		     match = cor*cgain
c
c	Independent (open-loop) quantization of gain and match (index):
c
	if (eng .le. 0.0) eng = 1.0
	cgain = cor/eng
	match = cor*cgain
c
c	Joint (closed-loop) quantization of gain and match (index):
c	(Beware that the match score can be negative!)
c
c>>>	cgain = gainencode2(cor, eng, cbgbits, cbgtype, i)
c>>>	match = cgain*(2.*cor - cgain*eng)
c
	return
	end

						