% Autokorrelationsfunktion und Autoleistungsspektrum fuer ungleichmaessig abgetastete Signale mit nicht ganz streng periodischen Anteilen
% mit Zeitquentisierung und FFT-Implementierung, Mittelung ueber einen Bereich von maximal M Werten, hier ohne Gewichtung, Gewichtung in snperfftw.m
% Slotkorrelation, direkte Spektralschaetzung und Zeitquantisierung ergeben leicht unterschiedliche Ergebnisse
function [tau,R,f,P]=snperfft(x,t,dt,NT,M)
  % x - Datensatz
  % t - Messzeitpunkte (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/(NT*dt);
  f=(mod((0:NT-1)+fix(NT/2),NT)-fix(NT/2))*df;
  xp=zeros(1,2^nextpow2(ceil(t(end)/dt)+M));
  np=zeros(1,2^nextpow2(ceil(t(end)/dt)+M));
  for i=1:N
    xp(floor(t(i)/dt)+1)=xp(floor(t(i)/dt)+1)+x(i);
    np(floor(t(i)/dt)+1)=np(floor(t(i)/dt)+1)+1;
  end
  X1=fft(xp);
  X0=fft(np);
  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,2^nextpow2(N+M))+1);
    R0x(mod(k,NT)+1)=R0x(mod(k,NT)+1)+R0(mod(k,2^nextpow2(N+M))+1);
  end
  P=real(fft(R1x./R0x))/NT;
  tau=(0:NT-1)*dt;
  R=NT*real(ifft(P));
end

