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)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 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