Signal- und Messdatenverarbeitung

Codes für komplexe, stochastische Einzelleistungssignale, mit Gewichtung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren von N komplexen Werten xi mit i=0…N−1 (zwei AR1-Prozesse als Beispiel, mit einem Erwartungswert verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)
N=1000
dt=1.0
t=arange(0,N)*dt
phi=0.95
mxs=3+1j
vxs=1
theta=sqrt((1-phi**2)*vxs/2.0)
x=zeros(N)+0j
x[0]=normal(real(mxs),sqrt(vxs/2.0))+1j*normal(imag(mxs),sqrt(vxs/2.0))
while random()>1/(1+exp(-abs(x[0]-mxs))):
  x[0]=normal(real(mxs),sqrt(vxs/2.0))+1j*normal(imag(mxs),sqrt(vxs/2.0))

for i in range(1,N):
  x[i]=phi*(x[i-1]-mxs)+normal(real(mxs),theta)+1j*normal(imag(mxs),theta)
  while random()>1/(1+exp(-abs(x[i]-mxs))):
    x[i]=phi*(x[i-1]-mxs)+normal(real(mxs),theta)+1j*normal(imag(mxs),theta)
  

plot(t,real(x),'o',t,imag(x),'o')
show()
N=1000;
dt=1.0;
t=(0:N-1)*dt;
phi=0.95;
mxs=3+1j;
vxs=1;
theta=sqrt((1-phi^2)*vxs/2);
x=zeros(1,N)+0j;
x(1)=real(mxs)+sqrt(vxs/2)*randn()+1j*(imag(mxs)+sqrt(vxs/2)*randn());
while rand>1/(1+exp(-abs(x(1)-mxs)))
x(1)=real(mxs)+sqrt(vxs/2)*randn()+1j*(imag(mxs)+sqrt(vxs/2)*randn());
end
for i=2:N
x(i)=phi*(x(i-1)-mxs)+real(mxs)+theta*randn()+1j*(imag(mxs)+theta*randn());
while rand>1/(1+exp(-abs(x(i)-mxs)))
x(i)=phi*(x(i-1)-mxs)+real(mxs)+theta*randn()+1j*(imag(mxs)+theta*randn());
end
end
plot(t,real(x),'o',t,imag(x),'o')
Generieren von N reellen Gewichten gi mit i=0…N−1, passend zu den xi mit i=0…N−1 (Inverse der verschobenen Sigmoidfunktion)
g=zeros(N)
for i in range(0,N):
  g[i]=(1+exp(-abs(x[i]-mxs)))**(1+phi)

plot(t,g,'o')
show()
g=zeros(1,N);
for i=1:N
g(i)=(1+exp(-abs(x(i)-mxs)))**(1+phi);
end
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
Varianz (mit Bessel-Korrektur für korrelierte Daten, erwartungstreu) s2=∑i=0N−1gi(xi−x‾)2∑i=0N−1gi+sx‾2=∑i=0N−1gixi2∑i=0N−1gi−x‾2+sx‾2
mit der erwartungstreuen Schätzung sx‾2 der Schätzvarianz des Mittelwertschätzers für korrelierte Daten aus der erwartungstreuen Schätzung der (Auto-)Kovarianzfunktion mit Bessel-Korrektur für korrelierte Daten (s. unten)
from numpy.fft import *
from numpy.linalg import *
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N
mxe=sum(g*x)/float(sum(g))
xe=append(g*(x-mxe),zeros(N))
ge=append(g,zeros(N))
X=ifft(abs(fft(xe))**2)
Y=real(ifft(abs(fft(ge))**2))
Cp=zeros(2*Nc-1)+0j
for k in range(-(Nc-1),Nc):
  Cp[k]=X[k]/Y[k]

D=sum(g)
A=zeros([2*Nc-1,2*Nc-1])
for k in range(-(Nc-1),Nc):
  Gk=real(ifft(conj(fft(ge*roll(ge,-k)))*fft(ge)))
  Hk=real(ifft(conj(fft(ge))*fft(ge*roll(ge,k))))
  for j in range(-(Nc-1),Nc):
    A[k,j]=int(k==j)-(Gk[j]+Hk[j])/float(Y[k]*D)+Y[j]/float(D**2)
  

Ainv=inv(A)
Ce=dot(Ainv,Cp)
sC=0+0j
for k in range(-(Nc-1),Nc):
  sC+=(N-abs(k))*Ce[k]

sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2+real(sC)/float(N**2)
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N
mxe=sum(g.*x)/sum(g);
xe=[g.*(x-mxe),zeros(1,N)];
ge=[g,zeros(1,N)];
X=ifft(abs(fft(xe)).^2);
Y=real(ifft(abs(fft(ge)).^2));
Cp=zeros(1,2*Nc-1)+0j;
for k=-(Nc-1):Nc-1
Cp(mod(k,2*Nc-1)+1)=X(mod(k,2*N)+1)/Y(mod(k,2*N)+1);
end
D=sum(g);
A=zeros(2*Nc-1,2*Nc-1);
for k=-(Nc-1):Nc-1
Gk=real(ifft(conj(fft(ge.*circshift(ge,[0;-k]))).*fft(ge)));
Hk=real(ifft(conj(fft(ge)).*fft(ge.*circshift(ge,[0;k]))));
for j=-(Nc-1):Nc-1
A(mod(k,2*Nc-1)+1,mod(j,2*Nc-1)+1)=(k==j)-(Gk(mod(j,2*N)+1)+Hk(mod(j,2*N)+1))/(Y(mod(k,2*N)+1)*D)+Y(mod(j,2*N)+1)/(D^2);
end
end
Ainv=inv(A);
Ce=Ainv*Cp';
sC=0+0j;
for k=-(Nc-1):Nc-1
sC=sC+(N-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2+real(sC)/(N^2)
(Auto-)Korrelationsfunktion (erwartungstreu) Rk=R(τk)=∑i=max(0,−k)N−1−max(0,k)gigi+kxi*xi+k∑i=max(0,−k)N−1−max(0,k)gigi+k
τk=kΔt
k=−(N−1)…N−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
R=zeros(2*N-1)+0j
tau=zeros(2*N-1)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=g[i]*g[i+k]*conj(x[i])*x[i+k]
    sum0+=g[i]*g[i+k]
  
  R[k]=sum1/float(sum0)

plot(tau,real(R),'o',tau,imag(R),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=roll(arange(-(N-1),N)*dt,N)
X1=fft(append(g*x,zeros(N)))
E1=dt**2*abs(X1)**2
X0=fft(append(g,zeros(N)))
E0=dt**2*abs(X0)**2
R=ifft(E1)/real(ifft(E0))
R=append(R[0:N],R[N+1:2*N])
plot(tau,real(R),'o',tau,imag(R),'o')
show()
 direkte Berechnung:
R=zeros(1,2*N-1)+0j;
tau=zeros(1,2*N-1);
for k=-(N-1):N-1
tau(mod(k,2*N-1)+1)=k*dt;
sum1=0+0j;
sum0=0;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+g(i)*g(i+k)*conj(x(i))*x(i+k);
sum0=sum0+g(i)*g(i+k);
end
R(mod(k,2*N-1)+1)=sum1/sum0;
end
plot(tau,real(R),'o',tau,imag(R),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
tau=circshift((-(N-1):N-1)*dt,[0;N]);
X1=fft([g.*x,zeros(1,N)]);
E1=dt^2*abs(X1).^2;
X0=fft([g,zeros(1,N)]);
E0=dt^2*abs(X0).^2;
R=ifft(E1)./real(ifft(E0));
R=[R(1:N),R(N+2:end)];
plot(tau,real(R),'o',tau,imag(R),'o')
(Auto-)Kovarianzfunktion (asymptotisch erwartungstreu) Ck=C(τk)=∑i=max(0,−k)N−1−max(0,k)gigi+k(xi−x‾)*(xi+k−x‾)∑i=max(0,−k)N−1−max(0,k)gigi+k
τk=kΔt
k=−(N−1)…N−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
C=zeros(2*N-1)+0j
tau=zeros(2*N-1)
mxe=sum(g*x)/float(sum(g))
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=g[i]*g[i+k]*conj(x[i]-mxe)*(x[i+k]-mxe)
    sum0+=g[i]*g[i+k]
  
  C[k]=sum1/float(sum0)

plot(tau,real(C),'o',tau,imag(C),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=roll(arange(-(N-1),N)*dt,N)
mxe=sum(g*x)/float(sum(g))
X1=fft(append(g*(x-mxe),zeros(N)))
E1=dt**2*abs(X1)**2
X0=fft(append(g,zeros(N)))
E0=dt**2*abs(X0)**2
C=ifft(E1)/real(ifft(E0))
C=append(C[0:N],C[N+1:2*N])
plot(tau,real(C),'o',tau,imag(C),'o')
show()
 direkte Berechnung:
C=zeros(1,2*N-1)+0j;
tau=zeros(1,2*N-1);
mxe=sum(g.*x)/sum(g);
for k=-(N-1):N-1
tau(mod(k,2*N-1)+1)=k*dt;
sum1=0+0j;
sum0=0;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+g(i)*g(i+k)*conj(x(i)-mxe)*(x(i+k)-mxe);
sum0=sum0+g(i)*g(i+k);
end
C(mod(k,2*N-1)+1)=sum1/sum0;
end
plot(tau,real(C),'o',tau,imag(C),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
tau=circshift((-(N-1):N-1)*dt,[0;N]);
mxe=sum(g.*x)/sum(g);
X1=fft([g.*(x-mxe),zeros(1,N)]);
E1=dt^2*abs(X1).^2;
X0=fft([g,zeros(1,N)]);
E0=dt^2*abs(X0).^2;
C=ifft(E1)./real(ifft(E0));
C=[C(1:N),C(N+2:end)];
plot(tau,real(C),'o',tau,imag(C),'o')
(Auto-)Kovarianzfunktion (mit Bessel-Korrektur für korrelierte Daten, erwartungstreu) s. H. Nobach: Practical Realization of Bessel's Correction for a Bias-Free Estimation of the Auto-Covariance and the Cross-Covariance Functions direkte Berechnung:
from numpy.linalg import *
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N
Cp=zeros(2*Nc-1)+0j
S=zeros(2*Nc-1)
tau=zeros(2*Nc-1)
mxe=sum(g*x)/float(sum(g))
for k in range(-(Nc-1),Nc):
  tau[k]=k*dt
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=g[i]*g[i+k]*conj(x[i]-mxe)*(x[i+k]-mxe)
    sum0+=g[i]*g[i+k]
  
  Cp[k]=sum1/float(sum0)
  S[k]=sum0

D=sum(g)
A=zeros([2*Nc-1,2*Nc-1])
for k in range(-(Nc-1),Nc):
  for j in range(-(Nc-1),Nc):
    s=0
    for i in range(maximum(0,maximum(-j,-k)),minimum(N,minimum(N-j,N-k))):
      s+=g[i]*g[i+j]*g[i+k]
    for i in range(maximum(0,maximum(-j,k-j)),minimum(N,minimum(N-j,N+k-j))):
      s+=g[i]*g[i+j]*g[i+j-k]
    A[k,j]=int(k==j)-s/float(S[k]*D)+S[j]/float(D**2)
  

Ainv=inv(A)
C=dot(Ainv,Cp)
plot(tau,real(C),'o',tau,imag(C),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
from numpy.linalg import *
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N
tau=roll(arange(-(Nc-1),Nc)*dt,Nc)
mxe=sum(g*x)/float(sum(g))
xe=append(g*(x-mxe),zeros(N))
ge=append(g,zeros(N))
X=ifft(abs(fft(xe))**2)
Y=real(ifft(abs(fft(ge))**2))
Cp=zeros(2*Nc-1)+0j
for k in range(-(Nc-1),Nc):
  Cp[k]=X[k]/Y[k]

D=sum(g)
A=zeros([2*Nc-1,2*Nc-1])
for k in range(-(Nc-1),Nc):
  Gk=real(ifft(conj(fft(ge*roll(ge,-k)))*fft(ge)))
  Hk=real(ifft(conj(fft(ge))*fft(ge*roll(ge,k))))
  for j in range(-(Nc-1),Nc):
    A[k,j]=int(k==j)-(Gk[j]+Hk[j])/float(Y[k]*D)+Y[j]/float(D**2)
  

Ainv=inv(A)
C=dot(Ainv,Cp)
plot(tau,real(C),'o',tau,imag(C),'o')
show()
 direkte Berechnung:
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N
Cp=zeros(1,2*Nc-1)+0j;
S=zeros(1,2*Nc-1);
tau=zeros(1,2*Nc-1);
mxe=sum(g.*x)/sum(g);
for k=-(Nc-1):Nc-1
tau(mod(k,2*Nc-1)+1)=k*dt;
sum1=0+0j;
sum0=0;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+g(i)*g(i+k)*conj(x(i)-mxe)*(x(i+k)-mxe);
sum0=sum0+g(i)*g(i+k);
end
Cp(mod(k,2*Nc-1)+1)=sum1/sum0;
S(mod(k,2*Nc-1)+1)=sum0;
end
D=sum(g);
A=zeros(2*Nc-1,2*Nc-1);
for k=-(Nc-1):Nc-1
for j=-(Nc-1):Nc-1
s=0;
for i=max(0,max(-j,-k))+1:min(N,min(N-j,N-k))
s=s+g(i)*g(i+j)*g(i+k);
end
for i=max(0,max(-j,k-j))+1:min(N,min(N-j,N+k-j))
s=s+g(i)*g(i+j)*g(i+j-k);
end
A(mod(k,2*Nc-1)+1,mod(j,2*Nc-1)+1)=(k==j)-s/(S(mod(k,2*Nc-1)+1)*D)+S(mod(j,2*Nc-1)+1)/(D^2);
end
end
Ainv=inv(A);
C=Ainv*Cp';
plot(tau,real(C),'o',tau,imag(C),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N
tau=circshift((-(Nc-1):Nc-1)*dt,[0;Nc]);
mxe=sum(g.*x)/sum(g);
xe=[g.*(x-mxe),zeros(1,N)];
ge=[g,zeros(1,N)];
X=ifft(abs(fft(xe)).^2);
Y=real(ifft(abs(fft(ge)).^2));
Cp=zeros(1,2*Nc-1)+0j;
for k=-(Nc-1):Nc-1
Cp(mod(k,2*Nc-1)+1)=X(mod(k,2*N)+1)/Y(mod(k,2*N)+1);
end
D=sum(g);
A=zeros(2*Nc-1,2*Nc-1);
for k=-(Nc-1):Nc-1
Gk=real(ifft(conj(fft(ge.*circshift(ge,[0;-k]))).*fft(ge)));
Hk=real(ifft(conj(fft(ge)).*fft(ge.*circshift(ge,[0;k]))));
for j=-(Nc-1):Nc-1
A(mod(k,2*Nc-1)+1,mod(j,2*Nc-1)+1)=(k==j)-(Gk(mod(j,2*N)+1)+Hk(mod(j,2*N)+1))/(Y(mod(k,2*N)+1)*D)+Y(mod(j,2*N)+1)/(D^2);
end
end
Ainv=inv(A);
C=Ainv*Cp';
plot(tau,real(C),'o',tau,imag(C),'o')
(Auto-)Kovarianzfunktion (mit lokaler Normierung, asymptotisch erwartungstreu) Ck=C(τk)=s2⋅∑i=max(0,−k)N−1−max(0,k)gigi+k(xi−x‾)*(xi+k−x‾)[∑i=max(0,−k)N−1−max(0,k)gigi+k|xi−x‾|2][∑i=max(0,−k)N−1−max(0,k)gigi+k|xi+k−x‾|2]
τk=kΔt
k=−(N−1)…N−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
C=zeros(2*N-1)+0j
tau=zeros(2*N-1)
mxe=sum(g*x)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0+0j
  sum2=0
  sum3=0
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=g[i]*g[i+k]*conj(x[i]-mxe)*(x[i+k]-mxe)
    sum2+=g[i]*g[i+k]*abs(x[i]-mxe)**2
    sum3+=g[i]*g[i+k]*abs(x[i+k]-mxe)**2
  
  C[k]=vxe*sum1/sqrt(sum2*sum3)

plot(tau,real(C),'o',tau,imag(C),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=roll(arange(-(N-1),N)*dt,N)
mxe=sum(g*x)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
X=fft(append(g*(x-mxe),zeros(N)))
Q=fft(append(g*abs(x-mxe)**2,zeros(N)))
O=fft(append(g,zeros(N)))
EX=dt**2*abs(X)**2
EA=dt**2*conj(Q)*O
EB=dt**2*conj(O)*Q
CEX=ifft(EX)/float(dt)
CEA=real(ifft(EA))/float(dt)
CEB=real(ifft(EB))/float(dt)
C=vxe*CEX/sqrt(CEA*CEB)
C=append(C[0:N],C[N+1:2*N])
plot(tau,real(C),'o',tau,imag(C),'o')
show()
 direkte Berechnung:
C=zeros(1,2*N-1)+0j;
tau=zeros(1,2*N-1);
mxe=sum(g.*x)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
for k=-(N-1):N-1
tau(mod(k,2*N-1)+1)=k*dt;
sum1=0+0j;
sum2=0;
sum3=0;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+g(i)*g(i+k)*conj(x(i)-mxe)*(x(i+k)-mxe);
sum2=sum2+g(i)*g(i+k)*abs(x(i)-mxe)^2;
sum3=sum3+g(i)*g(i+k)*abs(x(i+k)-mxe)^2;
end
C(mod(k,2*N-1)+1)=vxe*sum1/sqrt(sum2*sum3);
end
plot(tau,real(C),'o',tau,imag(C),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
tau=circshift((-(N-1):N-1)*dt,[0;N]);
mxe=sum(g.*x)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
X=fft([g.*(x-mxe),zeros(1,N)]);
Q=fft([g.*abs(x-mxe).^2,zeros(1,N)]);
O=fft([g,zeros(1,N)]);
EX=dt^2*abs(X).^2;
EA=dt^2*conj(Q).*O;
EB=dt^2*conj(O).*Q;
CEX=ifft(EX)/dt;
CEA=real(ifft(EA))/dt;
CEB=real(ifft(EB))/dt;
C=vxe*CEX./sqrt(CEA.*CEB);
C=[C(1:N),C(N+2:end)];
plot(tau,real(C),'o',tau,imag(C),'o')
(Auto-)Korrelationskoeffizientenfunktion (asymptotisch erwartungstreu) ρk=ρ(τk)=∑i=max(0,−k)N−1−max(0,k)gigi+k(xi−x‾)*(xi+k−x‾)s2⋅∑i=max(0,−k)N−1−max(0,k)gigi+k
τk=kΔt
k=−(N−1)…N−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
rho=zeros(2*N-1)+0j
tau=zeros(2*N-1)
mxe=sum(g*x)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=g[i]*g[i+k]*conj(x[i]-mxe)*(x[i+k]-mxe)
    sum0+=g[i]*g[i+k]
  
  rho[k]=sum1/sum0/vxe

plot(tau,real(rho),'o',tau,imag(rho),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=roll(arange(-(N-1),N)*dt,N)
mxe=sum(g*x)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
X1=fft(append(g*(x-mxe),zeros(N)))
E1=dt**2*abs(X1)**2
X0=fft(append(g,zeros(N)))
E0=dt**2*abs(X0)**2
rho=ifft(E1)/real(ifft(E0))/vxe
rho=append(rho[0:N],rho[N+1:2*N])
plot(tau,real(rho),'o',tau,imag(rho),'o')
show()
 direkte Berechnung:
rho=zeros(1,2*N-1)+0j;
tau=zeros(1,2*N-1);
mxe=sum(g.*x)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
for k=-(N-1):N-1
tau(mod(k,2*N-1)+1)=k*dt;
sum1=0+0j;
sum0=0;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+g(i)*g(i+k)*conj(x(i)-mxe)*(x(i+k)-mxe);
sum0=sum0+g(i)*g(i+k);
end
rho(mod(k,2*N-1)+1)=sum1/sum0/vxe;
end
plot(tau,real(rho),'o',tau,imag(rho),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
tau=circshift((-(N-1):N-1)*dt,[0;N]);
mxe=sum(g.*x)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
X1=fft([g.*(x-mxe),zeros(1,N)]);
E1=dt^2*abs(X1).^2;
X0=fft([g,zeros(1,N)]);
E0=dt^2*abs(X0).^2;
rho=ifft(E1)./real(ifft(E0))/vxe;
rho=[rho(1:N),rho(N+2:end)];
plot(tau,real(rho),'o',tau,imag(rho),'o')
(Auto-)Korrelationskoeffizientenfunktion (mit lokaler Normierung, asymptotisch erwartungstreu) ρk=ρ(τk)=∑i=max(0,−k)N−1−max(0,k)gigi+k(xi−x‾)*(xi+k−x‾)[∑i=max(0,−k)N−1−max(0,k)gigi+k|xi−x‾|2][∑i=max(0,−k)N−1−max(0,k)gigi+k|xi+k−x‾|2]
τk=kΔt
k=−(N−1)…N−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
rho=zeros(2*N-1)+0j
tau=zeros(2*N-1)
mxe=sum(g*x)/float(sum(g))
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0+0j
  sum2=0
  sum3=0
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=g[i]*g[i+k]*conj(x[i]-mxe)*(x[i+k]-mxe)
    sum2+=g[i]*g[i+k]*abs(x[i]-mxe)**2
    sum3+=g[i]*g[i+k]*abs(x[i+k]-mxe)**2
  
  rho[k]=sum1/sqrt(sum2*sum3)

plot(tau,real(rho),'o',tau,imag(rho),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=roll(arange(-(N-1),N)*dt,N)
mxe=sum(g*x)/float(sum(g))
X=fft(append(g*(x-mxe),zeros(N)))
Q=fft(append(g*abs(x-mxe)**2,zeros(N)))
O=fft(append(g,zeros(N)))
EX=dt**2*abs(X)**2
EA=dt**2*conj(Q)*O
EB=dt**2*conj(O)*Q
CEX=ifft(EX)/float(dt)
CEA=real(ifft(EA))/float(dt)
CEB=real(ifft(EB))/float(dt)
rho=CEX/sqrt(CEA*CEB)
rho=append(rho[0:N],rho[N+1:2*N])
plot(tau,real(rho),'o',tau,imag(rho),'o')
show()
 direkte Berechnung:
rho=zeros(1,2*N-1)+0j;
tau=zeros(1,2*N-1);
mxe=sum(g.*x)/sum(g);
for k=-(N-1):N-1
tau(mod(k,2*N-1)+1)=k*dt;
sum1=0+0j;
sum2=0;
sum3=0;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+g(i)*g(i+k)*conj(x(i)-mxe)*(x(i+k)-mxe);
sum2=sum2+g(i)*g(i+k)*abs(x(i)-mxe)^2;
sum3=sum3+g(i)*g(i+k)*abs(x(i+k)-mxe)^2;
end
rho(mod(k,2*N-1)+1)=sum1/sqrt(sum2*sum3);
end
plot(tau,real(rho),'o',tau,imag(rho),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
tau=circshift((-(N-1):N-1)*dt,[0;N]);
mxe=sum(g.*x)/sum(g);
X=fft([g.*(x-mxe),zeros(1,N)]);
Q=fft([g.*abs(x-mxe).^2,zeros(1,N)]);
O=fft([g,zeros(1,N)]);
EX=dt^2*abs(X).^2;
EA=dt^2*conj(Q).*O;
EB=dt^2*conj(O).*Q;
CEX=ifft(EX)/dt;
CEA=real(ifft(EA))/dt;
CEB=real(ifft(EB))/dt;
rho=CEX./sqrt(CEA.*CEB);
rho=[rho(1:N),rho(N+2:end)];
plot(tau,real(rho),'o',tau,imag(rho),'o')
(Auto-)Leistungsdichtespektrum (erwartungstreu) (wegen der nötigen Normierung nur über die Korrelationsfunktion, evtl. Einschränkung des verwendeten Bereichs der Korrelationsfunktion zur Reduktion der Schätzvarianz des Leistungsdichtespektrums, K≤2N−1)
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⌋ 		Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren:
from numpy.fft import *
X1=fft(append(g*x,zeros(N)))
E1=dt**2*abs(X1)**2
X0=fft(append(g,zeros(N)))
E0=dt**2*abs(X0)**2
RA=ifft(E1)/real(ifft(E0))
K=100
K=minimum(2*N-1,K)
R=zeros(K)+0j
for k in range(-(K//2),(K-1)//2+1):
  R[k]=RA[k]

f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
plot(f,S,'o')
show()
Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
K=100
K=minimum(2*N-1,K)
R=zeros(K)+0j
for k in range(-(K//2),(K-1)//2+1):
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=g[i]*g[i+k]*conj(x[i])*x[i+k]
    sum0+=g[i]*g[i+k]
  
  R[k]=sum1/float(sum0)

f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
plot(f,S,'o')
show()
 Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren:
X1=fft([g.*x,zeros(1,N)]);
E1=dt^2*abs(X1).^2;
X0=fft([g,zeros(1,N)]);
E0=dt^2*abs(X0).^2;
RA=ifft(E1)./real(ifft(E0));
K=100;
K=min(2*N-1,K);
R=zeros(1,K)+0j;
for k=-floor(K/2):floor((K-1)/2)
R(mod(k,K)+1)=RA(mod(k,2*N)+1);
end
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
plot(f,S,'o')
Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
K=100;
K=min(2*N-1,K);
R=zeros(1,K)+0j;
for k=-floor(K/2):floor((K-1)/2)
sum1=0+0j;
sum0=0;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+g(i)*g(i+k)*conj(x(i))*x(i+k);
sum0=sum0+g(i)*g(i+k);
end
R(mod(k,K)+1)=sum1/sum0;
end
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
plot(f,S,'o')