psearch.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.
%
% ******************************************************************
% PSEARCH
%
% PORTED TO MATLAB FROM CELP 3.2a C RELEASE
% 6-30-94
%
% ******************************************************************
%
% DESCRIPTION
%
% Find pitch VQ parameters
%
% DESIGN NOTES
%
% Pitch search is performed by closed-loop analysis using a modification
% of what is commonly called any one of the following: "self-excited",
% "adaptive code book" or "VQ" method.  This method was found to be
% superior to the conventional "filtering approach", especially for high
% pitched speakers.  The filtering and VQ methods are identical except when
% the delay is less than the frame length.  The pitch delay ranges from minptr
% to maxptr (i.e., 20 to 147 including noninteger) lags every odd subframe
% while even subframes are searched and coded within 2**pbits[1] (i.e., 32)
% lags relative to the previous subframe.  The delta search greatly reduces
% the computational complexity and data rate while causing no percievable
% loss in speech quality.
%
% The minimum squared prediction error (MSPE) search criteria is modified
% to check the error at submultiples of the delay to determine if it is
% within 1 dB of the MSPE.  The submultiple delay is selected if its error
% satisfies our modified criteria.  This results in a smooth "pitch" delay
% contour which produces higher quality speech and is crucial to the
% synthesizer's smoother in the presence of bit errors.
%
% 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.)
%
% 3. Kroon, Peter and Bishnu Atal, "On Improving the Performance of
%    Pitch Predictions in Speech Coding Systems," IEEE Speech Coding
%    Workshop, September 1989.
%
% 4. Marques, J.S., et al., "Pitch Prediction with Fractional Delays
%    in CELP Coding," European Conference on Speech Communication and
%    Technology, September, 1989.
%
% VARIABLES
%
% INPUTS
%   l          -     Length of error signal
%
% INTERNALS
%   i          -     General purpose loop variable
%   m          -     Integer pitch delay
%   lag        -     Pitch lag
%   start      -     Search starting point
%   first      -     First call flag, psearch
%   bigptr     -     Index of largest neighborhood match score
%   initbp     -     Initial value of bigptr
%   hiscore    -     Largest match score found
%   smin       -     Minimum submult index for neighborhood search
%   smax       -     Maximum submult index for neighborhood search
%   topptr     -     Index of largest match score
%   maxptr     -     Maximum pitch delay pointer
%   bufptr     -     Pointer to end of buffer
%   g          -     Vector of gains for different delays
%   match      -     Vector of match scores for different delays
%   emax       -     Maximum match score
%   fraction   -     Fractional pitch delays flag
%   neigh      -     Neighboring delays flag
%   whole      -     Integer pitch delays flag
%   sub        -     Search subframes flag
%   nrange     -     Neighborhood search range
%   v0         -     Excitation vector
%   v0shift    -     Delayed excitation vector
%   frac       -     Fractional pitch delay
%   imin       -     Minimum i value during neighborhood search
%   imax       -     Maximum i value during neighborhood search
%
% GLOBALS
%   idb        -     Dimension of d1a, d1b
%   no         -     Predictor order
%   nseg       -     Segment (subframe) counter, cumulative
%   pindex     -     Pitch gain index
%   tauptr     -     Pitch delay pointer
%   minptr     -     Minimum delay pointer
%   oldptr     -     Previous subframe pitch pointer for delta calcs
%   plevel1    -     Number of full search pitch delays
%   plevel2    -     Number of delta search pitch delays
%   pdelay     -     Adaptive codebook pitch delays
%   pstype     -     Pitch search type (hier, full, or intg)
%   submult    -     Pitch submultiple delay table
%   ptype      -     Pitch codebook index quantization type
%   d1b        -     Filter memory, 1/P(Z)
%   bb         -     Pitch lag search parameters, gain and index
%   h          -     Impulse response, perceptual weighting filter
%
% CONSTANTS
%   LEN        -     Length of truncated impulse response
%   MAXBUFPTR  -     Length of search vector
%   MAXL       -     Maximum codeword vector length
%   MMAX       -     Maximum delay pitch predictor
%   MAXLP      -     Maximum pitch prediction frame size
%   MAXPD      -     Maximum number of pitch delays
%   MAXNO      -     Maximum LPC filter order
%   MAXNP      -     Maximum pitch prediction order
%   MMIN       -     Minimum delay pitch predictor
%   TRUE       -     General purpose boolean flag
%   FALSE      -     General purpose boolean flag
%
% ******************************************************************

