Home > voicebox > lpcstable.m

lpcstable

PURPOSE ^

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

SYNOPSIS ^

function [m,a]=lpcstable(ar)

DESCRIPTION ^

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

 Usage: (1) [m,ar]=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]=lpcstable(ar)
0002 %LPCSTABLE Test AR coefficients for stability and stabilize if necessary [MA,A]=(AR)
0003 %
0004 % Usage: (1) [m,ar]=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: lpcstable.m 8558 2016-09-22 08:22:54Z 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 General Public License as published by
0020 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0027 %
0028 %   You can obtain a copy of the GNU General Public License from
0029 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0030 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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

Generated on Fri 22-Sep-2017 19:37:38 by m2html © 2003