arp2pacf.m

时间序列工具箱

function [PACF,ARPmx]= arp2pacf(ARP);
%	ARP2PACF conversion or Inverse Recursion
%       Calculates the Refelction coefficients from the AR parameters
% 	[PACF]= arp2pacf(ARP);
%	converts ARP to PACF
%
% Input:	ARP  AR-Parametermatrix
% Output:	PACF partitial autocorrelation function 
% 		(According to Haykin 1986 is it the PARCOR coeefficient or refection coefficient GAMMA)

%	Version 1.01
%	last revision 30.03.1998
%	Copyright (c) 1998 by Alois Schloegl
%	e-mail: [e-mail]	

%This library is free software; you can redistribute it and/or
%modify it under the terms of the GNU Library General Public
%License as published by the Free Software Foundation; either
%version 2 of the License, or (at your option) any later version.
%
%This library is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%Library General Public License for more details.
%
%You should have received a copy of the GNU Library General Public
%License along with this library; if not, write to the
%Free Software Foundation, Inc., 59 Temple Place - Suite 330,
%Boston, MA 02111-1307, USA.

% Check Input 
[nr,nc]=size(ARP);
if (nr~=1) & (nc~=1)
           error('ARP must be vector');
end;
if nr~=1, ARP=ARP.'; nc=nr;
end;

% Start
M=nc;
%A% K=1; %A% 
    K=M; %B% 
ARPmx=zeros(K,M);
ARPmx(K,1:M)=[ARP];
PACF(M)=ARP(M);

%M% PACF=zeros(size(ARP));
%M% arp=PACF;
%M% PACF(:,M)=ARP(:,M);
%M% PE(:,M)=ARP(:,M);

% Inverse Recursion
for m=M:-1:2,
        %A% k=1;K=1;    %A% 
            k=m-1,K=m;  %B% 
        ARPmx(k,1:m-1)=(ARPmx(K,1:m-1)+ARPmx(K,m)*conj(ARPmx(K,m-1:-1:1)))/(1-ARPmx(K,m)*ARPmx(K,m)');
        PACF(m-1)=ARPmx(k,m-1);  
        
%M%         arp(:,1:m-1)=arp(:,1:m-1) + conj(arp(:,m-1:-1:1) *. (arp(:,m) ./ (1-arp(:,m).^2))*ones(lr,m-1);
%M%         PACF(:,m-1)=arp(:,m-1);
%           PE(:,m-1)=1-arp(:,m-1).*arp(:,m-1).';
end;
%PACF=diag(ARPmx);

return;

%PE(:,1)=ones(nr,1);
%        r(:,1)=PACF(:,m);
%for m=2:M
%        r(:,m)=-PE(:,m-1)*PACF(:,m)-sum(arp,2)  ;
%        PE(:,m)=PE(:,m-1).*(1-arp(:,m-1).*arp(:,m-1).');
%end;