# v_lpcstable

## PURPOSE

V_LPCSTABLE Test AR coefficients for stability and stabilize if necessary [MA,A]=(AR)

## SYNOPSIS

function [m,a]=v_lpcstable(ar)

## DESCRIPTION

    V_LPCSTABLE Test AR coefficients for stability and stabilize if necessary [MA,A]=(AR)

     Usage: (1) [m,ar]=v_lpcstable(ar); % force ar polynolials to be stable

       Input:  ar(:,p+1)  Autoregressive coefficients
     Outputs:  m(:,1)    Mask identifying stable polynomials
               a(:,p+1)  Stabilized polynomials formed by reflecting unstable
                           poles in unit circle (with a(:,1)=1)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function [m,a]=v_lpcstable(ar)
    0002 %V_LPCSTABLE Test AR coefficients for stability and stabilize if necessary [MA,A]=(AR)
    0003 %
    0004 % Usage: (1) [m,ar]=v_lpcstable(ar); % force ar polynolials to be stable
    0005 %
    0006 %   Input:  ar(:,p+1)  Autoregressive coefficients
    0007 % Outputs:  m(:,1)    Mask identifying stable polynomials
    0008 %           a(:,p+1)  Stabilized polynomials formed by reflecting unstable
    0009 %                       poles in unit circle (with a(:,1)=1)
    0010
    0011 %      Copyright (C) Mike Brookes 2016
    0012 %      Version: $Id: v_lpcstable.m 10865 2018-09-21 17:22:45Z dmb $
    0013 %
    0014 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0015 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0016 %
    0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0018 %   This program is free software; you can redistribute it and/or modify
    0019 %   it under the terms of the GNU Lesser General Public License as published by
    0020 %   the Free Software Foundation; either version 3 of the License, or
    0021 %   (at your option) any later version.
    0022 %
    0023 %   This program is distributed in the hope that it will be useful,
    0024 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0025 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0026 %   GNU Lesser General Public License for more details.
    0027 %
    0028 %   You can obtain a copy of the GNU Lesser General Public License from
    0029 %   https://www.gnu.org/licenses/ .
    0030 %
    0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0032
    0033 [nf,p1] = size(ar);
    0034 mm=ar(:,1)~=1;
    0035 if any(mm)          % first ensure leading coefficient is always 1
    0036     ar(mm,:)=ar(mm,:)./ar(mm,ones(1,p1));
    0037 end
    0038 if p1==1
    0039     m=ones(nf,1); % 0'th order filter is always stable
    0040 elseif p1==2
    0041     m=abs(ar(:,2))<1; % 1'st order filter
    0042 else
    0043     rf = ar;
    0044     k = rf(:,p1); % check final coefficient in range
    0045     m=abs(k)<1;
    0046     if any(m)
    0047         d = (1-k(m).^2).^(-1);
    0048         wj=ones(1,p1-2);
    0049         rf(m,2:p1-1) = (rf(m,2:p1-1)-k(m,wj).*rf(m,p1-1:-1:2)).*d(:,wj);
    0050         for j = p1-2:-1:2
    0051             k(m) = rf(m,j+1);
    0052             m(m)=abs(k)<1;
    0053             if ~any(m), break, end
    0054             d = (1-k(m).^2).^(-1);
    0055             wj=ones(1,j-1);
    0056             rf(m,2:j) = (rf(m,2:j)-k(m,wj).*rf(m,j:-1:2)).*d(:,wj);
    0057         end
    0058     end
    0059 end
    0060 if nargout>1
    0061     a=ar;
    0062     if ~all(m)
    0063         for i=find(~m)'                 % unstable frames
    0064             z=roots(a(i,:));
    0065             k=abs(z)>1;                 % find any unstable roots
    0066             z(k)=conj(z(k)).^(-1);      % invert them
    0067             a(i,:)=real(poly(z));       % force a real polynomial
    0068         end
    0069     end
    0070 end