forked from ossadtchi/PSIICOS
-
Notifications
You must be signed in to change notification settings - Fork 0
/
CrossSpectralTimeseries.m
64 lines (50 loc) · 1.29 KB
/
CrossSpectralTimeseries.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
function [ A, key,Coh ] = CrossSpectralTimeseries(X,bInducedOnly)
if(nargin==1)
bInducedOnly = false;
end;
[Nch, Ns Ntr] = size(X);
if(bInducedOnly)
ERP = mean(X,3);
X = X-bsxfun(@minus, X,ERP);
end;
Xfft = fft(X,[],2);
h = zeros(1,Ns); % nx1 for nonempty. 0x0 for empty.
if Ns > 0 && 2*fix(Ns/2) == Ns
% even and nonempty
h([1 Ns/2+1]) = 1;
h(2:Ns/2) = 2;
elseif Ns>0
% odd and nonempty
h(1) = 1;
h(2:(Ns+1)/2) = 2;
end
HF = repmat(h,[Nch,1,Ntr]);
XH = ifft(Xfft.*HF,[],2);
Xph = XH; %./(abs(XH)+0.0001*mean(abs(XH(:))));
clear XH;
A = zeros(Nch*(Nch-1)/2,Ns);
Nch_2 = Nch/2;
XphConj = conj(Xph);
Ntr = size(Xph,3);
range = 1:Nch-1;
trs = 1:Ntr;
k = 1;
KEY = reshape(1:Nch*Nch,Nch,Nch);
% we will take the diagonal as well
A = zeros(Nch*(Nch+1)/2,Ns);
fprintf('Calculating vectorised form of the cross spectral matrix upper triangle ... \n');
fprintf('Reference sensor (max %d): ',Nch);
Coh = 1;
for i=1:Nch
mn = (mean( bsxfun(@times,(Xph(1:Nch,:,:)),XphConj(i,:,:)),3));
A(k:k+Nch-1,:) = mn;
key(k:k+Nch-1) = KEY(1:Nch,i);
k = k+Nch;
if i>1
for j=0:log10(i-1)
fprintf('\b'); % delete previous counter display
end
end
fprintf('%d', i);
end;
fprintf('\n');