Signal- und Messdatenverarbeitung

Codes für stochastisch abgetastete, komplexe, periodische Einzelsignale, mit 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) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)
T=1000.0
dr=5.0
mxs=3.0+1.0j
axs=1.0
t=[]
x=[]
te=0.0
while te<T:
  tp=exponential(1.0/(2*dr))
  te+=tp
  xe=axs*exp(0.1j*pi*te)+mxs
  if (te<T) and (random()<1/(1+exp(-real(xe-mxs)))):
    t.append(te)
    x.append(xe)
  

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;
t=[];
x=[];
te=0;
while te<T
tp=-log(1-rand())/(2*dr);
te=te+tp;
xe=axs*exp(0.1j*pi*te)+mxs;
if (te<T) && (rand<1/(1+exp(-real(xe-mxs))))
t(end+1)=te;
x(end+1)=xe;
end
end
plot(t,real(x),'o',t,imag(x),'o')
Generieren von reellen Gewichten gi passend zu den (ti,xi) (Inverse der verschobenen Sigmoidfunktion)
g=1+exp(-real(x-mxs))
plot(t,g,'o')
show()
g=1+exp(-real(x-mxs));
plot(t,g,'o')
Mittelwert x‾=∑i=0N−1gixi∑i=0N−1gi
sum(g*x)/float(sum(g))
sum(g.*x)/sum(g)
Varianz (ohne Bessel-Korrektur, asymptotisch erwartungstreu) s2=∑i=0N−1gi(xi−x‾)*(xi−x‾)∑i=0N−1gi=∑i=0N−1gi|xi−x‾|2∑i=0N−1gi
=∑i=0N−1gixi*xi∑i=0N−1gi−x‾*x‾=∑i=0N−1gi|xi|2∑i=0N−1gi−|x‾|2
sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2
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)gigjxi*xj∑i=0N−1∑j=0N−1bk(tj−ti)gigj
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]+=g[i]*g[j]*conj(x[i])*x[j]
    R0[k]+=g[i]*g[j]
  

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)+g(i)*g(j)*conj(x(i))*x(j);
R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+g(i)*g(j);
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)gigj(xi−x‾)*(xj−x‾)][∑i=0N−1∑j=0N−1bk(tj−ti)gigj|xi−x‾|2][∑i=0N−1∑j=0N−1bk(tj−ti)gigj|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=sum(g*x)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
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]+=g[i]*g[j]*conj(x[i]-mxe)*(x[j]-mxe)
    R2[k]+=g[i]*g[j]*abs(x[i]-mxe)**2
    R3[k]+=g[i]*g[j]*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=sum(g.*x)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
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)+g(i)*g(j)*conj(x(i)-mxe)*(x(j)-mxe);
R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+g(i)*g(j)*abs(x(i)-mxe)^2;
R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+g(i)*g(j)*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)=T2|Xj|2−∑i=0N−1gi2|xi|2(∑i=0N−1gi)2−∑i=0N−1gi2
Xj=∑i=0N−1gixi𝐞−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+=g[i]*x[i]*exp(-2j*pi*fp*t[i])

E=T**2*(abs(X)**2-sum((g*abs(x))**2))/float(sum(g)**2-sum(g**2))
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+g(i)*x(i)*exp(-2j*pi*fp*t(i));
end
E=T^2*(abs(X).^2-sum((g.*abs(x)).^2))/(sum(g)^2-sum(g.^2));
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)=T2|Xj|2−∑i=0N−1gi2|xi|2X0′2−∑i=0N−1gi2
Ej′=E′(fj)=T2|Xj′|2−∑i=0N−1gi2X0′2−∑i=0N−1gi2
Xj=∑i=0N−1gixi𝐞−2π𝐢fjti
Xj′=∑i=0N−1gi𝐞−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+=g[i]*x[i]*exp(-2j*pi*fp*t[i])
  Xp+=g[i]*exp(-2j*pi*fp*t[i])

