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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
.TH corr G "April 1993" "Scilab Group" "Scilab Function"
.so ../sci.an
.SH NAME
corr - correlation, covariance
.SH CALLING SEQUENCE
.nf
[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)
.fi
.SH PARAMETERS
.TP 10
x
: a real vector
.TP
y
: a real vector, default value x.
.TP
nlags
: integer, number of correlation coefficients desired.
.TP
xmacro
: a scilab external (see below).
.TP
ymacro
: a scilab external (see below), default value xmacro
.TP
n
: an integer, total size of the sequence (see below).
.TP
sect
: size of sections of the sequence (see below).
.TP
xi
: a real vector
.TP
yi
: a real vector,default value xi.
.TP
cov
: real vector, the correlation coefficients
.TP
Mean
: real number or vector, the mean of x and if given y
.SH DESCRIPTION
Computes
.nf
n - m
====
\\ 1
cov(m) = > (x(k) - xmean) (y(m+k) - ymean) * ---
/ n
====
k = 1
.fi
.LA $$ cov(m) = 1/n \sum_1^{n-m} (x(k) - E(x)) (y(m+k) - E(y)) $$
for m=0,..,\fVnlag-1\fR and two vectors \fVx=[x(1),..,x(n)]\fR
\fVy=[y(1),..,y(n)]\fR
.LP
Note that if x and y sequences are differents corr(x,y,...) is
different with corr(y,x,...)
.LP
Short sequences:
.LP
.Vb [cov,Mean]=corr(x,[y],nlags)
returns the first nlags
correlation coefficients and Mean = \fVmean(x)\fR
(mean of \fV[x,y]\fR if \fVy\fR is an argument).
The sequence \fVx\fR (resp. \fVy\fR) is assumed real, and \fVx\fR
and \fVy\fR are of same dimension n.
.LP
Long sequences:
.LP
.Vb [cov,Mean]=corr('fft',xmacro,[ymacro],n,sect)
.LP
Here \fVxmacro\fR is either
.TP
-
a function of type \fV[xx]=xmacro(sect,istart)\fR which returns
a vector \fVxx\fR of dimension \fVnsect\fR containing the part of the
sequence with indices from \fVistart\fR to \fVistart+sect-1\fR.
.TP
-
a fortran subroutine which performs the same calculation.
(See the source code of \fVdgetx\fR for an example).
\fVn\fR = total size of the sequence.
\fVsect\fR = size of sections of the sequence. \fVsect\fR must be a power
of 2. \fVcov\fR has dimension \fVsect\fR. Calculation is performed by FFT.
.LP
"Updating method":
.nf
[w,xu]=corr('updt',x1,[y1],w0)
[w,xu]=corr('updt',x2,[y2],w,xu)
...
wk=corr('updt',xk,[yk],w,xu)
.fi
With this calling sequence the calculation is updated at each
call to \fVcorr\fR.
.nf
w0 = 0*ones(1,2*nlags);
nlags = power of 2.
.fi
\fVx1,x2,...\fR are parts of \fVx\fR such that \fVx=[x1,x2,...]\fR and sizes
of \fVxi\fR a power of 2.
To get \fVnlags\fR coefficients a final fft must be performed
\fVc=fft(w,1)/n\fR; \fVcov=c(1nlags)\fR (\fVn\fR is the size of \fVx (y)\fR).
Caution: this calling sequence assumes that \fVxmean = ymean = 0\fR.
.fi
.SH EXAMPLE
.nf
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))]
.fi
.SH SEE ALSO
fft
|
