forked from grinsted/gwmcmc
-
Notifications
You must be signed in to change notification settings - Fork 0
/
eacorr.m
69 lines (62 loc) · 1.57 KB
/
eacorr.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
function [C,lags]=eacorr(m)
%% EACORR, Ensemble auto correlation
%
% Ensemble average auto correlation.
%
% USAGE: [C,lags]=eacorr(m)
%
% eacorr is designed to be applied on the MxWxT matrices output from gwmcmc,
% but can also be applied to 2d and 1d matrices.
%
% When m is a 3d matrix of size MxWxT:
% * the auto-correlation will be calculated along the T-dimension
% * the ensemble averaging will be done along the W-dimension
% * the auto-correlation will be calculated for each of the M-dimensions
%
% When m is a 1d or 2d matrix, then it will be assumed that there is just a
% single ensemble member, and the autocorrelation will be measured along
% the longest dimension.
%
% (The centering of the input series is done using the ensemble mean rather
% than the mean of each individual series.)
%
% example:
% eacorr(smooth(randn(1000,)))
%
% Aslak Grinsted 2015
sM=size(m);
sM(end+1:3)=1;
if sM(3)==1
if sM(1)>sM(2)
m=permute(m,[2 3 1]);
else
m=permute(m,[1 3 2]);
end
end
M=size(m,1);W=size(m,2);T=size(m,3);
lags=(0:T-1)';
nfft = 2^nextpow2(2*T-1);
C=nan(T,M);
for mix=1:M
c=zeros(T,1);
mm=mean(m(mix,:)); %center data using ensemble average!
for wix=1:W
d=m(mix,wix,:);
d=d(:)-mm;
r=ifft( abs(fft(d,nfft)).^2);
c=c+r(1:T);
end
C(:,mix)=c./(T-lags); %biased/unbiased
C(:,mix)=C(:,mix)./C(1,mix);
end
if isreal(m)
C=real(C);
end
if nargout==0
figure
plot(lags,C,'.-',lags([1 end]),[0 0],'k');
grid on
xlabel('lags')
ylabel('autocorrelation');
clearvars lags C
end