quikarma.m

biomedical stuffs...

function [arp,map, s2w]=quikarma(yn,p,q)
%	function [arp,map, s2w]=quikarma(yn,p,q)
%	Can be used for MA or ARMA models.
%	For AR models, this routine does NOT use the standard Yule-Walker
%  equations - it uses the higher-order ones - and should not be used.
%
%	This program is distributed as a supplement to the book
%	"Biomedical Signal Processing and Signal Modeling" by E. N. Bruce,
%	published by Wiley, 2000.  It is provided for educational use only.
%  While every effort has been made to insure its suitability to illustrate
%  principles described in the above book, no specific feature or capability 
%  is implied or guaranteed.
%
%Solve higher-order Yule-Walker equations for q+1 <= m <= q+2p
xn=detrend(yn);lnx=length(xn);
rxn=xcov(detrend(xn),'biased');AA=zeros(2*p,p);
arp=[];rn=xn;
if p ~= 0
for j=1:p
	AA(:,j)=rxn(lnx+1-j+q:lnx+1-j+q+2*p-1);
end
AA
arp=tls(AA,-rxn(lnx+q+1:lnx+q+2*p));
s2w=sum(arp.*rxn(lnx+1:lnx+p))+rxn(lnx);
%  Filter input by inverse AR filter to generate MA sequence
rn=filter([1;arp]',1,xn);
end %if p ~= 0
subplot(211),stairs(xn)
subplot(212),stairs(rn)
% Find MA coefficients by Durbin's method
if q ~= 0
ar4ma=ar(detrend(rn),4*q);
AB=zeros(4*q,q);
[arnum,arden]=th2tf(ar4ma);
for j=1:q
	AB(:,j)=[zeros(j-1,1);arden(1:4*q-j+1)'];
end
AB,arnum,arden
map=tls(AB,-arden(2:4*q+1)');
if p ~= 0
hn=dimpulse([1;map]',[1;arp]');
else, hn=[1;map]';
end %if p == 0
%hn(1:10)
cq0=hn(1)
for j=1:q,j,cq0=cq0+map(j)*hn(j+1);end
cq0
s2w=(sum(arp.*rxn(lnx+1:lnx+p))+rxn(lnx))/cq0;
end %if q ~= 0