cfind.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.
%
% ******************************************************************
% CFIND
%
% PORTED TO MATLAB FROM CELP 3.2a C RELEASE
% 6-16-94
%
% ******************************************************************
%
% DESCRIPTION
%
% Compute filter coefficients, cepstral coefficients, and filter
% coefficient autocorrelation lags.
%
% DESIGN NOTES
%
% See references.
%
% REFERENCES
%
% 1. "Distance Measures for Speech Processing", A.H. Gray
%    and J.D. Markel, IEEE Trans. on ASSP, Vol. ASSP-24,
%    no. 5, Oct. 1976
%
% 2. "Quantization and Bit Allocation in Speech Processing",
%    A.H. Gray and J.D. Markel, IEEE Trans. on ASSP, Vol. ASSP-24
%    no. 6, Dec. 1976
%
% 3. "A Note on Quantization and Bit Allocation in Speech Processing",
%    A.H. Gray and J.D. Markel, IEEE Trans. on ASSP, Vol. ASSP-25
%    no. 3, June 1977
%
% VARIABLES
%
% INPUTS
%   m          -     Predictor order
%   nf         -     Number of cepstral terms to compute
%   r          -     Autocorrelation sequence
%
% OUTPUTS
%   cep        -     Cepstral coefficients
%   ra         -     Filter autocorrelation lags
%   alpha      -
%   a          -     Filter coefficients
%   rc         -     Reflection coefficients
%
% INTERNALS
%   i          -     Counter index
%   z          -     Loop upper bound
%   q          -     Intermediate result of filter coefficient comps
%   at         -         "         "          "            "
%
% CONSTANTS
%   MAXNO      -     Filter order
%
% ******************************************************************

function [ cep, ra, alpha, a, rc ] = cfind( m, nf, r )

% DECLARE GLOBAL CONSTANTS
global MAXNO

% INITIALIZE LOCAL VARIABLES AND RETURN VECTORS
a = zeros( (2*MAXNO)+1, 1 );
rc = zeros( (2*MAXNO)+1, 1 );
ra = zeros( m+1, 1 );
cep = zeros( nf, 1 );
a(1) = 1.0;
if r(1) == 0
    fprintf( 'find: r(0) = 0.0, resetting to 1e-6\n' );
    r(1) = 1.0e-6;
end

% OBTAIN FILTER COEFFICIENTS
rc(1) = -r(2) / r(1);
a(2) = rc(1);
alpha = r(1) * ( 1.0 - ( a(2) * a(2) ) );
for i = 2:m
    q = r(i+1) + sum( a(2:i) .* r(i:-1:2) );
    q = -q / alpha;
    rc(i) = q;
    z = fix( i/2 );
    at = a(2:z+1) + ( q * a(i:-1:i-z+1) );
    a(i:-1:i-z+1) = a(i:-1:i-z+1) + ( q * a(2:z+1) );
    a(2:z+1) = at;
    a(i+1) = q;
    alpha = alpha * ( 1 - ( q * q ) );

    % TRAP UNSTABLE FILTER
    if alpha <= 0.00
        fprintf( 'find: unstable filter\n' );
    end
end

% EVALUATE CEPSTRUM
cep(1) = a(2);
for i = 2:m
    cep(i) = i * a(i+1);
    cep(i) = cep(i) - sum( cep(1:i-1) .* a(i:-1:2) );
end
if nf > m
    for i = m:nf-1
        cep(i+1) = -sum( cep(i:-1:i-m+1) .* a(2:m+1) );
    end
    cep = -cep ./ (1:nf)';
end

% EVALUATE FILTER POLYNOMIAL AUTOCORRELATION
for i = 1:m+1
    ra(i) = sum( a(1:m-i+2) .* a(i:m+1) );
end