Signal- und Messdatenverarbeitung

Codes für komplexe, stochastische Leistungssignalpaare, ohne Gewichtung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren von Nx komplexen Werten xi mit und Ny komplexen Werten yi mit i=0…Ny−1 (korrelierte und gegeneinander verschobene AR1-Prozesse als Beispiel, mit Erwartungswerten verschieden von null)
Nx=1000
Ny=1100
Nd=25
dt=1.0
tx=arange(0,Nx)*dt
ty=arange(0,Ny)*dt
phi=0.95
theta=sqrt(1-phi**2)
mxs=3+1.5j
mys=2+0.5j
vxs=1
vys=2
cc=0.5
c1=sqrt(0.5+0.5*sqrt(1-cc**2/(vxs*vys)))
c2=cc/sqrt(vxs*vys)/(2*c1)
c11=sqrt(vxs)*c1
c12=sqrt(vxs)*c2
c21=sqrt(vys)*c2
c22=sqrt(vys)*c1
Nuv=maximum(maximum(Nx,Nx+Nd),maximum(Ny,Ny-Nd))
u=zeros(Nuv)+0j
v=zeros(Nuv)+0j
x=zeros(Nx)+0j
y=zeros(Ny)+0j
for i in range(0,Nuv):
  if (i==0):
    u[i]=normal(0,sqrt(0.5))+1j*normal(0,sqrt(0.5))
    v[i]=normal(0,sqrt(0.5))+1j*normal(0,sqrt(0.5))
  
  else:
    u[i]=phi*u[i-1]+normal(0,sqrt(0.5)*theta)+1j*normal(0,sqrt(0.5)*theta)
    v[i]=phi*v[i-1]+normal(0,sqrt(0.5)*theta)+1j*normal(0,sqrt(0.5)*theta)
  
  if (i-maximum(0,Nd)>=0) and (i-maximum(0,Nd)<Nx):
    x[i-maximum(0,Nd)]=c11*u[i]+c12*v[i]+mxs
  
  if (i-maximum(0,-Nd)>=0) and (i-maximum(0,-Nd)<Ny):
    y[i-maximum(0,-Nd)]=c21*u[i]+c22*v[i]+mys
  

plot(tx,real(x),'o',tx,imag(x),'o',ty,real(y),'o',ty,imag(y),'o')
show()
Nx=1000;
Ny=1100;
Nd=25;
dt=1.0;
tx=(0:Nx-1)*dt;
ty=(0:Ny-1)*dt;
phi=0.95;
theta=sqrt(1-phi^2);
mxs=3+1.5j;
mys=2+0.5j;
vxs=1;
vys=2;
cc=0.5;
c1=sqrt(0.5+0.5*sqrt(1-cc^2/(vxs*vys)));
c2=cc/sqrt(vxs*vys)/(2*c1);
c11=sqrt(vxs)*c1;
c12=sqrt(vxs)*c2;
c21=sqrt(vys)*c2;
c22=sqrt(vys)*c1;
Nuv=max(max(Nx,Nx+Nd),max(Ny,Ny-Nd));
u=zeros(1,Nuv)+0j;
v=zeros(1,Nuv)+0j;
x=zeros(1,Nx)+0j;
y=zeros(1,Ny)+0j;
for i=1:Nuv
if (i==1)
u(i)=sqrt(0.5)*randn()+1j*sqrt(0.5)*randn();
v(i)=sqrt(0.5)*randn()+1j*sqrt(0.5)*randn();
else
u(i)=phi*u(i-1)+sqrt(0.5)*theta*randn()+1j*sqrt(0.5)*theta*randn();
v(i)=phi*v(i-1)+sqrt(0.5)*theta*randn()+1j*sqrt(0.5)*theta*randn();
end
if (i-max(0,Nd)>0) && (i-max(0,Nd)<=Nx)
x(i-max(0,Nd))=c11*u(i)+c12*v(i)+mxs;
end
if (i-max(0,-Nd)>0) && (i-max(0,-Nd)<=Ny)
y(i-max(0,-Nd))=c21*u(i)+c22*v(i)+mys;
end
end
plot(tx,real(x),'o',tx,imag(x),'o',ty,real(y),'o',ty,imag(y),'o')
Mittelwerte x‾=1Nx⋅∑i=0Nx−1xi
y‾=1Ny⋅∑i=0Ny−1yi
mxe=mean(x)
mye=mean(y)
mxe=mean(x);
mye=mean(y);
Varianzen (ohne Bessel-Korrektur, asymptotisch erwartungstreu) sx2=1Nx⋅∑i=0Nx−1(xi−x‾)*(xi−x‾)=1Nx⋅∑i=0Nx−1|xi−x‾|2
=(1Nx⋅∑i=0Nx−1xi*xi)−x‾*x‾=(1Nx⋅∑i=0Nx−1|xi|2)−|x‾|2
sy2=1Ny⋅∑i=0Ny−1(yi−y‾)*(yi−y‾)=1Ny⋅∑i=0Ny−1|yi−y‾|2
=(1Ny⋅∑i=0Ny−1yi*yi)−y‾*y‾=(1Ny⋅∑i=0Ny−1|yi|2)−|y‾|2
vxe=var(x)
vye=var(y)
vxe=var(x,1);
vye=var(y,1);
Varianzen (mit Bessel-Korrektur für korrelierte Datenreihen xi und yi, erwartungstreu) sx2=Cxx(0)
sy2=Cyy(0)
mit den erwartungstreuen Schätzungen Cxx(τk) und Cyy(τk) der (Auto-)Kovarianzfunktionen der beiden Zeitreihen xi und yi mit Bessel-Korrektur für korrelierte Daten (s. hier)
from numpy.fft import *
from numpy.linalg import *
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Cp=zeros(2*Nc-1)+0j
mxe=mean(x)
X=fft(append(x-mxe,zeros(Nx)))
E=dt**2*abs(X)**2
CE=ifft(E)/float(dt)
for k in range(-(Nc-1),Nc):
  Cp[k]=CE[k]/float(dt)/float(Nx-abs(k))

