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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
corr Scilab Group Scilab Function corr
NAME
corr - correlation, covariance
CALLING SEQUENCE
[cov,Mean]=corr(x,[y],nlags)
[cov,Mean]=corr('fft',xmacro,[ymacro],n,sect)
[w,xu]=corr('updt',x1,[y1],w0)
[w,xu]=corr('updt',x2,[y2],w,xu)
...
[wk]=corr('updt',xk,[yk],w,xu)
PARAMETERS
x : a real vector
y : a real vector, default value x.
nlags : integer, number of correlation coefficients desired.
xmacro : a scilab external (see below).
ymacro : a scilab external (see below), default value xmacro
n : an integer, total size of the sequence (see below).
sect : size of sections of the sequence (see below).
xi : a real vector
yi : a real vector,default value xi.
cov : real vector, the correlation coefficients
Mean : real number or vector, the mean of x and if given y
DESCRIPTION
Computes
n - m
====
\ 1
cov(m) = > (x(k) - xmean) (y(m+k) - ymean) * ---
/ n
====
k = 1
for m=0,..,nlag-1 and two vectors x=[x(1),..,x(n)] y=[y(1),..,y(n)]
Note that if x and y sequences are differents corr(x,y,...) is different
with corr(y,x,...)
Short sequences:
returns the first nlags correlation coefficients and Mean = mean(x) (mean
of [x,y] if y is an argument). The sequence x (resp. y) is assumed real,
and x and y are of same dimension n.
Long sequences:
Here xmacro is either
- a function of type [xx]=xmacro(sect,istart) which returns a vector
xx of dimension nsect containing the part of the sequence with
indices from istart to istart+sect-1.
- a fortran subroutine which performs the same calculation. (See the
source code of dgetx for an example). n = total size of the sequence.
sect = size of sections of the sequence. sect must be a power of 2.
cov has dimension sect. Calculation is performed by FFT.
"Updating method":
[w,xu]=corr('updt',x1,[y1],w0)
[w,xu]=corr('updt',x2,[y2],w,xu)
...
wk=corr('updt',xk,[yk],w,xu)
With this calling sequence the calculation is updated at each call to
corr.
w0 = 0*ones(1,2*nlags);
nlags = power of 2.
x1,x2,... are parts of x such that x=[x1,x2,...] and sizes of xi a
power of 2. To get nlags coefficients a final fft must be performed
c=fft(w,1)/n; cov=c(1nlags) (n is the size of x (y)). Caution: this
calling sequence assumes that xmean = ymean = 0.
EXAMPLE
x=%pi/10:%pi/10:102.4*%pi;
rand('seed');rand('normal');
y=[.8*sin(x)+.8*sin(2*x)+rand(x);.8*sin(x)+.8*sin(1.99*x)+rand(x)];
c=[];
for j=1:2,for k=1:2,c=[c;corr(y(k,:),y(j,:),64)];end;end;
c=matrix(c,2,128);cov=[];
for j=1:64,cov=[cov;c(:,(j-1)*2+1:2*j)];end;
rand('unif')
//
rand('normal');x=rand(1,256);y=-x;
deff('[z]=xx(inc,is)','z=x(is:is+inc-1)');
deff('[z]=yy(inc,is)','z=y(is:is+inc-1)');
[c,mxy]=corr(x,y,32);
x=x-mxy(1)*ones(x);y=y-mxy(2)*ones(y); //centring
c1=corr(x,y,32);c2=corr(x,32);
norm(c1+c2,1)
[c3,m3]=corr('fft',xx,yy,256,32);
norm(c1-c3,1)
[c4,m4]=corr('fft',xx,256,32);
norm(m3,1),norm(m4,1)
norm(c3-c1,1),norm(c4-c2,1)
x1=x(1:128);x2=x(129:256);
y1=y(1:128);y2=y(129:256);
w0=0*ones(1:64); //32 coeffs
[w1,xu]=corr('u',x1,y1,w0);w2=corr('u',x2,y2,w1,xu);
zz=real(fft(w2,1))/256;c5=zz(1:32);
norm(c5-c1,1)
[w1,xu]=corr('u',x1,w0);w2=corr('u',x2,w1,xu);
zz=real(fft(w2,1))/256;c6=zz(1:32);
norm(c6-c2,1)
rand('unif')
// test for Fortran or C external
//
deff('[y]=xmacro(sec,ist)','y=sin(ist:(ist+sec-1))');
x=xmacro(100,1);
[cc1,mm1]=corr(x,2^3);
[cc,mm]=corr('fft',xmacro,100,2^3);
[cc2,mm2]=corr('fft','corexx',100,2^3);
[maxi(abs(cc-cc1)),maxi(abs(mm-mm1)),maxi(abs(cc-cc2)),maxi(abs(mm-mm2))]
deff('[y]=ymacro(sec,ist)','y=cos(ist:(ist+sec-1))');
y=ymacro(100,1);
[cc1,mm1]=corr(x,y,2^3);
[cc,mm]=corr('fft',xmacro,ymacro,100,2^3);
[cc2,mm2]=corr('fft','corexx','corexy',100,2^3);
[maxi(abs(cc-cc1)),maxi(abs(mm-mm1)),maxi(abs(cc-cc2)),maxi(abs(mm-mm2))]
SEE ALSO
fft
|
