pctolsp2.m

实现fs1016w的CELP的低速率语音编解码功能的基于vc开发环境的原代码。

% 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.
%
% ******************************************************************
% PCTOLSP2
%
% PORTED TO MATLAB FROM CELP 3.2a C RELEASE
% 6-13-94
%
% ******************************************************************
%
% DESCRIPTION
%
% Compute LSP from predictor polynomial, A(Z)
%
% DESIGN NOTES
%
% Compute lsp frequencies by disection method as described in
% reference 1, below.
%
% LPC predictor coefficient convention is:
%
%          p+1       -(i-1)
%   A(z) = SUM   a  z          where a  = +1.0
%          i=1    i                   1
%
% Here, n=128, eps=1.e-06, nb=15 (this appears to be overkill)
%
% REFERENCES
%
% 1.  Line Spectrum Pair (LSP) and Speech Data Compression,
%     F.K. Soong and B-H Juang,
%     Proc. ICASSP 84, pp. 1.10.1-1.10.4
%
% VARIABLES
%
% INPUTS
%   a        -      LPC predictor polynomial
%   m        -      LPC predictor order
%
% OUTPUTS
%   freq     -      LSP frequencies
%   lspflag  -      Status of ill-conditioned lsp test
%
% INTERNALS
%   mp       -      m+1
%   mh       -      half m, or m/2
%   p        -      LSP conversion polynomial
%   q        -      LSP conversion polynomial
%   i        -      General purpose counter
%   j        -            "        "
%   mb       -            "        "
%   pxl      -      Intermediate conversion result
%   pxr      -            "        "        "
%   tpxr     -            "        "        "
%   pxm      -            "        "        "
%   qxm      -            "        "        "
%   qxl      -            "        "        "
%   tqxl     -            "        "        "
%   qxr      -            "        "        "
%   tqxr     -            "        "        "
%   fl       -            "        "        "
%   fm       -            "        "        "
%   fr       -            "        "        "
%   tfr      -            "        "        "
%   nf       -            "        "        "
%   jc       -            "        "        "
%   ang      -            "        "        "
%   tempfreq -      Temp storage for value swapping
%
% GLOBALS
%   lastfreq -      LSPs from previous frame
%
% CONSTANTS
%   EPS      -      Precision for computing polynomial zeros
%   NB       -      Iteration limit
%   MAXORD   -      Order of conversion polynomial
%   N        -      Maximum iterations to find P polynomial zeros
%   TRUE     -      Flag constant
%   FALSE    -      Flag constant
%   MAXNO    -      Predictor polynomial order
%
% ******************************************************************

function [ freq, lspflag ] = pctolsp2( a, m )

% DECLARE GLOBAL CONSTANTS
global TRUE FALSE MAXNO

% DECLARE GLOBAL VARIABLES (STATIC)
global lastfreq

% DEFINE CONSTANTS
EPS = 1.00e-6;
N = 128;
NB = 15;
MAXORD = 24;

% INIT LOCAL VARIABLES
p = zeros( MAXORD, 1 );
q = zeros( MAXORD, 1 );
freq = zeros( MAXNO+1, 1 );
mp = m + 1;
mh = fix ( m/2 );

% GENERATE P AND Q POLYNOMIALS
p( 1:mh ) = a( 2:mh+1 ) + a( m+1:-1:m-mh+2 );
q( 1:mh ) = a( 2:mh+1 ) - a( m+1:-1:m-mh+2 );

% COMPUTE P AT F=0
fl = 0.0;
pxl = 1.0 + sum( p(1:mh) );

