Signal- und Messdatenverarbeitung

Codes für stochastisch abgetastete, komplexe Einzelleistungssignale, 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 (zwei AR1-Prozess als Beispiel, mit Erwartungswert verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)
T=10000.0
Ti=10.0
dr=1.0
mxs=3.0+1.0j
vxs=1.0
t=[]
x=[]
te=0.0
xe=normal(real(mxs),sqrt(vxs/2.0))+1j*normal(imag(mxs),sqrt(vxs/2.0))
while te<T:
  tp=exponential(1.0/(2*dr))
  te+=tp
  phi=exp(-tp/Ti)
  theta=sqrt((1-phi**2)*vxs/2.0)
  xe=(xe-mxs)*phi+normal(real(mxs),theta)+1j*normal(imag(mxs),theta)
  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=10000;
Ti=10;
dr=1;
mxs=3+1j;
vxs=1;
t=[];
x=[];
te=0;
xe=real(mxs)+sqrt(vxs/2)*randn()+1j*(imag(mxs)+sqrt(vxs/2)*randn());
while te<T
tp=-log(1-rand())/(2*dr);
te=te+tp;
phi=exp(-tp/Ti);
theta=sqrt((1-phi^2)*vxs/2);
xe=(xe-mxs)*phi+real(mxs)+theta*randn()+1j*(imag(mxs)+theta*randn());
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 Leistungsdichtespektrum über Slotkorrelation (ohne Selbstprodukte durch bk(0)=0) 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|<Δt/20sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R1=zeros(K)+0j
R0=zeros(K)
i=1
while i<N:
  j=i-1
  while (j>=0) and (t[j]>t[i]-(K//2+0.5)*dt):
    k=int(round((t[j]-t[i])/float(dt)))
    R1[k]+=g[i]*g[j]*conj(x[i])*x[j]
    R0[k]+=g[i]*g[j]
    j-=1
  
  i+=1

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',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R1=zeros(1,K)+0j;
R0=zeros(1,K);
i=2;
while i<=N
j=i-1;
while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)*dt)
k=round((t(j)-t(i))/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);
j=j-1;
end
i=i+1;
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',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Slotkorrelation (ohne Selbstprodukte durch bk(0)=0, mit lokaler Normierung) 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|<Δt/20sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
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)
i=1
while i<N:
  j=i-1
  while (j>=0) and (t[j]>t[i]-(K//2+0.5)*dt):
    k=int(round((t[j]-t[i])/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
    j-=1
  
  i+=1

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',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
K=200;
K=min(ceil(2*T/dt)-1,K);
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);
i=2;
while i<=N
j=i-1;
while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)*dt)
k=round((t(j)-t(i))/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;
j=j-1;
end
i=i+1;
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',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum ü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≥⌈2T/Δ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)=RE,kT−|τk|
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑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(ceil(2*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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
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-abs(tau[k]))

plot(tau,real(R),'o',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*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=200;
K=min(ceil(2*T/dt)-1,K);
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-abs(tau(mod(k+K,K)+1)));
end
plot(tau,real(R),'o',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit Normierung) 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≥⌈2T/Δ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<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑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(ceil(2*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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*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=200;
K=min(ceil(2*T/dt)-1,K);
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',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit lokaler Normierung) 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≥⌈2T/Δ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<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑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(ceil(2*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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
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=ceil(2*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=200;
K=min(ceil(2*T/dt)-1,K);
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',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
(keine passende Lomb-Scargle-Methode in Matlab und Python; für Gewichtung s. The generalised Lomb-Scargle periodogram.)
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, ohne Normierung)
from numpy.fft import *
N=len(t)
dt=1.0
J=int(ceil(2*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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
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-abs(tau[k]))

plot(tau,real(R),'o',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*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=200;
K=min(ceil(2*T/dt)-1,K);
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-abs(tau(mod(k+K,K)+1)));
end
plot(tau,real(R),'o',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit Normierung)
from numpy.fft import *
N=len(t)
dt=1.0
J=int(ceil(2*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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*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=200;
K=min(ceil(2*T/dt)-1,K);
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',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit lokaler Normierung)
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(ceil(2*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=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+abs(mxs)**2,tau,imag(R),'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'o',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
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=ceil(2*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=200;
K=min(ceil(2*T/dt)-1,K);
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',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'o',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
(keine Gewichtung für Interpolation)