cgain.c
上传用户:szhypcb168
上传日期:2007-01-06
资源大小:2187k
文件大小:10k
- /*************************************************************************
- *
- * ROUTINE
- * cgain
- *
- * FUNCTION
- * Find codeword gain and error (TERNARY CODE BOOK ASSUMED!)
- *
- * SYNOPSIS
- * function cgain(ex, l, first, len, match)
- *
- * formal
- *
- * data I/O
- * name type type function
- * -------------------------------------------------------------------
- * ex[l] float i excitation vector (ternary codeword)
- * l int i size of ex (dimension of codeword)
- * first int i first call flag
- * len int i length to truncate impulse response
- * match float o negative partial squared error
- * cgain float fun optimal gain for ex
- *
- * external
- * data I/O
- * name type type function
- * -------------------------------------------------------------------
- * e0[] float i
- * h[]; float i
- *
- **************************************************************************
- *
- * DESCRIPTION
- *
- * Find the gain and error for each code word:
- * a. Filter code words through impulse response
- * of perceptual weighting filter (LPC filter with
- * bandwidth broadening).
- * b. Correlate filtered result with actual second error
- * signal (e0).
- * c. Compute MSPE gain and error for code book vector.
- *
- * Notes: Code words may contain many zeros (i.e., ex(1)=0). The
- * code book could be accessed by a pointer to nonzero samples.
- * Because the code book is static, it`s silly to test its
- * samples as in the code below.
- *
- * Proper selection of the convolution length (len) depends on
- * the perceptual weighting filter's expansion factor (gamma)
- * which controls the damping of the inpulse response.
- *
- * This is one of CELP's most computationally intensive
- * routines. Neglecting overhead, the approximate number of
- * DSP instructions (add, multiply, multiply accumulate, or
- * compare) per second (IPS) is:
- *
- * Code book size MIPS
- * -------------- ----
- * 64 1.1
- * 128 2.1
- * 256 4.2
- * 512 8.3
- *
- * C: convolution (recursive truncated end-point correction)
- * R: correlation
- * E: energy (recursive end-point correction)
- * G: gain quantization
- *
- *
- * celp code book search complexity (doesn't fully exploit ternary values!):
- *
- *#define j 10
- *
- * DSP chip instructions/operations:
- * int MUL=1; !multiply
- * int ADD=1; !add
- * int SUB=1; !subtract
- * int MAC=1; !multiply & accumulate
- * int MAD=2; !multiply & add
- * int CMP=1; !compare
- *
- * CELP algorithm parameters:
- * int L=60; !subframe length
- * int len=30; !length to truncate calculations (<= L)
- * int K=4; !number of subframes/frame
- * int shift=2; !shift between code words
- * int g_bits=5; !cbgain bit allocation
- * float p=0.77; !code book sparsity
- * float F=30.e-3; !time (seconds)/frame
- *
- * int N[j] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512};
- *
- * main()
- * {
- * int i;
- * float C, R, E, G, IPS;
- * printf("n N C R E G MIPSn");
- * for (i = 0; i < j; i++)
- * {
- * C = (335)*MAD + (N[i]-1)*shift*(1.0-p)*len*ADD;
- * R = N[i]*L*MAC;
- * E = L*MAC + (N[i]-1)*((1.0-p*p)*L*MAC + (p*p)*2*MAD);
- * G = N[i]*(g_bits*(CMP+MUL+ADD) + 3*MUL+1*SUB);
- * IPS = (C+R+E+G)*K/F;
- * printf(" %4d %f %f %f %f %fn", 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);
- * }
- * }
- *
- * N C R E G MIPS
- * 1 0.089333 0.008000 0.008000 0.002533 0.107867
- * 2 0.091173 0.016000 0.011573 0.005067 0.123813
- * 4 0.094853 0.032000 0.018719 0.010133 0.155706
- * 8 0.102213 0.064000 0.033011 0.020267 0.219491
- * 16 0.116933 0.128000 0.061595 0.040533 0.347062
- * 32 0.146373 0.256000 0.118763 0.081067 0.602203
- * 64 0.205253 0.512000 0.233100 0.162133 1.112486
- * 128 0.323013 1.024000 0.461773 0.324267 2.133053
- * 256 0.558533 2.048000 0.919118 0.648533 4.174185
- * 512 1.029573 4.096000 1.833810 1.297067 8.256450
- *
- **************************************************************************
- *
- * CALLED BY
- *
- * cbsearch
- *
- * CALLS
- *
- * gainencode2
- *
- ***************************************************************************
- *
- * REFERENCES
- *
- * Tremain, Thomas E., Joseph P. Campbell, Jr and Vanoy C. Welch,
- * "A 4.8 kbps Code Excited Linear Predictive Coder," Proceedings
- * of the Mobile Satellite Conference, 3-5 May 1988, pp. 491-496.
- *
- * Campbell, Joseph P. Jr., Vanoy C. Welch and Thomas E. Tremain,
- * "The New 4800 bps Voice Coding Standard," Military and
- * Government Speech Tech, 1989, p. 64-70.
- *
- * Lin, Daniel, "New Approaches to Stochastic Coding of Speech
- * Sources at Very Low Bit Rates," Signal Processing III: Theories
- * and Applications (Proceedings of EUSIPCO-86), 1986, p.445.
- *
- * Xydeas, C.S., M.A. Ireton and D.K. Baghbadrani, "Theory and
- * Real Time Implementation of a CELP Coder at 4.8 and 6.0 kbits/s
- * Using Ternary Code Excitation," Fifth International Conference on
- * Digital Processing of Signals in Communications, 1988, p. 167.
- *
- * Ess, Mike, "Simple Convolution on the Cray X-MP,"
- * Supercomputer, March 1988, p. 35
- *
- * Supercomputer, July 1988, p. 24
- *
- **************************************************************************/
- #include <math.h>
- #include "ccsub.h"
- extern float e0[MAXLP], h[MAXLP];
- float
- cgain(ex, l, first, len, match)
- int l, first, len;
- float ex[], *match;
- {
- register float cor;
- float cgain, gainencode2();
- static float y[MAXL], y59save, y60save, eng;
- int i, j;
- if (first)
- {
- /** For first code word, calculate and save convolution
- of code word with truncated (to len) impulse response:
- NOTES: A standard convolution of two L point sequences
- produces 2L-1 points, however, this convolution
- generates only the first L points.
- A "scalar times vector accumulation" method is used
- to exploit (skip) zero samples of the code words:
- min(L-i, len-1)
- y = SUM ex * h , where i = 0, ..., L-1 points
- i+j, t j=0 i j
- ex |0 1 . . . L-1|
- h |x x len-1...1 0| = y[0]
- h |x x len-1...1 0| = y[1]
- : :
- h |x x len-1...1 0| = y[L-1]
- ----------------------------------------------------- */
- for (i = 0; i < l; i++)
- y[i] = 0.0;
- for (i = 0; i < l; i++)
- {
- if (nint(ex[i]) != 0.0)
- for(j = 0; j < l-i && j < len; j++)
- y[i+j] += ex[i]*h[j];
- }
- }
- else
- {
- /** End correct the convolution sum on subsequent code words:
- (Do two end corrections for a shift by 2 code book)
- y = 0
- 0, 0
- y = y + ex * h where i = 0, ..., L-1 points
- i, m i-1, m-1 -m i and m = 0, ..., cbsize-1 code words
- NOTE: The data movements in the 2 loops with "y[i] = y[i - 1]"
- are performed many times and can be quite time consuming.
- Therefore, special precautions should be taken when
- implementing this. Some implementation suggestions:
- 1. Circular buffers with pointers to eliminate data moves.
- 2. Fast "block move" operation as offered on some DSPs.
- ----------------------------------------------------- */
- /* *First shift */
- if (nint(ex[1]) != 0)
- {
- /* *ternary stochastic code book (-1, 0, +1) */
- if (nint(ex[1]) == 1)
- for (i = len - 1; i > 0; i--)
- y[i - 1] += h[i];
- else
- for (i = len - 1; i > 0; i--)
- y[i - 1] -= h[i];
- }
- for (i = l - 1; i > 0; i--)
- y[i] = y[i - 1];
- y[0] = ex[1] * h[0];
- /* *Second shift */
- if (nint(ex[0]) != 0)
- {
- /* *ternary stochastic code book (-1, 0, +1) */
- if (nint(ex[0]) == 1)
- for (i = len - 1; i > 0; i--)
- y[i - 1] += h[i];
- else
- for (i = len - 1; i > 0; i--)
- y[i - 1] -= h[i];
- }
-
- for (i = l - 1; i > 0; i--)
- y[i] = y[i - 1];
- y[0] = ex[0] * h[0];
- }
- /** Calculate correlation and energy:
- e0 = spectrum & pitch prediction residual
- y = error weighting filtered code words
- /// CELP's computations are focused in this correlation ///
- - For a 512 code book this correlation takes 4 MIPS!
- - Decimation?, Down-sample & decimate?, FEC codes? */
- cor = 0.0;
- for (i = 0; i < l; i++)
- {
- cor += y[i] * e0[i];
- }
- /* *End correct energy on subsequent code words: */
- if (nint(ex[0]) == 0 && nint(ex[1]) == 0 && !first)
- eng = eng - y59save * y59save - y60save * y60save;
- else
- {
- eng = 0.0;
- for (i = 0; i < l; i++)
- eng += y[i] * y[i];
- }
- y59save = y[l - 2];
- y60save = y[l - 1];
- /* Compute gain and error:
- NOTE: Actual MSPE = e0.e0 - cgain(2*cor-cgain*eng)
- since e0.e0 is independent of the code word,
- minimizing MSPE is equivalent to maximizing:
- *match = cgain(2*cor-cgain*eng)
- If unquantized cgain is used, this simplifies:
- *match = cor*cgain */
- /* Independent (open-loop) quantization of gain and match (index): */
- if (eng <= 0.0) eng = 1.0;
- cgain = cor / eng;
- *match = cor * cgain;
- /* *Joint (closed-loop) quantization of gain and match (index): */
- /* *(Beware that the match score can be negative!) */
- /* cgain = gainencode2(cor, eng, &i); */
- /* *match = cgain * (2.0 * cor - cgain * eng); */
- return (cgain);
- }