% SEARCH FOR ZEROS OF P
nf = 1;
i = 1;
while i <= N
    mb = 0;
    fr = i * ( 0.5 / N );
    pxr = cos( mp * pi * fr );
    jc = mp - ( 2 * ( 1:mh )' );
    ang = pi * fr * jc;
    pxr = pxr + sum( cos(ang) .* p(1:mh) );
    tpxr = pxr;
    tfr = fr;
    if ( pxl * pxr ) <= 0.00
        mb = mb + 1;
        fm = fl + ( fr - fl ) / ( pxl - pxr ) * pxl;
        pxm = cos( mp * pi * fm );
        jc = mp - ( (1:mh)' * 2 );
        ang = pi * fm * jc;
        pxm = pxm + sum( cos(ang) .* p(1:mh) );
        if ( pxm * pxl ) > 0.00
            pxl = pxm;
            fl = fm;
        else
            pxr = pxm;
            fr = fm;
        end
        while (abs(pxm) > EPS) & ( mb < 4 )
            mb = mb + 1;
            fm = fl + ( fr - fl ) / ( pxl - pxr ) * pxl;
            pxm = cos( mp * pi * fm );
            jc = mp - ( (1:mh)' * 2 );
            ang = pi * fm * jc;
            pxm = pxm + sum( cos(ang) .* p(1:mh) );
            if ( pxm * pxl ) > 0.00
                pxl = pxm;
                fl = fm;
            else
                pxr = pxm;
                fr = fm;
            end
        end
        if ( (pxl-pxr) * pxl ) == 0
            freq( 1:m ) = ( 1:m ) * 0.04545;
            fprintf( 'pctolsp2: default lsps used, avoiding /0\n');
            return
        end
        freq(nf) = fl + (fr-fl) / (pxl-pxr) * pxl;
        nf = nf + 2;
        if nf > m-1
            break
        end
    end
    pxl = tpxr;
    fl = tfr;
    i = i + 1;
end

% SEARCH FOR ZEROS OF Q
freq(m+1) = 0.5;
fl = freq(1);
qxl = sin( pi * mp * fl );
jc = mp - ( (1:mh)' * 2 );
ang = pi * fl * jc;
qxl = qxl + sum( sin(ang) .* q(1:mh) );

i = 2;
while i < mp
    mb = 0;
    fr = freq(i+1);
    qxr = sin( mp * pi * fr );
    jc = mp - ( (1:mh)' * 2 );
    ang = pi * fr * jc;
    qxr = qxr + sum( sin(ang) .* q(1:mh) );
    tqxl = qxl;
    tfr = fr;
    tqxr = qxr;
    mb = mb + 1;
    fm = ( fl + fr ) * 0.5;
    qxm = sin( mp * pi * fm );
    jc = mp - ( (1:mh)' * 2 );
    ang = pi * fm * jc;
    qxm = qxm + sum( sin(ang) .* q(1:mh) );
    if ( qxm * qxl ) > 0.00
        qxl = qxm;
        fl = fm;
    else
        qxr = qxm;
        fr = fm;
    end
    while ( abs(qxm) > EPS*tqxl ) & ( mb < NB )
        mb = mb + 1;
        fm = ( fl + fr ) * 0.5;
        qxm = sin( mp * pi * fm );
        jc = mp - ( (1:mh)' * 2 );
        ang = pi * fm * jc;
        qxm = qxm + sum( sin(ang) .* q(1:mh) );
        if ( qxm * qxl ) > 0.00
            qxl = qxm;
            fl = fm;
        else
            qxr = qxm;
            fr = fm;
        end
    end
    if ( ( qxl - qxr ) * qxl ) == 0
        freq(1:m) = lastfreq(1:m);
        fprintf( 'pctolsp2: last lsps used, avoiding /0\n');
        return
    end
    freq(i) = fl + ( fr - fl ) / ( qxl - qxr ) * qxl;
    qxl = tqxr;
    fl = tfr;
    i = i + 2;
end

% TEST FOR ILL-CONDITIONED CASES AND TAKE CORRECTIVE ACTION IF REQUIRED
freq = freq(1:m);
lspflag = FALSE;
if ( any(freq==0.00) | any(freq==0.5) )
    lspflag = TRUE;
end

% REORDER LSPs IF NON-MONOTONIC
for i = 2:m
    if freq(i) < freq(i-1)
        lspflag = TRUE;
        fprintf( 'pctolsp2: non-monotonic lsps\n' );
        tempfreq = freq(i);
        freq(i) = freq(i-1);
        freq(i-1) = tempfreq;
    end
end

% IF NON-MONOTONIC AFTER 1ST PASS, RESET TO VALUES FROM PREVIOUS FRAME
for i = 2:m
    if freq(i) < freq(i-1)
        fprintf( 'pctolsp2: Reset to previous lsp values\n' );
        freq(1:m) = lastfreq;
        break;
    end
end

% SAVE CURRENT LSPS FOR NEXT FRAME ANALYSIS
lastfreq = freq(1:m);