Signal- und Messdatenverarbeitung

Codes für stochastisch abgetastete, komplexe, periodische Einzelsignale, ohne Gewichtung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren von Wertepaaren (ti,xi) mit den Messzeitpunkten ti und den komplexen Messwerten xi (harmonisches Signal als Beispiel, mit Erwartungswert verschieden von null)
T=1000.0
dr=5.0
mxs=3.0+1.0j
axs=1.5
t=[]
x=[]
te=0.0
while te<T:
  tp=exponential(1.0/dr)
  te+=tp
  if te<T:
    t.append(te)
    x.append(axs*exp(0.1j*pi*te)+mxs)
  

t=array(t)
x=array(x)
plot(t,real(x),'o',t,imag(x),'o')
show()
T=1000;
dr=5;
mxs=3+1j;
axs=1.5;
t=[];
x=[];
te=0;
while te<T
tp=-log(1-rand())/dr;
te=te+tp;
if te<T
t(end+1)=te;
x(end+1)=axs*exp(0.1j*pi*te)+mxs;
end
end
plot(t,real(x),'o',t,imag(x),'o')
Mittelwert x‾=1N⋅∑i=0N−1xi
mean(x)
mean(x)
Varianz (ohne Bessel-Korrektur, asymptotisch erwartungstreu) s2=1N⋅∑i=0N−1(xi−x‾)*(xi−x‾)=1N⋅∑i=0N−1|xi−x‾|2
=(1N⋅∑i=0N−1xi*xi)−x‾*x‾=(1N⋅∑i=0N−1|xi|2)−|x‾|2
var(x)
var(x,1)
Autokorrelationsfunktion und Leistungsspektrum über Slotkorrelation (ohne Selbstprodukte durch bk(0)=0, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Rk=R(τk)=∑i=0N−1∑j=0N−1bk(tj−ti)xi*xj∑i=0N−1∑j=0N−1bk(tj−ti)
bk(t)={1falls 0<|t+cT|<Δt/2 für irgendein ganzes c0sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=[T/Δt]
Pj=P(fj)=1K⋅FFT{Rk}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
(imaginäre Einheit 𝐢)
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
N=len(t)
dt=1.0
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R1=zeros(K)+0j
R0=zeros(K)
for i in range(1,N):
  for j in range(0,i):
    k=int(round((t[j]-t[i]-T*round((t[j]-t[i])/float(T)))/float(dt)))
    R1[k]+=conj(x[i])*x[j]
    R0[k]+=1.0
  

for k in range(-(K//2),1):
  R1[k]+=conj(R1[-k])
  R0[k]+=R0[-k]

for k in range(1,(K+1)//2):
  R1[k]=conj(R1[-k])
  R0[k]=R0[-k]

R=R1/R0
plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R1=zeros(1,K)+0j;
R0=zeros(1,K);
for i=2:N
for j=1:i-1
k=round((t(j)-t(i)-T*round((t(j)-t(i))/T))/dt);
R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+conj(x(i))*x(j);
R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+1;
end
end
R1(1)=R1(1)+conj(R1(1));
R0(1)=2*R0(1);
for k=fix((K+1)/2)+1:K
R1(k)=R1(k)+conj(R1(K-k+2));
R0(k)=R0(k)+R0(K-k+2);
end
for k=2:fix((K+1)/2)
R1(k)=conj(R1(K-k+2));
R0(k)=R0(K-k+2);
end
R=R1./R0;
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Slotkorrelation (ohne Selbstprodukte durch bk(0)=0, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Rk=R(τk)=s2[∑i=0N−1∑j=0N−1bk(tj−ti)(xi−x‾)*(xj−x‾)][∑i=0N−1∑j=0N−1bk(tj−ti)|xi−x‾|2][∑i=0N−1∑j=0N−1bk(tj−ti)|xj−x‾|2]+|x‾|2
bk(t)={1falls 0<|t+cT|<Δt/2 für irgendein ganzes c0sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=[T/Δt]
Pj=P(fj)=1K⋅FFT{Rk}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
(imaginäre Einheit 𝐢)
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
N=len(t)
dt=1.0
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
mxe=mean(x)
vxe=var(x)
R1=zeros(K)+0j
R2=zeros(K)
R3=zeros(K)
for i in range(1,N):
  for j in range(0,i):
    k=int(round((t[j]-t[i]-T*round((t[j]-t[i])/float(T)))/float(dt)))
    R1[k]+=conj(x[i]-mxe)*(x[j]-mxe)
    R2[k]+=abs(x[i]-mxe)**2
    R3[k]+=abs(x[j]-mxe)**2
  

for k in range(-(K//2),1):
  R1[k]+=conj(R1[-k])
  R2[k]+=R3[-k]
  if k%K!=(-k)%K:
    R3[k]+=R2[-k]
  else:
    R3[k]=R2[-k]

for k in range(1,(K+1)//2):
  R1[k]=conj(R1[-k])
  R2[k]=R3[-k]
  R3[k]=R2[-k]

R=vxe*R1/sqrt(R2*R3)+abs(mxe)**2
plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
mxe=mean(x);
vxe=var(x,1);
R1=zeros(1,K)+0j;
R2=zeros(1,K);
R3=zeros(1,K);
for i=2:N
for j=1:i-1
k=round((t(j)-t(i)-T*round((t(j)-t(i))/T))/dt);
R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+conj(x(i)-mxe)*(x(j)-mxe);
R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+abs(x(i)-mxe)^2;
R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+abs(x(j)-mxe)^2;
end
end
R1(1)=R1(1)+conj(R1(1));
R2(1)=R2(1)+R3(1);
R3(1)=R2(1);
for k=fix((K+1)/2)+1:K
R1(k)=R1(k)+conj(R1(K-k+2));
R2(k)=R2(k)+R3(K-k+2);
if mod(k,K)~=mod(K-k+2,K)
R3(k)=R3(k)+R2(K-k+2);
else
R3(k)=R2(K-k+2);
end
end
for k=2:fix((K+1)/2)
R1(k)=conj(R1(K-k+2));
R2(k)=R3(K-k+2);
R3(k)=R2(K-k+2);
end
R=vxe*R1./sqrt(R2.*R3)+abs(mxe)^2;
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, ohne Normierung) Ej=E(fj)=T2N(N−1){|Xj|2−∑i=0N−1|xi|2}
Xj=∑i=0N−1xi𝐞−2π𝐢fjti
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J=[T/Δt]
RE,k=RE(τk)=1Δt⋅IFFT{Ej}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej𝐞2π𝐢fjτk/J
Rk=R(τk)=1TRE,k
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=J=[T/Δt]
Pj=P(fj)=1K⋅FFT{Rk}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
N=len(t)
dt=1.0
J=int(round(T/float(dt)))
fp=roll(arange(-(J//2),(J+1)//2)/float(J*dt),(J+1)//2)
X=zeros(size(fp))+0j
for i in range(0,N):
  X+=x[i]*exp(-2j*pi*fp*t[i])

E=T**2*(abs(X)**2-sum(abs(x)**2))/float(N**2-N)
RE=ifft(E)/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)+0j
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(T)

plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
J=round(T/dt);
fp=circshift((-fix(J/2):fix((J-1)/2))/(J*dt),[0;fix((J+1)/2)]);
X=zeros(size(fp))+0j;
for i=1:N
X=X+x(i)*exp(-2j*pi*fp*t(i));
end
E=T^2*(abs(X).^2-sum(abs(x).^2))/(N^2-N);
RE=ifft(E)/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros([1,K])+0j;
for k=-fix(K/2):fix((K-1)/2)
R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/T;
end
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Ej=E(fj)=T2N(N−1){|Xj|2−∑i=0N−1|xi|2}
Ej′=E′(fj)=T2N(N−1){|Xj′|2−N}
Xj=∑i=0N−1xi𝐞−2π𝐢fjti
Xj′=∑i=0N−1𝐞−2π𝐢fjti
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J=[T/Δt]
RE,k=RE(τk)=1Δt⋅IFFT{Ej}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej𝐞2π𝐢fjτk/J
RE,k′=RE′(τk)=1Δt⋅IFFT{Ej′}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej′𝐞2π𝐢fjτk/J
Rk=R(τk)=RE,kRE,k′
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=J=[T/Δt]
Pj=P(fj)=1K⋅FFT{Rk}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
N=len(t)
dt=1.0
J=int(round(T/float(dt)))
fp=roll(arange(-(J//2),(J+1)//2)/float(J*dt),(J+1)//2)
X=zeros(size(fp))+0j
Xp=zeros(size(fp))+0j
for i in range(0,N):
  X+=x[i]*exp(-2j*pi*fp*t[i])
  Xp+=exp(-2j*pi*fp*t[i])

E=T**2*(abs(X)**2-sum(abs(x)**2))/float(N**2-N)
Ep=T**2*(abs(Xp)**2-N)/float(N**2-N)
RE=ifft(E)/float(dt)
REp=real(ifft(Ep))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)+0j
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(REp[k])

plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
J=round(T/dt);
fp=circshift((-fix(J/2):fix((J-1)/2))/(J*dt),[0;fix((J+1)/2)]);
X=zeros(size(fp))+0j;
Xp=zeros(size(fp))+0j;
for i=1:N
X=X+x(i)*exp(-2j*pi*fp*t(i));
Xp=Xp+exp(-2j*pi*fp*t(i));
end
E=T^2*(abs(X).^2-sum(abs(x).^2))/(N^2-N);
Ep=T^2*(abs(Xp).^2-N)/(N^2-N);
RE=ifft(E)/dt;
REp=real(ifft(Ep))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros([1,K])+0j;
for k=-fix(K/2):fix((K-1)/2)
R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/REp(mod(k+length(REp),length(REp))+1);
end
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Ej=E(fj)=T2N(N−1){|Xj|2−∑i=0N−1|xi−x‾|2}
Ej′=E′(fj)=T2N(N−1){Xj′′*Xj′−∑i=0N−1|xi−x‾|2}
Ej′′=E′′(fj)=T2N(N−1){Xj′*Xj′′−∑i=0N−1|xi−x‾|2}
Xj=∑i=0N−1(xi−x‾)𝐞−2π𝐢fjti
Xj′=∑i=0N−1𝐞−2π𝐢fjti
Xj′′=∑i=0N−1|xi−x‾|2𝐞−2π𝐢fjti
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J=[T/Δt]
RE,k=RE(τk)=1Δt⋅IFFT{Ej}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej𝐞2π𝐢fjτk/J
RE,k′=RE′(τk)=1Δt⋅IFFT{Ej′}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej′𝐞2π𝐢fjτk/J
RE,k′′=RE′′(τk)=1Δt⋅IFFT{Ej′′}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej′′𝐞2π𝐢fjτk/J
Rk=R(τk)=s2RE,kRE,k′RE,k′′+|x‾|2
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=J=[T/Δt]
Pj=P(fj)=1K⋅FFT{Rk}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
N=len(t)
dt=1.0
mxe=mean(x)
vxe=var(x)
J=int(round(T/float(dt)))
fp=roll(arange(-(J//2),(J+1)//2)/float(J*dt),(J+1)//2)
X=zeros(size(fp))+0j
Xp=zeros(size(fp))+0j
Xpp=zeros(size(fp))+0j
for i in range(0,N):
  X+=(x[i]-mxe)*exp(-2j*pi*fp*t[i])
  Xp+=exp(-2j*pi*fp*t[i])
  Xpp+=abs(x[i]-mxe)**2*exp(-2j*pi*fp*t[i])

E=T**2*(abs(X)**2-sum(abs(x-mxe)**2))/float(N**2-N)
Ep=T**2*(conj(Xpp)*Xp-sum(abs(x-mxe)**2))/float(N**2-N)
Epp=T**2*(conj(Xp)*Xpp-sum(abs(x-mxe)**2))/float(N**2-N)
RE=ifft(E)/float(dt)
REp=real(ifft(Ep))/float(dt)
REpp=real(ifft(Epp))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)+0j
for k in range(-(K//2),(K+1)//2):
  R[k]=vxe*RE[k]/sqrt(REp[k]*REpp[k])+abs(mxe)**2

plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
mxe=mean(x);
vxe=var(x,1);
J=round(T/dt);
fp=circshift((-fix(J/2):fix((J-1)/2))/(J*dt),[0;fix((J+1)/2)]);
X=zeros(size(fp))+0j;
Xp=zeros(size(fp))+0j;
Xpp=zeros(size(fp))+0j;
for i=1:N
X=X+(x(i)-mxe)*exp(-2j*pi*fp*t(i));
Xp=Xp+exp(-2j*pi*fp*t(i));
Xpp=Xpp+abs(x(i)-mxe)^2*exp(-2j*pi*fp*t(i));
end
E=T^2*(abs(X).^2-sum(abs(x-mxe).^2))/(N^2-N);
Ep=T^2*(conj(Xpp).*Xp-sum(abs(x-mxe).^2))/(N^2-N);
Epp=T^2*(conj(Xp).*Xpp-sum(abs(x-mxe).^2))/(N^2-N);
RE=ifft(E)/dt;
REp=real(ifft(Ep))/dt;
REpp=real(ifft(Epp))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros([1,K])+0j;
for k=-fix(K/2):fix((K-1)/2)
R(mod(k+K,K)+1)=vxe*RE(mod(k+length(RE),length(RE))+1)/sqrt(REp(mod(k+length(REp),length(REp))+1)*REpp(mod(k+length(REpp),length(REpp))+1))+abs(mxe)^2;
end
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
(keine Lomb-Scargle-Methode für komplexe Werte in Matlab und Python)
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, ohne Normierung)
from numpy.fft import *
N=len(t)
dt=1.0
J=int(round(T/float(dt)))
x1=zeros(J)+0j
for i in range(0,N):
  j=int(floor((t[i]-T*floor(t[i]/float(T)))/float(dt)))
  x1[j]+=x[i];

X=fft(x1)
E=T**2*(abs(X)**2-sum(abs(x)**2))/float(N**2-N)
RE=ifft(E)/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)+0j
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(T)

plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
J=round(T/dt);
x1=zeros([1,J])+0j;
for i=1:N
j=fix((t(i)-T*fix(t(i)/T))/dt)+1;
x1(j)=x1(j)+x(i);
end
X=fft(x1);
E=T^2*(abs(X).^2-sum(abs(x).^2))/(N^2-N);
RE=ifft(E)/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros([1,K])+0j;
for k=-fix(K/2):fix((K-1)/2)
R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/T;
end
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)
from numpy.fft import *
N=len(t)
dt=1.0
J=int(round(T/float(dt)))
x0=zeros(J)
x1=zeros(J)+0j
for i in range(0,N):
  j=int(floor((t[i]-T*floor(t[i]/float(T)))/float(dt)))
  x0[j]+=1.0;
  x1[j]+=x[i];

X=fft(x1)
Xp=fft(x0)
E=T**2*(abs(X)**2-sum(abs(x)**2))/float(N**2-N)
Ep=T**2*(abs(Xp)**2-N)/float(N**2-N)
RE=ifft(E)/float(dt)
REp=real(ifft(Ep))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)+0j
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(REp[k])

plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
J=round(T/dt);
x0=zeros([1,J]);
x1=zeros([1,J])+0j;
for i=1:N
j=fix((t(i)-T*fix(t(i)/T))/dt)+1;
x0(j)=x0(j)+1;
x1(j)=x1(j)+x(i);
end
X=fft(x1);
Xp=fft(x0);
E=T^2*(abs(X).^2-sum(abs(x).^2))/(N^2-N);
Ep=T^2*(abs(Xp).^2-N)/(N^2-N);
RE=ifft(E)/dt;
REp=real(ifft(Ep))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros([1,K])+0j;
for k=-fix(K/2):fix((K-1)/2)
R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/REp(mod(k+length(REp),length(REp))+1);
end
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)
from numpy.fft import *
N=len(t)
dt=1.0
mxe=mean(x)
vxe=var(x)
J=int(round(T/float(dt)))
x0=zeros(J)
x1=zeros(J)+0j
x2=zeros(J)
for i in range(0,N):
  j=int(floor((t[i]-T*floor(t[i]/float(T)))/float(dt)))
  x0[j]+=1.0;
  x1[j]+=x[i]-mxe;
  x2[j]+=abs(x[i]-mxe)**2;

X=fft(x1)
Xp=fft(x0)
Xpp=fft(x2)
E=T**2*(abs(X)**2-sum(abs(x-mxe)**2))/float(N**2-N)
Ep=T**2*(conj(Xpp)*Xp-sum(abs(x-mxe)**2))/float(N**2-N)
Epp=T**2*(conj(Xp)*Xpp-sum(abs(x-mxe)**2))/float(N**2-N)
RE=ifft(E)/float(dt)
REp=real(ifft(Ep))/float(dt)
REpp=real(ifft(Epp))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)+0j
for k in range(-(K//2),(K+1)//2):
  R[k]=vxe*RE[k]/sqrt(REp[k]*REpp[k])+abs(mxe)**2

plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
mxe=mean(x);
vxe=var(x,1);
J=round(T/dt);
x0=zeros([1,J]);
x1=zeros([1,J])+0j;
x2=zeros([1,J]);
for i=1:N
j=fix((t(i)-T*fix(t(i)/T))/dt)+1;
x0(j)=x0(j)+1;
x1(j)=x1(j)+x(i)-mxe;
x2(j)=x2(j)+abs(x(i)-mxe)^2;
end
X=fft(x1);
Xp=fft(x0);
Xpp=fft(x2);
E=T^2*(abs(X).^2-sum(abs(x-mxe).^2))/(N^2-N);
Ep=T^2*(conj(Xpp).*Xp-sum(abs(x-mxe).^2))/(N^2-N);
Epp=T^2*(conj(Xp).*Xpp-sum(abs(x-mxe).^2))/(N^2-N);
RE=ifft(E)/dt;
REp=real(ifft(Ep))/dt;
REpp=real(ifft(Epp))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros([1,K])+0j;
for k=-fix(K/2):fix((K-1)/2)
R(mod(k+K,K)+1)=vxe*RE(mod(k+length(RE),length(RE))+1)/sqrt(REp(mod(k+length(REp),length(REp))+1)*REpp(mod(k+length(REpp),length(REpp))+1))+abs(mxe)^2;
end
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)
from numpy.fft import *
N=len(t)
dt=1.0
J=int(round(T/float(dt)))
xp=x[N-1]*ones(J)
for i in range(0,N-1):
  for j in range(int(floor(t[i]/float(dt)))+1,int(floor(t[i+1]/float(dt)))+1):
    if j<J:
      xp[j]=x[i]
    
  

Xp=fft(xp)
Ep=dt**2*abs(Xp)**2
REp=ifft(Ep)/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)+0j
c=exp(-dd/float(dt))/(1-exp(-dd/float(dt)))**2
for k in range(-(K//2),(K+1)//2):
  if k==0:
    R[k]=REp[k]/float(T)
  else:
    R[k]=((2*c+1)*REp[k]-c*REp[k-1]-c*REp[k+1])/float(T)
  

plot(tau,real(R),'o',tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
P=real(fft(R))/float(K)
plot(f,P,'o')
show()
N=length(t);
dt=1.0;
J=round(T/dt);
xp=x(N)*ones([1,J]);
for i=1:N-1
for j=fix(t(i)/dt)+2:fix(t(i+1)/dt)+1
if j<=J
xp(j)=x(i);
end
end
end
Xp=fft(xp);
Ep=dt^2*abs(Xp).^2;
REp=ifft(Ep)/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros([1,K])+0j;
c=exp(-dd/dt)/(1-exp(-dd/dt))^2;
for k=-fix(K/2):fix((K-1)/2)
if k==0
R(mod(k+K,K)+1)=REp(mod(k+length(REp),length(REp))+1)/T;
else
R(mod(k+K,K)+1)=((2*c+1)*REp(mod(k+length(REp),length(REp))+1)-c*REp(mod(k-1+length(REp),length(REp))+1)-c*REp(mod(k+1+length(REp),length(REp))+1))/T;
end
end
plot(tau,real(R),'o',tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
P=real(fft(R))/K;
plot(f,P,'o')