function psearch( l )

% DECLARE GLOBAL VARIABLES
global idb no nseg pindex tauptr minptr oldptr submult
global plevel1 plevel2 bb d1b h pdelay ptype pstype

% DECLARE GLOBAL CONSTANTS
global MAXLP MAXPD MMAX MAXNO MAXNP MAXL MMIN TRUE FALSE

% INITIALIZE LOCAL CONSTANTS
MAXBUFPTR = MMAX + MAXNO + 2*MAXLP + MAXNP - 1;
LEN = 30;

% INITIALIZE LOCAL VARIABLES
v0 =  zeros( MAXBUFPTR, 1 );
v0shift = zeros( MAXLP, 1 );
g = zeros( MAXPD, 1 );
match = g;
nrange = 0;
bufptr = MMAX + no + 2*l + MAXNP - 1;

% CHOOSE TYPE OF PITCH DELAY SEARCH:
% 1. TWO STAGE HIERARCHICAL SEARCH OF INTEGER AND
%    NEIGHBORING NONINTEGER DELAYS
% 2. INTEGER ONLY SEARCH
% 3. FULL EXHAUSTIVE SEARCH
if strcmp( pstype, 'hier' ) == 1
    whole = 1;
    fraction = 0;
    sub = 1;
    neigh = 1;
    nrange = 3;
elseif strcmp( pstype, 'intg' ) == 1
    whole = 1;
    fraction = 0;
    sub = 1;
    neigh = 0;
elseif strcmp( pstype, 'full' ) == 1
    whole = 1;
    fraction = 1;
    sub = 1;
    neigh = 0;
else
    error( 'psearch: incorrect pitch search type (pstype)\n' );
    quit;
end

% CHECK ILL-CONDITIONED CASES
if LEN  > l
    error( 'psearch:  impulse response too long\n' );
end
if MAXLP < MAXL
    error( 'psearch:  MAXLP < MAXL \n' );
end

% UPDATE ADAPTIVE CODEBOOK (PITCH MEMORY)
v0low = bufptr - idb - l + 1;
v0high = v0low + idb - 1;
v0( v0low:v0high ) = d1b( 1:idb );

% RUN CODEBOOK SEARCH
if nseg == 1
    % SET INITIAL CONDITIONS
    bb(3) = 0.00;
    bb(1) = MMIN;
