% 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 (602) 965-1837.  For questions on rules,
% royalties, or patents associated with the standard, please contact the DoD.
%
% ALL DERIVATIVE WORKS MUST INCLUDE THIS COPYRIGHT NOTICE.
%
% ******************************************************************
% PGAIN
%
% PORTED TO MATLAB FROM CELP 3.2a C RELEASE
% 7-1-94
%
% ******************************************************************
%
% DESCRIPTION
%
% Find pitch gain and error
%
% DESIGN NOTES
%
% For each lag:
%
% a.  Filter first error signal (v0) through truncated
%     impulse response of perceptual weighting filter
%     (LPC filter with bandwidth broadening).
% b.  Correlate filtered result with actual first error
%     signal (e0).
% c.  Compute first order pitch filter coefficient (Pgain)
%     and error (er) for each lag.
%
% 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:
%
%   C      : convolution (recursive truncated end-point correction)
%   C'     : convolution (recursive truncated end-point correction)
%   R  = E : full correlation & energy
%   R' = E': delta correlation & energy
%   G      : gain quantization
%   G'     : delta gain quantization
%
%       IPS = 2.34 M (for integer delays)
%
%   i.e.,  L = 60, N = 128 pitch lags, N'= 32 delta delays
%          K = K'= 2 pitch updates/frame, and F=30 ms frame rate:
%          C = 9450, C'= 3690, R = E = 7680, R'= E'= 1920
%
%   IPS = 2.2 M
%
% Pitch search complexity for integer delays:
%
%   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,
%    "An Expandable Error-Protected 4800 bps CELP Coder (U.S. Federal
%    Standard 4800 bps Voice Coder)," Proceedings of ICASSP, 1989.
%    (and Proceedings of Speech Tech, 1989.)
%
% VARIABLES
%
% INPUTS
%   ex         -     Excitation vector
%   l          -     Size of excitation
%   first      -     First call flag
%   m          -     Pitch lag
%   len        -     Length to truncate impulse response
%
% OUTPUTS
%   match      -     Negative partial squared error
%   Pgain      -     Optimal gain for excitation (ex)
%
% INTERNALS
%   i          -     General purpose counter
%   jmax       -     Upper bound on j
%   imin       -     Lower bound on i during convolutions
%   imax       -     Upper bound on i during convolutions
%   eng        -     Energy
%   cor        -     Correlation
%   y2         -     Error weighting filtered reconstructed pitch prediction signal
%
% GLOBALS
%   h          -     Impulse response, perceptual weighting filter
%   e0         -     First error signal
%   Ypg        -     Convoluation: h * ex
%
% CONSTANTS
%   MAXLP      -     Maximum pitch prediction frame size
%
% ******************************************************************

function [ Pgain, match ] = pgain( ex, l, first, m, len )

% DECLARE GLOBAL VARIABLES
global e0 h Ypg

% DECLARE GLOBAL CONSTANTS
global MAXLP

% INITIALIZE LOCAL VARIABLES
y2 = zeros( MAXLP, 1 );

% RUN GAIN COMPUTATIONS
if first == 1

    % CALCULATE AND SAVE CONVOLUTION OF TRUNCATED (TO LEN)
    % IMPLUSE RESPONSE FOR FIRST LAG OF T (=MMIN) SAMPLES:
    %
    %           MIN(i, len-1)
    %      y     =  SUM  h * ex       , WHERE i = 0, ..., L-1 POINTS
    %       i, t    j=0   j    i-j
    %
    %                     h |0 1...len-1 x x|
    %      ex |L-1  . . .  1 0|               = y[0]
    %        ex |L-1  . . .  1 0|             = y[1]
    %                          :                :
    %                    ex |L-1  . . .  1 0| = y[L-1]
    for i = 0:l-1
	jmax = min( i, len-1 );
	Ypg(i+1) = sum( h(1:jmax+1) .* ex(i+1:-1:i-jmax+1) );
    end
else

    % END CORRECT THE CONVOLUTION SUM ON SUBSEQUENT PITCH LAGS
    %      y  =  0
    %   0, t
    %   y     =  y        + ex  * h   WHERE i = 1, ..., L POINTS
    %   i, m     i-1, m-1   -m    i  AND   m = t+1, ..., tmax LAGS
    Ypg(len-1:-1:1) = Ypg(len-1:-1:1) + ( ex(1) * h(len:-1:2) );
    Ypg(l:-1:2) = Ypg(l-1:-1:1);
    Ypg(1) = ex(1) * h(1);
end

% FOR LAGS (M) SHORTER THAN FRAME SIZE (L), REPLICATE THE SHORT
% ADAPTIVE CODEWORD TO THE FULL CODEWORD LENGTH BY OVERLAPPING
% AND ADDING THE CONVOLUTION
y2(1:l) = Ypg(1:l);
if m < l
    % ADD IN 2ND CONVOLUTION
    y2(m+1:l) = Ypg(m+1:l) + Ypg(1:l-m);
    if m < fix(l/2)
	% ADD IN 3RD CONVOLUTION
	imin = ( 2 * m ) + 1;
	imax = l;
	y2( imin:imax ) = y2( imin:imax ) + Ypg( 1:l-(2*m) );
    end
end

% CALCULATE CORRELATION AND ENERGY
% E0 = R(N)   = SPECTRUM PREDICTION RESIDUAL
% Y2 = R(N-M) = ERROR WEIGHTING FILTER RECONSTRUCTED PITCH PREDICTION
%               SIGNAL (M = CORRELATION LAG)
cor = sum( y2(1:l) .* e0(1:l) );
eng = sum( y2(1:l) .* y2(1:l) );

% COMPUTE GAIN AND ERROR
% ACTUAL MSPE = E0.E0 - PGAIN(2*COR-PGAIN*ENG) SINCE E0.E0 IS INDEPENDENT
% OF THE CODE WORD, MINIMIZING MSPE IS EQUIVALENT TO MAXIMIZING:
%
% MATCH = PGAIN( 2*COR - PGAIN*ENG )  (1)
%
% IF UNQUANTIZED PGAIN IS USED, THIS SIMPLIFIES:
%
% MATCH = COR * PGAIN  (2)
%
% NOTES IN THE CELP 3.2A SOURCE CODE FROM NSA INDICATE THAT INFERIOR
% RESULTS WERE OBTAINED WHEN QUANTIZED PGAIN WAS USED IN EQUATION (1).
% ALSO, NOTE THAT WHEN DELAY IS LESS THAN THE FRAME LENGTH, "MATCH" IS
% ONLY AN APPROXIMATION TO THE ACTUAL ERROR.
%
% INDEPENDENT (OPEN-LOOP) QUANTIZATION OF GAIN AND MATCH (INDEX):
if eng <= 0.0
    eng = 1.0;
end
Pgain = cor / eng;
match = cor * Pgain;