V_LPCRA2AR Convert inverse filter autocorrelation coefs to AR filter. AR=(RA)
Usage: (1) ar0=poly([0.5 0.2]); % ar0=[1 -0.7 0.1]
ra0=v_lpcar2ra(ar0); % ra0=[1.5 -0.77 0.1]
ar1=v_lpcra2ar(ra0); % ar1 = ar0
(2) ar0=poly([0.5 0.2]); % ar0=[1 -0.7 0.1]
arx=xcorr(ar0); % arx=[0.1 -0.77 1.5 -0.77 0.1]
ar1=v_lpcra2ar(arx(length(ar0):end)); % ar1 = ar0
Inputs: ra(n,p+1) each row is the second half of the autocorrelation of
the coefficients of a stable AR filter of order p
tol tolerance relative to ra(1) [1e-8]
Outputs: ar(n,p+1) each row gives coefficients of an AR filter of order p
This routine uses a Newton-Raphson iteration described in [1] to invert
the cross-correlation function (as in the second usage example above).
Refs:
[1] G. Wilson. Factorization of the covariance generating function of a pure moving average process.
SIAM Journal on Numerical Analysis, 6 (1): 1�7, 1969.0001 function ar=v_lpcra2ar(ra,tol) 0002 %V_LPCRA2AR Convert inverse filter autocorrelation coefs to AR filter. AR=(RA) 0003 % 0004 % Usage: (1) ar0=poly([0.5 0.2]); % ar0=[1 -0.7 0.1] 0005 % ra0=v_lpcar2ra(ar0); % ra0=[1.5 -0.77 0.1] 0006 % ar1=v_lpcra2ar(ra0); % ar1 = ar0 0007 % 0008 % (2) ar0=poly([0.5 0.2]); % ar0=[1 -0.7 0.1] 0009 % arx=xcorr(ar0); % arx=[0.1 -0.77 1.5 -0.77 0.1] 0010 % ar1=v_lpcra2ar(arx(length(ar0):end)); % ar1 = ar0 0011 % 0012 % Inputs: ra(n,p+1) each row is the second half of the autocorrelation of 0013 % the coefficients of a stable AR filter of order p 0014 % tol tolerance relative to ra(1) [1e-8] 0015 % 0016 % Outputs: ar(n,p+1) each row gives coefficients of an AR filter of order p 0017 % 0018 % This routine uses a Newton-Raphson iteration described in [1] to invert 0019 % the cross-correlation function (as in the second usage example above). 0020 % 0021 % Refs: 0022 % [1] G. Wilson. Factorization of the covariance generating function of a pure moving average process. 0023 % SIAM Journal on Numerical Analysis, 6 (1): 1�7, 1969. 0024 0025 % Copyright (C) Mike Brookes 2015 0026 % Version: $Id: v_lpcra2ar.m 10865 2018-09-21 17:22:45Z dmb $ 0027 % 0028 % VOICEBOX is a MATLAB toolbox for speech processing. 0029 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0030 % 0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0032 % This program is free software; you can redistribute it and/or modify 0033 % it under the terms of the GNU General Public License as published by 0034 % the Free Software Foundation; either version 2 of the License, or 0035 % (at your option) any later version. 0036 % 0037 % This program is distributed in the hope that it will be useful, 0038 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0039 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0040 % GNU General Public License for more details. 0041 % 0042 % You can obtain a copy of the GNU General Public License from 0043 % http://www.gnu.org/copyleft/gpl.html or by writing to 0044 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 persistent p0 i1 i2 j1 j2 0047 if nargin<2 0048 tol=1e-8; % tolerance in ra/ra(1) 0049 end 0050 imax=20; 0051 [nf,pp]=size(ra); 0052 if ~numel(p0) || p0~=pp % create index lists for hankel and toeplitz matrices 0053 p0=pp; 0054 ix=zeros(1,(pp)*(pp+1)/2); 0055 nn=1:pp; 0056 ix(1+(nn-1).*nn/2)=1; 0057 j1=cumsum(ix); 0058 i1=cumsum(pp-1+(j1*(1-pp)+pp).*ix)-pp+1; 0059 j2=pp+1-j1; 0060 i2=cumsum(((pp+1)*j1-1).*ix-pp-1)+pp+1; 0061 end 0062 ar=zeros(nf,pp); % space for output filter coefficients 0063 t1=zeros(pp,pp); % space for hankel coefficient matrix 0064 t2=t1; % space for toeplitz lower triangular coefficient matrix 0065 ax0=zeros(1,pp); % temporary filter coefficient row vector 0066 for n=1:nf % process the input matrix one row at a time 0067 xa=ra(n,:); % pick out row n 0068 ax=ax0; % initialize ax to have all roots at zero 0069 ax(1)=sqrt(xa(1)+2*sum(xa(2:end))); 0070 i=imax; % maximum number of iterations 0071 while(i>0) 0072 t1(i1)=ax(j1); % t1=hankel(ax) 0073 t2(i2)=ax(j2); % t2=toeplitz(ax,[ax(1) zeros(1,p)]) 0074 ct=ax*t1; 0075 ax=(xa+ct)/(t1+t2); 0076 i=min(i-1,i*(max(abs(ct-xa))>tol*xa(1))+1); % do one final iteration after tolerance reached 0077 end 0078 ar(n,:)=ax; 0079 end