E=T**2*(abs(X)**2-sum((g*abs(x))**2))/(abs(Xp[0])**2-sum(g**2))
Ep=T**2*(abs(Xp)**2-sum(g**2))/(abs(Xp[0])**2-sum(g**2))
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+g(i)*x(i)*exp(-2j*pi*fp*t(i));
Xp=Xp+g(i)*exp(-2j*pi*fp*t(i));
end
E=T^2*(abs(X).^2-sum((g.*abs(x)).^2))/(abs(Xp(1))^2-sum(g.^2));
Ep=T^2*(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2));
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)=T2|Xj|2−∑i=0N−1gi2|xi−x‾|2X0′2−∑i=0N−1gi2
Ej′=E′(fj)=T2Xj′′*Xj′−∑i=0N−1gi2|xi−x‾|2X0′2−∑i=0N−1gi2
Ej′′=E′′(fj)=T2Xj′*Xj′′−∑i=0N−1gi2|xi−x‾|2X0′2−∑i=0N−1gi2
Xj=∑i=0N−1gi(xi−x‾)𝐞−2π𝐢fjti
Xj′=∑i=0N−1gi𝐞−2π𝐢fjti
Xj′′=∑i=0N−1gi|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=sum(g*x)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
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+=g[i]*(x[i]-mxe)*exp(-2j*pi*fp*t[i])
  Xp+=g[i]*exp(-2j*pi*fp*t[i])
  Xpp+=g[i]*abs(x[i]-mxe)**2*exp(-2j*pi*fp*t[i])

E=T**2*(abs(X)**2-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
Ep=T**2*(conj(Xpp)*Xp-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
Epp=T**2*(conj(Xp)*Xpp-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
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=sum(g.*x)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
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+g(i)*(x(i)-mxe)*exp(-2j*pi*fp*t(i));
Xp=Xp+g(i)*exp(-2j*pi*fp*t(i));
Xpp=Xpp+g(i)*abs(x(i)-mxe)^2*exp(-2j*pi*fp*t(i));
end
E=T^2*(abs(X).^2-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
Ep=T^2*(conj(Xpp).*Xp-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
Epp=T^2*(conj(Xp).*Xpp-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
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 passende Lomb-Scargle-Methode in Matlab und Python; für Gewichtung s. The generalised Lomb-Scargle periodogram.)
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]+=g[i]*x[i];

X=fft(x1)
E=T**2*(abs(X)**2-sum((g*abs(x))**2))/float(sum(g)**2-sum(g**2))
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)+g(i)*x(i);
end
X=fft(x1);
E=T^2*(abs(X).^2-sum((g.*abs(x)).^2))/(sum(g)^2-sum(g.^2));
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]+=g[i];
  x1[j]+=g[i]*x[i];

X=fft(x1)
Xp=fft(x0)
E=T**2*(abs(X)**2-sum((g*abs(x))**2))/(abs(Xp[0])**2-sum(g**2))
Ep=T**2*(abs(Xp)**2-sum(g**2))/(abs(Xp[0])**2-sum(g**2))
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)+g(i);
x1(j)=x1(j)+g(i)*x(i);
end
X=fft(x1);
Xp=fft(x0);
E=T^2*(abs(X).^2-sum((g.*abs(x)).^2))/(abs(Xp(1))^2-sum(g.^2));
Ep=T^2*(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2));
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=sum(g*x)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
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]+=g[i];
  x1[j]+=g[i]*(x[i]-mxe);
  x2[j]+=g[i]*abs(x[i]-mxe)**2;

X=fft(x1)
Xp=fft(x0)
Xpp=fft(x2)
E=T**2*(abs(X)**2-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
Ep=T**2*(conj(Xpp)*Xp-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
Epp=T**2*(conj(Xp)*Xpp-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
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=sum(g.*x)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
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)+g(i);
x1(j)=x1(j)+g(i)*(x(i)-mxe);
x2(j)=x2(j)+g(i)*abs(x(i)-mxe)^2;
end
X=fft(x1);
Xp=fft(x0);
Xpp=fft(x2);
E=T^2*(abs(X).^2-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
Ep=T^2*(conj(Xpp).*Xp-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
Epp=T^2*(conj(Xp).*Xpp-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
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 Gewichtung für Interpolation)