else
    % FIND ALLOWABLE POINTER RANGE (MINPTR TO MAXPTR)
    if rem( nseg, 2 ) == 0
        % USE DELTA DELAY CODING ON EVEN SUBFRAMES
        minptr = oldptr - ( fix( plevel2/2 ) - 1 );
        maxptr = oldptr + fix( plevel2/2 );
        if minptr < 0
            minptr = 0;
            maxptr = plevel2 - 1;
        end
        if maxptr > ( plevel1 - 1 )
            maxptr = plevel1 - 1;
            minptr = plevel1 - plevel2;
        end
    else
        % USE FULL RANGE CODING ON ODD SUBFRAMES
        minptr = 0;
        maxptr = plevel1 - 1;
    end
    start = bufptr - l + 1;

    % FIND GAIN AND MATCH SCORE FOR INTEGER PITCH DELAYS, USING
    % END-POINT CORRECTION ON UNITY SPACED DELAYS
    if whole == 1
        first = TRUE;
        for i = minptr:maxptr
            m = fix( pdelay(i+1) );
            frac = pdelay(i+1) - m;
            if abs(frac) < 1.0e-4
                lag = start - m;
                [ g(i+1), match(i+1) ] = pgain( v0(lag:lag+l-1), l, first, m, LEN );
                first = FALSE;
            else
                match( i+1 ) = 0.00;
            end
        end
    end

    % FIND GAIN AND MATCH SCORE FOR FRACTIONAL DELAYS
    if fraction == 1
        for i = minptr:maxptr
            m = fix( pdelay(i+1) );
            frac = pdelay(i+1) - m;
            if abs(frac) >= 1.0e-4
                [ v0, v0shift ] = delay( v0, start, l, frac, m );
                [ g(i+1), match(i+1) ] = pgain( v0shift, l, TRUE, 70, LEN );
            end
        end
    end

    % FIND POINTER TO TOP MATCH SCORE (MSPE)
    % SEARCH FOR BEST MATCH SCORE (MAX -ERROR TERM)
    [ emax, topptr ] = max( match(minptr+1:maxptr+1) );
    topptr = topptr + minptr - 1;

    % FOR FULL SEARCH (ODD) SUBFRAMES:
    % SELECT SHORTEST DELAY OF 2, 3, OR 4 SUBMULTIPLES IF ITS MATCH
    % IS WITHIN 1 DB OF MSPE TO FAVOR SMOOTH PITCH
    tauptr = topptr;
    if sub == 1
        if rem( nseg, 2 ) ~= 0

            % REPEAT FOR EACH SUBMULTIPLE: 2, 3 AND 4
            for i = 1:submult( topptr+1, 1 )

                % FIND BEST NEIGHBORHOOD MATCH FOR GIVEN SUBMULTIPLE
                bigptr =  submult( topptr+1, i+1 );
                initbp = bigptr;
                smin = max( bigptr - 8, minptr );
                smax = min( bigptr + 8, maxptr );
                [ hiscore, bigptr ] = max( match(smin+1:smax+1) );
                bigptr = bigptr + smin - 1;
                if match(initbp+1) >= hiscore
                    bigptr = initbp;
                end

                % SELECT SUBMULTIPLE MATCH IF WITHIN 1 DB OF MSPE MATCH
                if match( bigptr+1 ) >= ( 0.88 * match( topptr+1 ) )
                    tauptr = bigptr;
                end
            end
        end
    end

    % SEARCH TAUPTR'S NEIGHBORING DELAYS (TO BE USED WITH EARLIER STAGES
    % OF SEARCHING). FIND GAIN AND MATCH SCORE FOR NEIGHBORING DELAYS
    % AND FIND BEST NEIGHBORHOOD MATCH (COULD USE END-POINT CORRECTION
    % ON UNITY SPACED DELAYS
    if neigh == 1
        bigptr = tauptr;
        imin = max( tauptr-nrange, minptr );
        imax = min( tauptr+nrange, maxptr );
        for i = imin:imax
            if i ~= tauptr
                m = fix( pdelay( i+1 ) );
                frac = pdelay( i+1 ) - m;
                lag = start - m;
                if abs(frac) < 1.0e-4
                    [ g(i+1), match(i+1) ] = pgain( v0(lag:lag+l-1), l, TRUE, m, LEN );
                else
                    [ v0, v0shift ] = delay( v0, start, l, frac, m );
                    [ g(i+1), match(i+1) ] = pgain( v0shift, l, TRUE, 70, LEN );
                end
                if match(i+1) > match(tauptr+1)
                    bigptr = i;
                end
            end
        end
        tauptr = bigptr;
    end

    % GIVEN A POINTER TO DELAY (TAUPTR), RECOMPUTE ITS GAIN TO CORRECT
    % ERRORS ACCUMULATED IN RECURSIONS AND ERRORS DUE TO TRUNCATION
    m = fix( pdelay( tauptr+1 ) );
    frac = pdelay( tauptr+1 ) - m;
    lag = start - m;

    % TWO CASES: INTEGER DELAY OR FRACTIONAL DELAY
    if abs(frac) < 1.0e-4

        % INTEGER
        [ g(tauptr+1), match(tauptr+1) ] = pgain( v0(lag:lag+l-1), l, TRUE, m, l );
    else

        % FRACTIONAL
        [ v0, v0shift ] = delay( v0, start, l, frac, m );
        [ g(tauptr+1), match(tauptr+1) ] = pgain( v0shift, l, TRUE, 70, l );
    end

    % PLACE PITCH PARAMETERS IN COMMON BB STRUCTURE
    bb(3) = g( tauptr+1 );
    bb(1) = pdelay( tauptr+1 );

    % SAVE PITCH POINTER TO DETERMINE DELTA DELAY
    oldptr = tauptr;
end

% QUANTIZE PITCH, BB(3)
if strcmp( ptype, 'none' ) ~= 1
    [ pindex, bb(3) ] = pencode( bb(3) );
else
    error( 'psearch: no pitch quantization!' );
end