1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
|
// [ir,r,cl] = cra(z,M,n)
// M: No. of lags
// n: Order of pre whitening filter
// z is of the form:
// z = [y u]
// y and u are column matrices
function [ir,r,cl] = cra(z,varargin)
len = length(varargin);
if len==0,M = 20;n = 10;
elseif len==1,M = varargin(1);n=10;
else M = varargin(1);n = varargin(2);
end;
[ro,co] = size(z);
a1 = armax1(n,-1,0,z(:,2)',zeros(1,ro));
[A,B,D] = arma2p(a1);
a11 = coeff(A);
z22(1,:) = filt(a11,1,z(:,1)');
z22(2,:) = filt(a11,1,z(:,2)');
Covr = covf(z22',M+1);
scir=Covr(4,1);
sccf=sqrt(Covr(1,1)*Covr(4,1));
ir = Covr(2,:)'/scir;
r(:,1) = -M:M;
r(:,2) = [Covr(1,M+1:-1:1) Covr(1,2:M+1)]'
r(:,3) = [Covr(4,M+1:-1:1) Covr(4,2:M+1)]'
rhoyu = Covr(3,:)'/sccf;
rhouy = Covr(2,:)'/sccf;
r(:,4) = [rhoyu(M+1:-1:1); rhouy(2:M+1)];
sdreu=2.58*sqrt(r(:,3)'*r(:,2))/scir/sqrt(ro)*ones(2*M+1,1);
cl=sdreu(1);
timeax=[0:length(ir)-1];
plot(timeax,ir,'bo');
plot2d3(timeax,ir,style = 2);
plot(timeax,cl*ones(1,length(ir)),'b-.');
plot(timeax,-cl*ones(1,length(ir)),'b-.');
xtitle('Impulse response estimate','Lags');
xgrid();
endfunction;
|
