-
Notifications
You must be signed in to change notification settings - Fork 30
/
Copy pathiact.m
40 lines (34 loc) · 791 Bytes
/
iact.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
function [tau,m] = iact(dati)
%IACT estimates the integrated autocorrelation time
% using Sokal's adaptive truncated periodogram estimator.
% Originally contributed by Antonietta Mira by name sokal.m
% Marko Laine <[email protected]>
% $Revision: 1.2 $ $Date: 2012/09/27 11:47:37 $
if length(dati) == prod(size(dati))
dati = dati(:);
end
[mx,nx] = size(dati);
tau = zeros(1,nx);
m = zeros(1,nx);
x = fft(dati);
xr = real(x);
xi = imag(x);
xr = xr.^2+xi.^2; %%%controllare questo
xr(1,:)=0;
xr=real(fft(xr));
var=xr(1,:)./length(dati)/(length(dati)-1);
for j = 1:nx
if var(j) == 0
continue
end
xr(:,j)=xr(:,j)./xr(1,j);
sum=-1/3;
for i=1:length(dati)
sum=sum+xr(i,j)-1/6;
if sum<0
tau(j)=2*(sum+(i-1)/6);
m(j)=i;
break
end
end
end