function [ir,r,cl]=cra(z,M,n)
%CRA	Performs correlation analysis
%
%	IR = cra(Z)
%
%	Z: The data, entered with two columns Z = [y u]
% 	IR: The estimated impulse response (IR(1) corresponds to g(0))
%
%	[IR,R,CL] = cra(Z,M,NA) gives access to
%	M: The number of lags for which the functions are computed (def 20)
%	NA: The order of the whitening filter.(Def 10). With NA=0, no prewhite-
%	    ning is performed. Then the covariance functions of the original
%	    data are obtained.
%	R: The covariance/correlation information
%	   R(:,1) contains the lag indices
%	   R(:,2) contains the covariance function of y (poss. prewhitened)
%	   R(:,3) contains the covariance function of u (poss. prewhitened)
%	   R(:,4) contains the correlation function between (poss prewhitened)
%	     u and y (positive lags correponds to an influence from u to y)
%	   CL is the 99 % significance level for the impulse response
%	   These functions are also plotted by the cra. The plots can be 
%	   redisplayed by cra(R);
%


%	L. Ljung 10-2-90,1-25-92
%	Copyright (c) 1990-92 by the MathWorks, Inc.
%	All Rights Reserved.

[Ncap,nz]=size(z);redisp=0;
if nz==4, if z((Ncap+1)/2,1)==0, redisp=1;end,end
if nz>2 & ~redisp,error('This routine is for two data records only. To check correlations between more signals, use cra for two signals at a time.'),end
if redisp, r=z; n=0;end
clf reset
if ~redisp
if nargin<3, n=[];end
if nargin<2, M=[];end
if isempty(n),n=10;end
if isempty(M),M=20;end
if nz==1, error('For a single signal, use covf instead!'),end
if n>0
	a=th2poly(ar(z(:,2),n,'fb0'));
	z(:,1)=filter(a,1,z(:,1));
	z(:,2)=filter(a,1,z(:,2));
end
M=M-1;
Rcap=covf(z,M+1);
r(:,1)=[-M:1:M]';
r(M+1:2*M+1,2:3)=Rcap([1 4],:)';
r(1:M,2:3)=Rcap([1 4],M+1:-1:2)';
scir=Rcap(4,1); sccf=sqrt(Rcap(1,1)*Rcap(4,1));
r(M+1:2*M+1,4)=Rcap(2,:)'/sccf;
r(1:M,4)=Rcap(3,M+1:-1:2)'/sccf;
ir=r(M+1:2*M+1,4)*sccf/scir;
end
subplot(221)
plot(r(:,1),r(:,2)),
if n>0,title('Covf for filtered y'),else title('Covf for y'),end
subplot(222)
plot(r(:,1),r(:,3)),
if n>0,title('Covf for prewhitened u'),else title('Covf for u'),end
subplot(223)
plot(r(:,1),r(:,4))
if n>0,title('Correlation from u to y (prewh)'),else title('Correlation from u to y'),end
if ~redisp,
	sdreu=2.58*sqrt(r(:,3)'*r(:,2))/scir/sqrt(Ncap)*ones(2*M+1,1);
	subplot(224),plot(r(:,1),r(:,4)*sccf/scir)
	hold on,plot(r(:,1),sdreu,'g-.',r(:,1),-sdreu,'g-.'),hold off
	cl=sdreu(1);
	title('Impulse response estimate')
end
set(gcf,'NextPlot','replace');