A=zeros([2*Nc-1,2*Nc-1])
for i in range(-(Nc-1),Nc):
  for j in range(-(Nc-1),Nc):
    A[i,j]=int(i==j)+(Nx-abs(j))/float(Nx**2)-2*(Nx-maximum(abs(i),maximum(abs(j),minimum(Nx,abs(i-j)))))/float(Nx*(Nx-abs(i)))
  

Ainv=inv(A)
Ce=dot(Ainv,Cp)
vxe=real(Ce[0])
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cp=zeros(2*Nc-1)+0j
mye=mean(y)
Y=fft(append(y-mye,zeros(Ny)))
E=dt**2*abs(Y)**2
CE=ifft(E)/float(dt)
for k in range(-(Nc-1),Nc):
  Cp[k]=CE[k]/float(dt)/float(Ny-abs(k))

A=zeros([2*Nc-1,2*Nc-1])
for i in range(-(Nc-1),Nc):
  for j in range(-(Nc-1),Nc):
    A[i,j]=int(i==j)+(Ny-abs(j))/float(Ny**2)-2*(Ny-maximum(abs(i),maximum(abs(j),minimum(Ny,abs(i-j)))))/float(Ny*(Ny-abs(i)))
  

Ainv=inv(A)
Ce=dot(Ainv,Cp)
vye=real(Ce[0])
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Cp=zeros(1,2*Nc-1)+0j;
mxe=mean(x);
X=fft([x-mxe,zeros(1,Nx)]);
E=dt^2*abs(X).^2;
CE=ifft(E)/dt;
for k=-(Nc-1):Nc-1
Cp(mod(k,2*Nc-1)+1)=CE(mod(k,2*Nx)+1)/dt/(Nx-abs(k));
end
A=zeros(2*Nc-1,2*Nc-1);
for i=-(Nc-1):Nc-1
for j=-(Nc-1):Nc-1
A(mod(i,2*Nc-1)+1,mod(j,2*Nc-1)+1)=(i==j)+(Nx-abs(j))/Nx^2-2*(Nx-max(abs(i),abs(j)),min(Nx,abs(i-j))))/(Nx*(Nx-abs(i)));
end
end
Ainv=inv(A);
Ce=Ainv*Cp';
vxe=real(Ce(1));
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cp=zeros(1,2*Nc-1)+0j;
mye=mean(y);
Y=fft([y-mye,zeros(1,Ny)]);
E=dt^2*abs(Y).^2;
CE=ifft(E)/dt;
for k=-(Nc-1):Nc-1
Cp(mod(k,2*Nc-1)+1)=CE(mod(k,2*Ny)+1)/dt/(Ny-abs(k));
end
A=zeros(2*Nc-1,2*Nc-1);
for i=-(Nc-1):Nc-1
for j=-(Nc-1):Nc-1
A(mod(i,2*Nc-1)+1,mod(j,2*Nc-1)+1)=(i==j)+(Ny-abs(j))/Ny^2-2*(Ny-max(max(abs(i),abs(j)),min(Ny,abs(i-j))))/(Ny*(Ny-abs(i)));
end
end
Ainv=inv(A);
Ce=Ainv*Cp';
vye=real(Ce(1));
Kreuzkovarianz (ohne Bessel-Korrektur, asymptotisch erwartungstreu) cyx=cxy*=1min(Nx,Ny)⋅∑i=0min(Nx,Ny)−1(xi−x‾)*(yi−y‾)
cov(y[0:minimum(Nx,Ny)],x[0:minimum(Nx,Ny)],ddof=0)[0,1]
sum(conj(x(1:min(Nx,Ny))-mean(x)).*(y(1:min(Nx,Ny))-mean(y)))/min(Nx,Ny)
Kreuzkovarianz (mit Bessel-Korrektur für korrelierte Datenreihen xi und yi, erwartungstreu) cyx=Cyx(0)
mit der erwartungstreuen Schätzung Cyx(τk) der Kreuzkovarianzfunktion mit Bessel-Korrektur für korrelierte Daten (s. unten)
from numpy.fft import *
from numpy.linalg import *
Nxc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Nyc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cyxp=zeros(Nxc+Nyc-1)+0j
mxe=mean(x)
mye=mean(y)
X=fft(append(x-mxe,zeros(Ny)))
Y=fft(append(y-mye,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
CEyx=ifft(Eyx)/float(dt)
for k in range(-(Nxc-1),Nyc):
  Cyxp[k]=CEyx[k]/float(dt)/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

A=zeros([Nxc+Nyc-1,Nxc+Nyc-1])
for i in range(-(Nxc-1),Nyc):
  for j in range(-(Nxc-1),Nyc):
    A[i,j]=int(i==j)-(minimum(Nx,minimum(Ny-i,Ny-j))-maximum(0,maximum(-i,-j)))/float(Ny*(minimum(Nx,Ny-i)-maximum(0,-i)))-(minimum(Nx,Ny-j,maximum(0,Nx+i-j))-maximum(0,maximum(-j,minimum(Nx,i-j))))/float(Nx*(minimum(Nx,Ny-i)-maximum(0,-i)))+(minimum(Nx,Ny-j)-maximum(0,-j))/float(Nx*Ny)
  

Ainv=inv(A)
Cyxe=dot(Ainv,Cyxp)
Cyxe[0]
Nxc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Nyc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cyxp=zeros(1,Nxc+Nyc-1)+0j;
mxe=mean(x);
mye=mean(y);
X=fft([x-mxe,zeros(1,Ny)]);
Y=fft([y-mye,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
CEyx=ifft(Eyx)/dt;
for k=-(Nxc-1):Nyc-1
Cyxp(mod(k,Nxc+Nyc-1)+1)=CEyx(mod(k,Nx+Ny)+1)/dt/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
A=zeros(Nxc+Nyc-1,Nxc+Nyc-1);
for i=-(Nxc-1):Nyc-1
for j=-(Nxc-1):Nyc-1
A(mod(i,Nxc+Nyc-1)+1,mod(j,Nxc+Nyc-1)+1)=(i==j)-(min(min(Nx,Ny-i),Ny-j)-max(max(0,-i),-j))/(Ny*(min(Nx,Ny-i)-max(0,-i)))-(min(min(Nx,Ny-j),max(0,Nx+i-j))-max(max(0,-j),min(Nx,i-j)))/(Nx*(min(Nx,Ny-i)-max(0,-i)))+(min(Nx,Ny-j)-max(0,-j))/(Nx*Ny);
end
end
Ainv=inv(A);
Cyxe=Ainv*Cyxp';
Cyxe(1)
(Kreuz-)Korrelationskoeffizient (asymptotisch erwartungstreu) γyx=γxy*=[∑i=0min(Nx,Ny)−1(xi−x‾)*(yi−y‾)]⋅NxNymin(Nx,Ny)⋅(∑i=0Nx−1|xi−x‾|2)(∑i=0Ny−1|yi−y‾|2)
=[∑i=0min(Nx,Ny)−1(xi−x‾)*(yi−y‾)]⋅NxNymin(Nx,Ny)⋅[(∑i=0Nx−1|xi|2)−Nx|x‾|2][(∑i=0Ny−1|yi|2)−Ny|y‾|2]
corrcoef(y[0:minimum(Nx,Ny)],x[0:minimum(Nx,Ny)],ddof=0)[0,1]
sum(conj(x(1:min(Nx,Ny))-mean(x)).*(y(1:min(Nx,Ny))-mean(y)))/(min(Nx,Ny)*sqrt(var(x,1)*var(y,1)))
(Kreuz-)Korrelationskoeffizient (mit lokaler Normierung, asymptotisch erwartungstreu) γyx=γxy*=∑i=0min(Nx,Ny)−1(xi−x‾)*(yi−y‾)[∑i=0min(Nx,Ny)−1|xi−x‾|2][∑i=0min(Nx,Ny)−1|yi−y‾|2]
mxe=mean(x)
mye=mean(y)
sum(conj(x[0:minimum(Nx,Ny)]-mxe)*(y[0:minimum(Nx,Ny)]-mye))/sqrt(sum(abs(x[0:minimum(Nx,Ny)]-mxe)**2)*sum(abs(y[0:minimum(Nx,Ny)]-mye)**2))
mxe=mean(x);
mye=mean(y);
sum(conj(x(1:min(Nx,Ny))-mxe).*(y(1:min(Nx,Ny))-mye))/sqrt(sum(abs(x(1:min(Nx,Ny))-mxe).^2)*sum(abs(y(1:min(Nx,Ny))-mye).^2))
Kreuzkorrelationsfunktion (erwartungstreu) Ryx,k=Ryx(τk)=Rxy*(−τk)=1min(Nx,Nx+k,Ny,Ny−k)⋅∑i=max(0,−k)min(Nx,Ny−k)−1xi*yi+k
τk=kΔt
k=−(Nx−1)…Ny−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
Ryx=zeros(Nx+Ny-1)+0j
tau=zeros(Nx+Ny-1)
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i])*y[i+k]
  
  Ryx[k]=sum1/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
Ryx=zeros(Nx+Ny-1)+0j
tau=roll(arange(-(Nx-1),Ny)*dt,Ny)
X=fft(append(x,zeros(Ny)))
Y=fft(append(y,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
REyx=ifft(Eyx)/float(dt)
for k in range(-(Nx-1),Ny):
  Ryx[k]=REyx[k]/float(dt)/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
show()
 direkte Berechnung:
Ryx=zeros(1,Nx+Ny-1)+0j;
tau=zeros(1,Nx+Ny-1);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+conj(x(i))*y(i+k);
end
Ryx(mod(k,Nx+Ny-1)+1)=sum1/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Ryx=zeros(1,Nx+Ny-1)+0j;
tau=circshift((-(Nx-1):Ny-1)*dt,[0;Ny]);
X=fft([x,zeros(1,Ny)]);
Y=fft([y,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
REyx=ifft(Eyx)/dt;
for k=-(Nx-1):Ny-1
Ryx(mod(k,Nx+Ny-1)+1)=REyx(mod(k,Nx+Ny)+1)/dt/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Kreuzkovarianzfunktion (ohne Bessel-Korrektur, asymptotisch erwartungstreu) Cyx,k=Cyx(τk)=Cxy*(−τk)=1min(Nx,Nx+k,Ny,Ny−k)⋅∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)
τk=kΔt
k=−(Nx−1)…Ny−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
Cyx=zeros(Nx+Ny-1)+0j
tau=zeros(Nx+Ny-1)
mxe=mean(x)
mye=mean(y)
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i]-mxe)*(y[i+k]-mye)
  
  Cyx[k]=sum1/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
Cyx=zeros(Nx+Ny-1)+0j
tau=roll(arange(-(Nx-1),Ny)*dt,Ny)
mxe=mean(x)
mye=mean(y)
X=fft(append(x-mxe,zeros(Ny)))
Y=fft(append(y-mye,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
CEyx=ifft(Eyx)/float(dt)
for k in range(-(Nx-1),Ny):
  Cyx[k]=CEyx[k]/float(dt)/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
 direkte Berechnung:
Cyx=zeros(1,Nx+Ny-1)+0j;
tau=zeros(1,Nx+Ny-1);
mxe=mean(x);
mye=mean(y);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+conj(x(i)-mxe)*(y(i+k)-mye);
end
Cyx(mod(k,Nx+Ny-1)+1)=sum1/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Cyx=zeros(1,Nx+Ny-1)+0j;
tau=circshift((-(Nx-1):Ny-1)*dt,[0;Ny]);
mxe=mean(x);
mye=mean(y);
X=fft([x-mxe,zeros(1,Ny)]);
Y=fft([y-mye,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
CEyx=ifft(Eyx)/dt;
for k=-(Nx-1):Ny-1
Cyx(mod(k,Nx+Ny-1)+1)=CEyx(mod(k,Nx+Ny)+1)/dt/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Kreuzkovarianzfunktion (mit Bessel-Korrektur für korrelierte Datenreihen xi und yi, 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 *
Nxc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Nyc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cyxp=zeros(Nxc+Nyc-1)+0j
tau=zeros(Nxc+Nyc-1)
mxe=mean(x)
mye=mean(y)
for k in range(-(Nxc-1),Nyc):
  tau[k]=k*dt
  sum1=0+0j
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i]-mxe)*(y[i+k]-mye)
  
  Cyxp[k]=sum1/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

A=zeros([Nxc+Nyc-1,Nxc+Nyc-1])
for i in range(-(Nxc-1),Nyc):
  for j in range(-(Nxc-1),Nyc):
    A[i,j]=int(i==j)-(minimum(Nx,minimum(Ny-i,Ny-j))-maximum(0,maximum(-i,-j)))/float(Ny*(minimum(Nx,Ny-i)-maximum(0,-i)))-(minimum(Nx,Ny-j,maximum(0,Nx+i-j))-maximum(0,maximum(-j,minimum(Nx,i-j))))/float(Nx*(minimum(Nx,Ny-i)-maximum(0,-i)))+(minimum(Nx,Ny-j)-maximum(0,-j))/float(Nx*Ny)
  

Ainv=inv(A)
Cyx=dot(Ainv,Cyxp)
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
from numpy.linalg import *
Nxc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Nyc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cyxp=zeros(Nxc+Nyc-1)+0j
tau=roll(arange(-(Nxc-1),Nyc)*dt,Nyc)
mxe=mean(x)
mye=mean(y)
X=fft(append(x-mxe,zeros(Ny)))
Y=fft(append(y-mye,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
CEyx=ifft(Eyx)/float(dt)
for k in range(-(Nxc-1),Nyc):
  Cyxp[k]=CEyx[k]/float(dt)/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

A=zeros([Nxc+Nyc-1,Nxc+Nyc-1])
for i in range(-(Nxc-1),Nyc):
  for j in range(-(Nxc-1),Nyc):
    A[i,j]=int(i==j)-(minimum(Nx,minimum(Ny-i,Ny-j))-maximum(0,maximum(-i,-j)))/float(Ny*(minimum(Nx,Ny-i)-maximum(0,-i)))-(minimum(Nx,Ny-j,maximum(0,Nx+i-j))-maximum(0,maximum(-j,minimum(Nx,i-j))))/float(Nx*(minimum(Nx,Ny-i)-maximum(0,-i)))+(minimum(Nx,Ny-j)-maximum(0,-j))/float(Nx*Ny)
  

Ainv=inv(A)
Cyx=dot(Ainv,Cyxp)
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
 direkte Berechnung:
Nxc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Nyc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cyxp=zeros(1,Nxc+Nyc-1)+0j;
tau=zeros(1,Nxc+Nyc-1);
mxe=mean(x);
mye=mean(y);
for k=-(Nxc-1):Nyc-1
tau(mod(k,Nxc+Nyc-1)+1)=k*dt;
sum1=0+0j;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+conj(x(i)-mxe)*(y(i+k)-mye);
end
Cyxp(mod(k,Nxc+Nyc-1)+1)=sum1/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
A=zeros(Nxc+Nyc-1,Nxc+Nyc-1);
for i=-(Nxc-1):Nyc-1
for j=-(Nxc-1):Nyc-1
A(mod(i,Nxc+Nyc-1)+1,mod(j,Nxc+Nyc-1)+1)=(i==j)-(min(min(Nx,Ny-i),Ny-j)-max(max(0,-i),-j))/(Ny*(min(Nx,Ny-i)-max(0,-i)))-(min(min(Nx,Ny-j),max(0,Nx+i-j))-max(max(0,-j),min(Nx,i-j)))/(Nx*(min(Nx,Ny-i)-max(0,-i)))+(min(Nx,Ny-j)-max(0,-j))/(Nx*Ny);
end
end
Ainv=inv(A);
Cyx=Ainv*Cyxp';
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Nxc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
Nyc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
Cyxp=zeros(1,Nxc+Nyc-1)+0j;
tau=circshift((-(Nxc-1):Nyc-1)*dt,[0;Nyc]);
mxe=mean(x);
mye=mean(y);
X=fft([x-mxe,zeros(1,Ny)]);
Y=fft([y-mye,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
CEyx=ifft(Eyx)/dt;
for k=-(Nxc-1):Nyc-1
Cyxp(mod(k,Nxc+Nyc-1)+1)=CEyx(mod(k,Nx+Ny)+1)/dt/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
A=zeros(Nxc+Nyc-1,Nxc+Nyc-1);
for i=-(Nxc-1):Nyc-1
for j=-(Nxc-1):Nyc-1
A(mod(i,Nxc+Nyc-1)+1,mod(j,Nxc+Nyc-1)+1)=(i==j)-(min(min(Nx,Ny-i),Ny-j)-max(max(0,-i),-j))/(Ny*(min(Nx,Ny-i)-max(0,-i)))-(min(min(Nx,Ny-j),max(0,Nx+i-j))-max(max(0,-j),min(Nx,i-j)))/(Nx*(min(Nx,Ny-i)-max(0,-i)))+(min(Nx,Ny-j)-max(0,-j))/(Nx*Ny);
end
end
Ainv=inv(A);
Cyx=Ainv*Cyxp';
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Kreuzkovarianzfunktion (mit lokaler Normierung, asymptotisch erwartungstreu) Cyx,k=Cyx(τk)=Cxy*(−τk)=[∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)]⋅sx2sy2[∑i=max(0,−k)min(Nx,Ny−k)−1|xi−x‾|2][∑i=max(0,k)min(Nx+k,Ny)−1|yi−y‾|2]
=[∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)]⋅(∑i=0Nx−1|xi−x‾|2)(∑i=0Ny−1|yi−y‾|2)NxNy[∑i=max(0,−k)min(Nx,Ny−k)−1|xi−x‾|2][∑i=max(0,k)min(Nx+k,Ny)−1|yi−y‾|2]
=[∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)]⋅[(∑i=0Nx−1|xi|2)−Nx|x‾|2][(∑i=0Ny−1|yi|2)−Ny|y‾|2]NxNy[∑i=max(0,−k)min(Nx,Ny−k)−1|xi−x‾|2][∑i=max(0,k)min(Nx+k,Ny)−1|yi−y‾|2]
τk=kΔt
k=−(Nx−1)…Ny−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
Cyx=zeros(Nx+Ny-1)+0j
tau=zeros(Nx+Ny-1)
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  sum2=0+0j
  sum3=0+0j
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i]-mxe)*(y[i+k]-mye)
    sum2+=abs(x[i]-mxe)**2
    sum3+=abs(y[i+k]-mye)**2
  
  Cyx[k]=sqrt(vxe*vye)*sum1/sqrt(sum2*sum3)

plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
Cyx=zeros(Nx+Ny-1)+0j
tau=roll(arange(-(Nx-1),Ny)*dt,Ny)
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
X=fft(append(x-mxe,zeros(Ny)))
Y=fft(append(y-mye,zeros(Nx)))
QX=fft(append(abs(x-mxe)**2,zeros(Ny)))
QY=fft(append(abs(y-mye)**2,zeros(Nx)))
OX=fft(append(ones(Nx),zeros(Ny)))
OY=fft(append(ones(Ny),zeros(Nx)))
Eyx=dt**2*conj(X)*Y
EA=dt**2*conj(QX)*OY
EB=dt**2*conj(OX)*QY
CEyx=ifft(Eyx)/float(dt)
CEA=ifft(EA)/float(dt)
CEB=ifft(EB)/float(dt)
for k in range(-(Nx-1),Ny):
  Cyx[k]=sqrt(vxe*vye)*CEyx[k]/sqrt(CEA[k]*CEB[k])

plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
 direkte Berechnung:
Cyx=zeros(1,Nx+Ny-1)+0j;
tau=zeros(1,Nx+Ny-1);
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
sum2=0+0j;
sum3=0+0j;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+conj(x(i)-mxe)*(y(i+k)-mye);
sum2=sum2+abs(x(i)-mxe)^2;
sum3=sum3+abs(y(i+k)-mye)^2;
end
Cyx(mod(k,Nx+Ny-1)+1)=sqrt(vxe*vye)*sum1/sqrt(sum2*sum3);
end
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Cyx=zeros(1,Nx+Ny-1)+0j;
tau=circshift((-(Nx-1):Ny-1)*dt,[0;Ny]);
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
X=fft([x-mxe,zeros(1,Ny)]);
Y=fft([y-mye,zeros(1,Nx)]);
QX=fft([abs(x-mxe).^2,zeros(1,Ny)]);
QY=fft([abs(y-mye).^2,zeros(1,Nx)]);
OX=fft([ones(1,Nx),zeros(1,Ny)]);
OY=fft([ones(1,Ny),zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
EA=dt^2*conj(QX).*OY;
EB=dt^2*conj(OX).*QY;
CEyx=ifft(Eyx)/dt;
CEA=ifft(EA)/dt;
CEB=ifft(EB)/dt;
for k=-(Nx-1):Ny-1
Cyx(mod(k,Nx+Ny-1)+1)=sqrt(vxe*vye)*CEyx(mod(k,Nx+Ny)+1)/dt/sqrt(CEA(mod(k,Nx+Ny)+1)*CEB(mod(k,Nx+Ny)+1));
end
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Kreuzkorrelationskoeffizientenfunktion (asymptotisch erwartungstreu) ρyx,k=ρyx(τk)=ρxy*(−τk)=1min(Nx,Nx+k,Ny,Ny−k)⋅sx2sy2⋅∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)
=[∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)]⋅NxNymin(Nx,Nx+k,Ny,Ny−k)⋅(∑i=0Nx−1|xi−x‾|2)(∑i=0Ny−1|yi−y‾|2)
=[∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)]⋅NxNymin(Nx,Nx+k,Ny,Ny−k)⋅[(∑i=0Nx−1|xi|2)−Nx|x‾|2][(∑i=0Ny−1|yi|2)−Ny|y‾|2]
τk=kΔt
k=−(Nx−1)…Ny−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
rhoyx=zeros(Nx+Ny-1)+0j
tau=zeros(Nx+Ny-1)
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i]-mxe)*(y[i+k]-mye)
  
  rhoyx[k]=sum1/(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k)))*sqrt(vxe*vye))

plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
rhoyx=zeros(Nx+Ny-1)+0j
tau=roll(arange(-(Nx-1),Ny)*dt,Ny)
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
X=fft(append(x-mxe,zeros(Ny)))
Y=fft(append(y-mye,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
CEyx=ifft(Eyx)/float(dt)
for k in range(-(Nx-1),Ny):
  rhoyx[k]=CEyx[k]/float(dt)/(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k)))*sqrt(vxe*vye))

plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
show()
 direkte Berechnung:
rhoyx=zeros(1,Nx+Ny-1)+0j;
tau=zeros(1,Nx+Ny-1);
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+conj(x(i)-mxe)*(y(i+k)-mye);
end
rhoyx(mod(k,Nx+Ny-1)+1)=sum1/(min(min(Nx,Nx+k),min(Ny,Ny-k))*sqrt(vxe*vye));
end
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
rhoyx=zeros(1,Nx+Ny-1)+0j;
tau=circshift((-(Nx-1):Ny-1)*dt,[0;Ny]);
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
X=fft([x-mxe,zeros(1,Ny)]);
Y=fft([y-mye,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
CEyx=ifft(Eyx)/dt;
for k=-(Nx-1):Ny-1
rhoyx(mod(k,Nx+Ny-1)+1)=CEyx(mod(k,Nx+Ny)+1)/dt/(min(min(Nx,Nx+k),min(Ny,Ny-k))*sqrt(vxe*vye));
end
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Kreuzkorrelationskoeffizientenfunktion (mit lokaler Normierung, asymptotisch erwartungstreu) ρyx,k=ρyx(τk)=ρxy*(−τk)=∑i=max(0,−k)min(Nx,Ny−k)−1(xi−x‾)*(yi+k−y‾)[∑i=max(0,−k)min(Nx,Ny−k)−1|xi−x‾|2][∑i=max(0,k)min(Nx+k,Ny)−1|yi−y‾|2]
τk=kΔt
k=−(Nx−1)…Ny−1
(außerhalb dieses Bereichs nicht bekannt) 		direkte Berechnung:
rhoyx=zeros(Nx+Ny-1)+0j
tau=zeros(Nx+Ny-1)
mxe=mean(x)
mye=mean(y)
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  sum2=0+0j
  sum3=0+0j
  for i in range(maximum(0,-k),minimum(minimum(Nx,Ny-k)):
    sum1+=conj(x[i]-mxe)*(y[i+k]-mye)
    sum2+=abs(x[i]-mxe)**2
    sum3+=abs(y[i+k]-mye)**2
  
  rhoyx[k]=sum1/sqrt(sum2*sum3)

plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
rhoyx=zeros(Nx+Ny-1)+0j
tau=roll(arange(-(Nx-1),Ny)*dt,Ny)
mxe=mean(x)
mye=mean(y)
X=fft(append(x-mxe,zeros(Ny)))
Y=fft(append(y-mye,zeros(Nx)))
QX=fft(append(abs(x-mxe)**2,zeros(Ny)))
QY=fft(append(abs(y-mye)**2,zeros(Nx)))
OX=fft(append(ones(Nx),zeros(Ny)))
OY=fft(append(ones(Ny),zeros(Nx)))
Eyx=dt**2*conj(X)*Y
EA=dt**2*conj(QX)*OY
EB=dt**2*conj(OX)*QY
CEyx=ifft(Eyx)/float(dt)
CEA=ifft(EA)/float(dt)
CEB=ifft(EB)/float(dt)
for k in range(-(Nx-1),Ny):
  rhoyx[k]=CEyx[k]/sqrt(CEA[k]*CEB[k])

plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
show()
 direkte Berechnung:
rhoyx=zeros(1,Nx+Ny-1)+0j;
tau=zeros(1,Nx+Ny-1);
mxe=mean(x);
mye=mean(y);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
sum2=0+0j;
sum3=0+0j;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+conj(x(i)-mxe)*(y(i+k)-mye);
sum2=sum2+abs(x(i)-mxe)^2;
sum3=sum3+abs(y(i+k)-mye)^2;
end
rhoyx(mod(k,Nx+Ny-1)+1)=sum1/sqrt(sum2*sum3);
end
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
rhoyx=zeros(1,Nx+Ny-1)+0j;
tau=circshift((-(Nx-1):Ny-1)*dt,[0;Ny]);
mxe=mean(x);
mye=mean(y);
X=fft([x-mxe,zeros(1,Ny)]);
Y=fft([y-mye,zeros(1,Nx)]);
QX=fft([abs(x-mxe).^2,zeros(1,Ny)]);
QY=fft([abs(y-mye).^2,zeros(1,Nx)]);
OX=fft([ones(1,Nx),zeros(1,Ny)]);
OY=fft([ones(1,Ny),zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
EA=dt^2*conj(QX).*OY;
EB=dt^2*conj(OX).*QY;
CEyx=ifft(Eyx)/dt;
CEA=ifft(EA)/dt;
CEB=ifft(EB)/dt;
for k=-(Nx-1):Ny-1
rhoyx(mod(k,Nx+Ny-1)+1)=CEyx(mod(k,Nx+Ny)+1)/dt/sqrt(CEA(mod(k,Nx+Ny)+1)*CEB(mod(k,Nx+Ny)+1));
end
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Kreuzleistungsdichtespektrum (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≤2Nx−1,K≤2Ny−1)
Syx,j=Syx(fj)=Δt⋅FFT{Ryx,k}=Δt⋅∑k=−⌊K/2⌋⌊(K−1)/2⌋Ryx,k𝐞−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 des primären Spektrums:
from numpy.fft import *
X=fft(append(x,zeros(Ny)))
Y=fft(append(y,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
REyx=ifft(Eyx)/float(dt)
K=200
K=minimum(2*minimum(Nx,Ny)-1,K)
Ryx=zeros(K)+0j
for k in range(-(K//2),(K-1)//2+1):
  Ryx[k]=REyx[k]/float(dt)/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
plot(f,real(Syx),'o',f,imag(Syx),'o')
show()
Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
K=200
K=minimum(2*minimum(Nx,Ny)-1,K)
Ryx=zeros(K)+0j
for k in range(-(K//2),(K-1)//2+1):
  sum1=0+0j
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i])*y[i+k]
  
  Ryx[k]=sum1/float(minimum(Nx,minimum(Nx+k,minimum(Ny,Ny-k))))

f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
plot(f,real(Syx),'o',f,imag(Syx),'o')
show()
 Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung des primären Spektrums:
X=fft([x,zeros(1,Ny)]);
Y=fft([y,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
REyx=ifft(Eyx)/dt;
K=200;
K=min(2*min(Nx,Ny)-1,K);
Ryx=zeros(1,K)+0j;
for k=-floor(K/2):floor((K-1)/2)
Ryx(mod(k,K)+1)=REyx(mod(k,Nx+Ny)+1)/dt/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
plot(f,real(Syx),'o',f,imag(Syx),'o')
Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
K=200;
K=min(2*min(Nx,Ny)-1,K);
Ryx=zeros(1,K)+0j;
for k=-floor(K/2):floor((K-1)/2)
sum1=0+0j;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+conj(x(i))*y(i+k);
end
Ryx(mod(k,K)+1)=sum1/min(min(Nx,Nx+k),min(Ny,Ny-k));
end
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
plot(f,real(Syx),'o',f,imag(Syx),'o')