% Autokorrelationsfunktion und Autoleistungsspektrum fuer ungleichmaessig abgetastete Signale mit wirklich streng periodischen Anteilen
% mit Zeitquentisierung und FFT-Implementierung, hier ohne Gewichtung, Gewichtung in sperfftw.m
% Slotkorrelation, direkte Spektralschaetzung und Zeitquantisierung ergeben leicht unterschiedliche Ergebnisse
function [tau,R,f,P]=sperfft(x,t,dt,NT)
% x - Datensatz
% t - Messzeitpunkte (gleiche Laenge wie x)
% dt - Zeitschritt (fuer Korrelationsfunktion und Spektrum)
% NT - angenommene Periodenlaenge (Vielfache des Zeitschritts)
N=length(x);
df=1/(NT*dt);
f=(mod((0:NT-1)+fix(NT/2),NT)-fix(NT/2))*df;
xp=zeros(1,NT);
np=zeros(1,NT);
for i=1:N
k=mod(floor(t(i)/dt),NT)+1;
xp(k)=xp(k)+x(i);
np(k)=np(k)+1;
X1=fft(xp);
X0=fft(np);
P1=abs(X1).^2/(NT^2);
P0=abs(X0).^2/(NT^2);
R1=NT*real(ifft(P1));
R0=NT*real(ifft(P0));
P=real(fft(R1./R0))/NT;
tau=(0:NT-1)*dt;
R=NT*real(ifft(P));
end