cgain.m

FS1016源代码

% MATLAB SIMULATION OF NSA FS-1016 CELP v3.2
% COPYRIGHT (C) 1995-99 ANDREAS SPANIAS AND TED PAINTER
%
% This Copyright applies only to this particular MATLAB implementation
% of the FS-1016 CELP coder.  The MATLAB software is intended only for educational
% purposes.  No other use is intended or authorized.  This is not a public
% domain program and distribution to individuals or networks is strictly
% prohibited.  Be aware that use of the standard in any form is goverened
% by rules of the US DoD.  Therefore patents and royalties may apply to
% authors, companies, or committees associated with this standard, FS-1016.  For
% questions regarding the MATLAB implementation please contact Andreas
% Spanias at  [phone].  For questions on rules,
% royalties, or patents associated with the standard, please contact the DoD.
%
% ALL DERIVATIVE WORKS MUST INCLUDE THIS COPYRIGHT NOTICE.
%
% ******************************************************************
% CGAIN
%
% PORTED TO MATLAB FROM CELP 3.2a C RELEASE
% 7-11-94
%
% ******************************************************************
%
% DESCRIPTION
%
% Find codeword gain and error (Ternary codebook assumed)
%
% DESIGN NOTES
%
% 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.
%
% 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 impulse 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):
%
%     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
%
% REFERENCES
%
% 1. 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.
%
% 2. 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.
%
% 3. 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.
%
% 4. 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.
%
% 5. Ess, Mike, "Simple Convolution on the Cray X-MP,"
%    Supercomputer, March 1988, p. 35
%
% 6. Supercomputer, July 1988, p. 24
%
% VARIABLES
%
% INPUTS
%   ex         -     Excitation vector (ternary codeword)
%   l          -     Size of ex (dimension of codeword)
%   first      -     First call flag
%   len        -     Length of truncated impulse response
%
% OUTPUTS
%   Cgain      -     Optimal gain for codeword ex
%   match      -     Negative partial squared error
%
% INTERNALS
%   i          -     Convolution output sequence index, y(i)
%   jmin       -     Lower bound on convolution summation index
%   jmax       -     Upper bound on convolution summation index
%
% GLOBALS
%   Ycg        -     Output of perceptual weighting filter
%   y59save    -     2nd to last sample, Ycg
%   y60save    -     Last sample, Ycg
%   Engcg      -     Codeword energy
%   h          -     Impulse response, perceptual weighting filter
%   e0         -     First error sequence
%
% CONSTANTS
%   TRUE       -     Boolean logic constant
%   FALSE      -     Boolean logic constant
%
% ******************************************************************

function [ Cgain, match ] = cgain( ex, l, first, len )

% DECLARE GLOBAL VARIABLES
global Ycg y59save y60save Engcg h e0

% DECLARE GLOBAL CONSTANTS
global TRUE FALSE

if first == TRUE
    % FOR FIRST CODE WORD, CALCULATE AND SAVE CONVOLUTION OF CODE WORD WITH
    % TRUNCATED (TO LEN) IMPULSE RESPONSE.
    %
    % 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]
    Ycg = zeros( l, 1 );
    for i = 0:l-1
	    if round( ex(i+1) ) ~= 0.00
	        jmin = 1;  jmax = min( l-i, len );
	        Ycg(i+jmin:i+jmax) = Ycg(i+jmin:i+jmax) + ...
                                 ( ex(i+1) * h(jmin:jmax) );
        end
    end
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
    %
    % 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 round( ex(2) ) ~= 0

        % TERNARY STOCHASTIC CODEBOOK (-1,0,1)
        if round( ex(2) ) == 1
            Ycg(len-1:-1:1) = Ycg(len-1:-1:1) + h(len:-1:2);
        else
            Ycg(len-1:-1:1) = Ycg(len-1:-1:1) - h(len:-1:2);
        end
    end
    Ycg(l:-1:2) = Ycg(l-1:-1:1);
    Ycg(1) = ex(2) * h(1);

    % SECOND SHIFT
    if round( ex(1) ) ~= 0

        % TERNARY STOCHASTIC CODEBOOK (-1,0,1)
        if round( ex(1) ) == 1
            Ycg(len-1:-1:1) = Ycg(len-1:-1:1) + h(len:-1:2);
        else
            Ycg(len-1:-1:1) = Ycg(len-1:-1:1) - h(len:-1:2);
        end
    end
    Ycg(l:-1:2) = Ycg(l-1:-1:1);
    Ycg(1) = ex(1) * h(1);
end

% CALCULATE CORRELATION AND ENERGY
% E0 = SPECTRUM AND PITCH PREDICTION RESIDUAL
% YCG = ERROR WEIGHTING FILTERED CODEWORDS
%
% NOTE:  CELP'S COMPUTATIONS ARE FOCUSED IN THIS CORRELATION -
%        FOR A 512 ELEMENT CODEBOOK THIS CORRELATION TAKES 4 MIPS.
%
cor = sum( Ycg .* e0 );

% END CORRECT ENERGY ON SUBSEQUENT CODE WORDS
if ( round(ex(1)) == 0 ) & ( round(ex(2)) == 0 ) & ( first ~= TRUE )
    Engcg = Engcg - ( y59save * y59save ) - ( y60save * y60save );
else
    Engcg = sum( Ycg .* Ycg );
end
y59save = Ycg(l-1);
y60save = Ycg(l);

% COMPUTE GAIN AND ERROR
%
%   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
%
if Engcg <= 0.00
    Engcg = 1.00;
end

% INDEPENDENT (OPEN-LOOP) QUANTIZATION OF GAIN AND MATCH (INDEX)
% IS USED HERE
Cgain = cor / Engcg;
match = cor * Cgain;