% Autokorrelationsfunktion und Autoleistungsspektrum fuer ungleichmaessig abgetastete Signale mit nicht ganz streng periodischen Anteilen
% mit direkter Spektralschaetzung (als Referenz - langsam), Mittelung ueber einen Bereich von maximal M Werten, mit Gewichtung
% Slotkorrelation, direkte Spektralschaetzung und Zeitquantisierung ergeben leicht unterschiedliche Ergebnisse
function [tau,R,f,P]=snperdirw(x,t,w,dt,NT,M)
% x - Datensatz
% t - Messzeitpunkte (gleiche Laenge wie x)
% w - Gewichte (gleiche Laenge wie x)
% dt - Zeitschritt (fuer Korrelationsfunktion und Spektrum)
% NT - angenommene Periodenlaenge (Vielfache des Zeitschritts)
% M - maximale Bereich der Mittelung (z.B. =NT fuer Mittelung ueber eine Periode)
N=length(x);
df=1/((N+M)*dt);
fi=(mod((0:N+M-1)+fix((N+M)/2),N+M)-fix((N+M)/2))*df;
X1=zeros(1,N+M)+0j;
X0=zeros(1,N+M)+0j;
for i=1:N
E=exp(-2j*pi*fi*t(i));
X1=X1+w(i)*x(i)*E;
X0=X0+w(i)*E;
P1=abs(X1).^2/(M^2);
P0=abs(X0).^2/(M^2);
R1=M*real(ifft(P1));
R0=M*real(ifft(P0));
R1x=zeros(1,NT);
R0x=zeros(1,NT);
for k=-M:M-1
R1x(mod(k,NT)+1)=R1x(mod(k,NT)+1)+R1(mod(k,N+M)+1);
R0x(mod(k,NT)+1)=R0x(mod(k,NT)+1)+R0(mod(k,N+M)+1);
end
df=1/(NT*dt);
f=(mod((0:NT-1)+fix(NT/2),NT)-fix(NT/2))*df;
P=real(fft(R1x./R0x))/NT;
tau=(0:NT-1)*dt;
R=NT*real(ifft(P));
end