-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathCUM4X.M
169 lines (135 loc) · 5.59 KB
/
CUM4X.M
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
166
167
168
169
function y_cum = cum4x (w,x,y,z, maxlag, nsamp, overlap, flag, k1, k2)
%CUM4X Fourth-order cross-cumulants.
% y_cum = cum4x (w,x,y,z, maxlag, samp_seg, overlap, flag, k1, k2)
% w,x,y,z - data vectors/matrices with identical dimensions
% if w,x,y,z are matrices, rather than vectors, columns are
% assumed to correspond to independent realizations,
% overlap is set to 0, and samp_seg to the row dimension.
% maxlag - maximum lag to be computed [default = 0]
% samp_seg - samples per segment [default = data_length]
% overlap - percentage overlap of segments [default = 0]
% overlap is clipped to the allowed range of [0,99].
% flag : 'biased', biased estimates are computed [default]
% 'unbiased', unbiased estimates are computed.
% k1,k2 : the fixed lags in C4(m,k1,k2); defaults to 0
% y_cum: estimated fourth-order cross cumulant,
% c4(t1,t2,t3) := cum( w^*(t), x(t+t1), y(t+t2), z^*(t+t3) )
% Copyright (c) 1991-2001 by United Signals & Systems, Inc.
% $Revision: 1.4 $
% A. Swami January 20, 1995
% RESTRICTED RIGHTS LEGEND
% Use, duplication, or disclosure by the Government is subject to
% restrictions as set forth in subparagraph (c) (1) (ii) of the
% Rights in Technical Data and Computer Software clause of DFARS
% 252.227-7013.
% Manufacturer: United Signals & Systems, Inc., P.O. Box 2374,
% Culver City, California 90231.
%
% This material may be reproduced by or for the U.S. Government pursuant
% to the copyright license under the clause at DFARS 252.227-7013.
% c4(t1,t2,t3) := cum( x^*(t), w(t+t1), y(t+t2), z^*(t+t3) )
% cum(w,x,y,z) := E(wxyz) - E(wx)E(yz) - E(wy)E(xz) - E(wz)E(xy)
% and, w,x,y,z are assumed to be zero-mean.
% ---- Parameter checks are done in CUMEST ----------------------
[lx,nrecs] = size(w);
if ([lx,nrecs] ~= size(x) || [lx,nrecs] ~= size(y) || [lx,nrecs] ~= size(z))
error('w,x,y,z should have identical dimensions')
end
if (lx == 1)
lx = nrecs; nrecs = 1;
end
if (exist('maxlag') ~= 1) maxlag = 0; end
if (maxlag < 0)
error ('"maxlag" must be non-negative ')
end
if (exist('nsamp') ~= 1) nsamp = lx; end
if (nrecs > 1) nsamp = lx; end
if (nsamp <= 0 || nsamp > lx) nsamp = lx; end
if (exist('overlap') ~= 1) overlap = 0; end
if (nrecs > 1) overlap = 0; end
overlap = max(0,min(overlap,99));
if (exist('flag') ~= 1) flag = 'biased' ; end
if (flag(1:1) ~= 'b' && flag(1:1) ~= 'B')
flag = 'unbiased';
else, flag = 'biased';
end
overlap0 = overlap;
overlap = fix(overlap/100 * nsamp);
nadvance = nsamp - overlap;
if (nrecs == 1)
nrecs = fix( (lx - overlap)/nadvance ) ;
end
if (exist('k1') ~= 1) k1 = 0; end
if (exist('k2') ~= 1) k2 = 0; end
% ------ scale factors for unbiased estimates --------------------
nlags = 2 * maxlag + 1;
zlag = 1 + maxlag;
tmp = zeros(nlags,1);
if (flag(1:1) == 'b' || flag(1:1) == 'B')
scale = ones(nlags,1) / nsamp;
sc1 = 1/nsamp;
sc2 = sc1; sc12 = sc1;
else
ind = [-maxlag:maxlag]';
kmin = min(0,min(k1,k2));
kmax = max(0,max(k1,k2));
scale = nsamp - max(ind,kmax) + min(ind,kmin);
scale = ones(nlags,1) ./ scale;
sc1 = 1/(nsamp-abs(k1));
sc2 = 1/(nsamp-abs(k2));
sc12 = 1/(nsamp-abs(k1-k2));
disp(scale)
end
% ----------- estimate second- and fourth-order moments; combine ------
y_cum = zeros(2*maxlag+1,1) ;
rind = -maxlag:maxlag;
ind = 1:nsamp;
for i=1:nrecs
tmp = y_cum * 0;
R_zy = 0; R_wy = 0; M_wz = 0;
ws = w(ind); ws = ws(:) - mean(ws);
xs = x(ind); xs = xs(:) - mean(xs);
ys = y(ind); ys = ys(:) - mean(ys); cys = conj(ys);
zs = z(ind); zs = zs(:) - mean(zs);
ziv = xs * 0;
% create the "IV" matrix: offset for second lag
if (k1 >= 0)
ziv(1:nsamp-k1) = ws(1:nsamp-k1) .* cys(k1+1: nsamp);
R_wy = R_wy + ws(1:nsamp-k1)' * ys(k1+1:nsamp);
else
ziv(-k1+1:nsamp) = ws(-k1+1:nsamp) .* cys(1:nsamp+k1);
R_wy = R_wy + ws(-k1+1:nsamp)' * ys(1:nsamp+k1);
end
% create the "IV" matrix: offset for third lag
if (k2 >= 0)
ziv(1:nsamp-k2) = ziv(1:nsamp-k2) .* zs(k2+1: nsamp);
ziv(nsamp-k2+1:nsamp) = zeros(k2,1);
M_wz = M_wz + ws(1:nsamp-k2).' * zs(k2+1:nsamp);
else
ziv(-k2+1:nsamp) = ziv(-k2+1:nsamp) .* zs(1:nsamp+k2);
ziv(1:-k2) = zeros(-k2,1);
M_wz = M_wz + ws(-k2+1:nsamp).' * zs(1:nsamp-k2);
end
if (k1-k2 >= 0)
R_zy = R_zy + zs(1:nsamp-k1+k2)' * ys(k1-k2+1:nsamp);
else
R_zy = R_zy + zs(-k1+k2+1:nsamp)' * ys(1:nsamp-k2+k1);
end
tmp(zlag) = tmp(zlag) + ziv' * xs ;
for k = 1:maxlag
tmp(zlag-k) = tmp(zlag-k) + ziv([k+1:nsamp])' * xs([1:nsamp-k]);
tmp(zlag+k) = tmp(zlag+k) + ziv([1:nsamp-k])' * xs([k+1:nsamp]);
end
y_cum = y_cum + tmp .* scale ; % fourth-order moment estimates done
R_wx = cum2x(ws, xs, maxlag, nsamp, overlap0, flag);
R_zx = cum2x(zs, xs, maxlag+abs(k2), nsamp, overlap0, flag);
M_yx = cum2x(cys, xs, maxlag+abs(k1), nsamp, overlap0, flag);
y_cum = y_cum - R_zy * R_wx * sc12 ...
- R_wy * R_zx(rind - k2 + maxlag+abs(k2)+ 1) * sc1 ...
- M_wz'* M_yx(rind - k1 + maxlag+abs(k1)+ 1) * sc2 ;
ind = ind + nadvance ;
end
y_cum = y_cum / nrecs;
return
% mods Sept 7, 1998
% fixed conjugation errors in computing R_wy, R_zy, M_yx