Signal- und Messdatenverarbeitung

Codes für komplexe, stochastische Leistungssignalpaare, mit 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 i=0…Nx−1 und Ny komplexen Werten yi mit i=0…Ny−1 (korrelierte und gegeneinander verschobene AR1-Prozesse als Beispiel, mit Erwartungswerten verschieden von null) mit einer von den Werten abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktionen als Beispiel)
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):
  good=False
  while (not good):
    good=True
    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 random()>1/(1+exp(-abs(x[i-maximum(0,Nd)]-mxs))):
        good=False
      
    
    if (i-maximum(0,-Nd)>=0) and (i-maximum(0,-Nd)<Ny):
      y[i-maximum(0,-Nd)]=c21*u[i]+c22*v[i]+mys
      if random()>1/(1+exp(-abs(y[i-maximum(0,-Nd)]-mys))):
        good=False
      
    
  

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
good=false;
while (~good)
good=true;
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;
if (rand>1/(1+exp(-abs(x(i-max(0,Nd))-mxs))))
good=false;
end
end
if (i-max(0,-Nd)>0) && (i-max(0,-Nd)<=Ny)
y(i-max(0,-Nd))=c21*u(i)+c22*v(i)+mys;
if (rand>1/(1+exp(-abs(y(i-max(0,-Nd))-mys))))
good=false;
end
end
end
end
plot(tx,real(x),'o',tx,imag(x),'o',ty,real(y),'o',ty,imag(y),'o')
Generieren von Nx Gewichten gi mit i=0…Nx−1, passend zu den xi mit i=0…Nx−1 und Ny Gewichten hi mit i=0…Ny−1, passend zu den yi mit i=0…Ny−1 (Inverse der verschobenen Sigmoidfunktion)
g=zeros(Nx)
for i in range(0,Nx):
  g[i]=(1+exp(-abs(x[i]-mxs)))**(1+phi)
  if (i+Nd>=0) and (i+Nd<Ny):
    g[i]*=(1+exp(-abs(y[i+Nd]-mys)))**(1+phi)
  

h=zeros(Ny)
for i in range(0,Ny):
  h[i]=(1+exp(-abs(y[i]-mys)))**(1+phi)
  if (i-Nd>=0) and (i-Nd<Nx):
    h[i]*=(1+exp(-abs(x[i-Nd]-mxs)))**(1+phi)
  

plot(tx,g,'o',ty,h,'o')
show()
g=zeros(1,Nx);
for i=1:Nx
g(i)=(1+exp(-abs(x(i)-mxs)))^(1+phi);
if ((i+Nd>0) && (i+Nd<=Ny))
g(i)=g(i)*(1+exp(-abs(y(i+Nd)-mys)))^(1+phi);
end
end
h=zeros(1,Ny);
for i=1:Ny
h(i)=(1+exp(-abs(y(i)-mys)))^(1+phi);
if ((i-Nd>0) && (i-Nd<=Nx))
h(i)=h(i)*(1+exp(-abs(x(i-Nd)-mxs)))^(1+phi);
end
end
plot(tx,g,'o',ty,h,'o')
Mittelwerte x‾=∑i=0Nx−1gixi∑i=0Nx−1gi
y‾=∑i=0Ny−1hiyi∑i=0Ny−1hi
mxe=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
Varianzen (ohne Bessel-Korrektur, asymptotisch erwartungstreu) sx2=∑i=0Nx−1gi(xi−x‾)*(xi−x‾)∑i=0Nx−1gi=∑i=0Nx−1gi|xi−x‾|2∑i=0Nx−1gi
=∑i=0Nx−1gixi*xi∑i=0Nx−1gi−x‾*x‾=∑i=0Nx−1gi|xi|2∑i=0Nx−1gi−|x‾|2
sy2=∑i=0Ny−1hi(yi−y‾)*(yi−y‾)∑i=0Ny−1hi=∑i=0Ny−1hi|yi−y‾|2∑i=0Ny−1hi
=∑i=0Ny−1hiyi*yi∑i=0Ny−1hi−y‾*y‾=∑i=0Ny−1hi|yi|2∑i=0Ny−1hi−|y‾|2
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(h*abs(y)**2)/float(sum(h))-abs(sum(h*y)/float(sum(h)))**2
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(h.*abs(y).^2)/sum(h)-abs(sum(h.*y)/sum(h))^2;
Varianzen (mit Bessel-Korrektur für korrelierte Datenreihen xi und yi, erwartungstreu) sx2=∑i=0Nx−1gi(xi−x‾)2∑i=0Nx−1gi+sx‾2=∑i=0Nx−1gixi2∑i=0Nx−1gi−x‾2+sx‾2
sy2=∑i=0Ny−1hi(yi−y‾)2∑i=0Ny−1hi+sy‾2=∑i=0Ny−1hiyi2∑i=0Ny−1hi−y‾2+sy‾2
mit den erwartungstreuen Schätzungen sx‾2 und sy‾2 der Schätzvarianzen der Mittelwertschätzer für korrelierte Daten aus den erwartungstreuen Schätzungen der (Auto-)Kovarianzfunktionen 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
mxe=sum(g*x)/float(sum(g))
xe=append(g*(x-mxe),zeros(Nx))
ge=append(g,zeros(Nx))
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+=(Nx-abs(k))*Ce[k]

vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2+real(sC)/float(Nx**2)
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
mye=sum(h*y)/float(sum(h))
ye=append(h*(y-mye),zeros(Ny))
he=append(h,zeros(Ny))
X=ifft(abs(fft(ye))**2)
Y=real(ifft(abs(fft(he))**2))
Cp=zeros(2*Nc-1)+0j
for k in range(-(Nc-1),Nc):
  Cp[k]=X[k]/Y[k]

D=sum(h)
A=zeros([2*Nc-1,2*Nc-1])
for k in range(-(Nc-1),Nc):
  Gk=real(ifft(conj(fft(he*roll(he,-k)))*fft(he)))
  Hk=real(ifft(conj(fft(he))*fft(he*roll(he,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+=(Ny-abs(k))*Ce[k]

vye=sum(h*abs(y)**2)/float(sum(h))-abs(sum(h*y)/float(sum(h)))**2+real(sC)/float(Ny**2)
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
mxe=sum(g.*x)/sum(g);
xe=[g.*(x-mxe),zeros(1,Nx)];
ge=[g,zeros(1,Nx)];
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*Nx)+1)/Y(mod(k,2*Nx)+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*Nx)+1)+Hk(mod(j,2*Nx)+1))/(Y(mod(k,2*Nx)+1)*D)+Y(mod(j,2*Nx)+1)/(D^2);
end
end
Ainv=inv(A);
Ce=Ainv*Cp';
sC=0+0j;
for k=-(Nc-1):Nc-1
sC=sC+(Nx-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2+real(sC)/(Nx^2);
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
mye=sum(h.*y)/sum(h);
ye=[h.*(y-mye),zeros(1,Ny)];
he=[h,zeros(1,Ny)];
X=ifft(abs(fft(ye)).^2);
Y=real(ifft(abs(fft(he)).^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*Ny)+1)/Y(mod(k,2*Ny)+1);
end
D=sum(h);
A=zeros(2*Nc-1,2*Nc-1);
for k=-(Nc-1):Nc-1
Gk=real(ifft(conj(fft(he.*circshift(he,[0;-k]))).*fft(he)));
Hk=real(ifft(conj(fft(he)).*fft(he.*circshift(he,[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*Ny)+1)+Hk(mod(j,2*Ny)+1))/(Y(mod(k,2*Ny)+1)*D)+Y(mod(j,2*Ny)+1)/(D^2);
end
end
Ainv=inv(A);
Ce=Ainv*Cp';
sC=0+0j;
for k=-(Nc-1):Nc-1
sC=sC+(Ny-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vye=sum(h.*abs(y).^2)/sum(h)-abs(sum(h.*y)/sum(h))^2+real(sC)/(Ny^2);
Kreuzkovarianz (ohne Bessel-Korrektur, asymptotisch erwartungstreu) cyx=cxy*=∑i=0min(Nx,Ny)−1gihi(xi−x‾)*(yi−y‾)∑i=0min(Nx,Ny)−1gihi
sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*conj(x[0:minimum(Nx,Ny)]-sum(g*x)/float(sum(g)))*(y[0:minimum(Nx,Ny)]-sum(h*y)/float(sum(h))))/float(sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]))
sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*conj(x(1:min(Nx,Ny))-sum(g.*x)/sum(g)).*(y(1:min(Nx,Ny))-sum(h.*y)/sum(h)))/sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)))
Kreuzkovarianz (mit Bessel-Korrektur für korrelierte Datenreihen xi und yi, erwartungstreu) cyx=γyxsx2sy2
mit den erwartungstreuen Schätzungen der Varianzen sx2 und sy2 für korrelierte Daten zuvor und der folgenden asymptotisch erwartungstreuen Schätzung des (Kreuz-)Korrelationskoeffizienten γyx
from numpy.fft import *
from numpy.linalg import *
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
mxe=sum(g*x)/float(sum(g))
xe=append(g*(x-mxe),zeros(Nx))
ge=append(g,zeros(Nx))
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+=(Nx-abs(k))*Ce[k]

vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2+real(sC)/float(Nx**2)
Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
mye=sum(h*y)/float(sum(h))
ye=append(h*(y-mye),zeros(Ny))
he=append(h,zeros(Ny))
X=ifft(abs(fft(ye))**2)
Y=real(ifft(abs(fft(he))**2))
Cp=zeros(2*Nc-1)+0j
for k in range(-(Nc-1),Nc):
  Cp[k]=X[k]/Y[k]

D=sum(h)
A=zeros([2*Nc-1,2*Nc-1])
for k in range(-(Nc-1),Nc):
  Gk=real(ifft(conj(fft(he*roll(he,-k)))*fft(he)))
  Hk=real(ifft(conj(fft(he))*fft(he*roll(he,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+=(Ny-abs(k))*Ce[k]

vye=sum(h*abs(y)**2)/float(sum(h))-abs(sum(h*y)/float(sum(h)))**2+real(sC)/float(Ny**2)
sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*conj(x[0:minimum(Nx,Ny)]-mxe)*(y[0:minimum(Nx,Ny)]-mye))*sqrt(sum(g)*sum(h))/(sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)])*sqrt(sum(g*abs(x-mxe)**2)*sum(h*abs(y-mye)**2)))*sqrt(vxe*vye)
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Nx
mxe=sum(g.*x)/sum(g);
xe=[g.*(x-mxe),zeros(1,Nx)];
ge=[g,zeros(1,Nx)];
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*Nx)+1)/Y(mod(k,2*Nx)+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*Nx)+1)+Hk(mod(j,2*Nx)+1))/(Y(mod(k,2*Nx)+1)*D)+Y(mod(j,2*Nx)+1)/(D^2);
end
end
Ainv=inv(A);
Ce=Ainv*Cp';
sC=0+0j;
for k=-(Nc-1):Nc-1
sC=sC+(Nx-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2+real(sC)/(Nx^2);
Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als Ny
mye=sum(h.*y)/sum(h);
ye=[h.*(y-mye),zeros(1,Ny)];
he=[h,zeros(1,Ny)];
X=ifft(abs(fft(ye)).^2);
Y=real(ifft(abs(fft(he)).^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*Ny)+1)/Y(mod(k,2*Ny)+1);
end
D=sum(h);
A=zeros(2*Nc-1,2*Nc-1);
for k=-(Nc-1):Nc-1
Gk=real(ifft(conj(fft(he.*circshift(he,[0;-k]))).*fft(he)));
Hk=real(ifft(conj(fft(he)).*fft(he.*circshift(he,[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*Ny)+1)+Hk(mod(j,2*Ny)+1))/(Y(mod(k,2*Ny)+1)*D)+Y(mod(j,2*Ny)+1)/(D^2);
end
end
Ainv=inv(A);
Ce=Ainv*Cp';
sC=0+0j;
for k=-(Nc-1):Nc-1
sC=sC+(Ny-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vye=sum(h.*abs(y).^2)/sum(h)-abs(sum(h.*y)/sum(h))^2+real(sC)/(Ny^2);
sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*conj(x(1:min(Nx,Ny))-mxe).*(y(1:min(Nx,Ny))-mye))*sqrt(sum(g)*sum(h))/(sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)))*sqrt(sum(g.*abs(x-mxe).^2)*sum(h.*abs(y-mye).^2)))*sqrt(vxe*vye)
(Kreuz-)Korrelationskoeffizient (asymptotisch erwartungstreu) γyx=γxy*=[∑i=0min(Nx,Ny)−1gihi(xi−x‾)*(yi−y‾)]⋅(∑i=0Nx−1gi)(∑i=0Ny−1hi)[∑i=0min(Nx,Ny)−1gihi]⋅(∑i=0Nx−1gi|xi−x‾|2)(∑i=0Ny−1hi|yi−y‾|2)
mxe=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*conj(x[0:minimum(Nx,Ny)]-mxe)*(y[0:minimum(Nx,Ny)]-mye))*sqrt(sum(g)*sum(h))/(sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)])*sqrt(sum(g*abs(x-mxe)**2)*sum(h*abs(y-mye)**2)))
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*conj(x(1:min(Nx,Ny))-mxe).*(y(1:min(Nx,Ny))-mye))*sqrt(sum(g)*sum(h))/(sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)))*sqrt(sum(g.*abs(x-mxe).^2)*sum(h.*abs(y-mye).^2)))
(Kreuz-)Korrelationskoeffizient (mit lokaler Normierung, asymptotisch erwartungstreu) γyx=γxy*=∑i=0min(Nx,Ny)−1gihi(xi−x‾)*(yi−y‾)[∑i=0min(Nx,Ny)−1gihi|xi−x‾|2][∑i=0min(Nx,Ny)−1gihi|yi−y‾|2]
mxe=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*conj(x[0:minimum(Nx,Ny)]-mxe)*(y[0:minimum(Nx,Ny)]-mye))/sqrt(sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*abs(x[0:minimum(Nx,Ny)]-mxe)**2)*sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*abs(y[0:minimum(Nx,Ny)]-mye)**2))
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*conj(x(1:min(Nx,Ny))-mxe).*(y(1:min(Nx,Ny))-mye))/sqrt(sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*abs(x(1:min(Nx,Ny))-mxe).^2)*sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*abs(y(1:min(Nx,Ny))-mye).^2))
Kreuzkorrelationsfunktion (erwartungstreu) Ryx,k=Ryx(τk)=Rxy*(−τk)=∑i=max(0,−k)min(Nx,Ny−k)−1gihi+kxi*yi+k∑i=max(0,−k)min(Nx,Ny−k)−1gihi+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
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*conj(x[i])*y[i+k]
    sum0+=g[i]*h[i+k]
  
  Ryx[k]=sum1/float(sum0)

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)
X1=fft(append(g*x,zeros(Ny)))
Y1=fft(append(h*y,zeros(Nx)))
Eyx1=dt**2*conj(X1)*Y1
REyx1=ifft(Eyx1)/float(dt)
X0=fft(append(g,zeros(Ny)))
Y0=fft(append(h,zeros(Nx)))
Eyx0=dt**2*conj(X0)*Y0
REyx0=real(ifft(Eyx0))/float(dt)
for k in range(-(Nx-1),Ny):
  Ryx[k]=REyx1[k]/REyx0[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;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*conj(x(i))*y(i+k);
sum0=sum0+g(i)*h(i+k);
end
Ryx(mod(k,Nx+Ny-1)+1)=sum1/sum0;
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]);
X1=fft([g.*x,zeros(1,Ny)]);
Y1=fft([h.*y,zeros(1,Nx)]);
Eyx1=dt^2*conj(X1).*Y1;
REyx1=ifft(Eyx1)/dt;
X0=fft([g,zeros(1,Ny)]);
Y0=fft([h,zeros(1,Nx)]);
Eyx0=dt^2*conj(X0).*Y0;
REyx0=real(ifft(Eyx0))/dt;
for k=-(Nx-1):Ny-1
Ryx(mod(k,Nx+Ny-1)+1)=REyx1(mod(k,Nx+Ny)+1)/REyx0(mod(k,Nx+Ny)+1);
end
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Kreuzkovarianzfunktion (ohne Bessel-Korrektur, asymptotisch erwartungstreu) Cyx,k=Cyx(τk)=Cxy*(−τk)=∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k(xi−x‾)*(yi+k−y‾)∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k
τ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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*conj(x[i]-mxe)*(y[i+k]-mye)
    sum0+=g[i]*h[i+k]
  
  Cyx[k]=sum1/float(sum0)

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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
X1=fft(append(g*(x-mxe),zeros(Ny)))
Y1=fft(append(h*(y-mye),zeros(Nx)))
Eyx1=dt**2*conj(X1)*Y1
REyx1=ifft(Eyx1)/float(dt)
X0=fft(append(g,zeros(Ny)))
Y0=fft(append(h,zeros(Nx)))
Eyx0=dt**2*conj(X0)*Y0
REyx0=real(ifft(Eyx0))/float(dt)
for k in range(-(Nx-1),Ny):
  Cyx[k]=REyx1[k]/real(REyx0[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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*conj(x(i)-mxe)*(y(i+k)-mye);
sum0=sum0+g(i)*h(i+k);
end
Cyx(mod(k,Nx+Ny-1)+1)=sum1/sum0;
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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
X1=fft([g.*(x-mxe),zeros(1,Ny)]);
Y1=fft([h.*(y-mye),zeros(1,Nx)]);
Eyx1=dt^2*conj(X1).*Y1;
REyx1=ifft(Eyx1)/dt;
X0=fft([g,zeros(1,Ny)]);
Y0=fft([h,zeros(1,Nx)]);
Eyx0=dt^2*conj(X0).*Y0;
REyx0=real(ifft(Eyx0))/dt;
for k=-(Nx-1):Ny-1
Cyx(mod(k,Nx+Ny-1)+1)=REyx1(mod(k,Nx+Ny)+1)/real(REyx0(mod(k,Nx+Ny)+1));
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
S=zeros(Nxc+Nyc-1)
tau=zeros(Nxc+Nyc-1)
mxe=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
for k in range(-(Nxc-1),Nyc):
  tau[k]=k*dt
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*conj(x[i]-mxe)*(y[i+k]-mye)
    sum0+=g[i]*h[i+k]
  
  Cyxp[k]=sum1/float(sum0)
  S[k]=sum0

Dx=sum(g)
Dy=sum(h)
A=zeros([Nxc+Nyc-1,Nxc+Nyc-1])
for k in range(-(Nxc-1),Nyc):
  for j in range(-(Nxc-1),Nyc):
    s1=0
    for i in range(maximum(0,maximum(-j,-k)),minimum(Nx,minimum(Ny-j,Ny-k))):
      s1+=g[i]*h[i+j]*h[i+k]
    s2=0
    for i in range(maximum(0,maximum(-j,k-j)),minimum(Nx,minimum(Ny-j,Nx+k-j))):
      s2+=g[i]*h[i+j]*g[i+j-k]
    A[k,j]=int(k==j)-s1/float(S[k]*Dy)-s2/float(S[k]*Dx)+S[j]/float(Dx*Dy)
  

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
tau=roll(arange(-(Nxc-1),Nyc)*dt,Nyc)
mxe=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
xe=append(g*(x-mxe),zeros(Ny))
ge=append(g,zeros(Ny))
ye=append(h*(y-mye),zeros(Nx))
he=append(h,zeros(Nx))
X=ifft(conj(fft(xe))*fft(ye))
Y=real(ifft(conj(fft(ge))*fft(he)))
Cyxp=zeros(Nxc+Nyc-1)+0j
for k in range(-(Nxc-1),Nyc):
  Cyxp[k]=X[k]/Y[k]

Dx=sum(g)
Dy=sum(h)
A=zeros([Nxc+Nyc-1,Nxc+Nyc-1])
for k in range(-(Nxc-1),Nyc):
  Gk=real(ifft(conj(fft(ge*roll(he,-k)))*fft(he)))
  Hk=real(ifft(conj(fft(ge))*fft(he*roll(ge,k))))
  for j in range(-(Nxc-1),Nyc):
    A[k,j]=int(k==j)+Y[j]/float(Dx*Dy)-Gk[j]/float(Y[k]*Dy)-Hk[j]/float(Y[k]*Dx)
  

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;
S=zeros(1,Nxc+Nyc-1);
tau=zeros(1,Nxc+Nyc-1);
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
for k=-(Nxc-1):Nyc-1
tau(mod(k,Nxc+Nyc-1)+1)=k*dt;
sum1=0+0j;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*conj(x(i)-mxe)*(y(i+k)-mye);
sum0=sum0+g(i)*h(i+k);
end
Cyxp(mod(k,Nxc+Nyc-1)+1)=sum1/sum0;
S(mod(k,Nxc+Nyc-1)+1)=sum0;
end
Dx=sum(g);
Dy=sum(h);
A=zeros(Nxc+Nyc-1,Nxc+Nyc-1);
for k=-(Nxc-1):Nyc-1
for j=-(Nxc-1):Nyc-1
s1=0;
for i=max(0,max(-j,-k))+1:min(Nx,min(Ny-j,Ny-k))
s1=s1+g(i)*h(i+j)*h(i+k);
end
s2=0;
for i=max(0,max(-j,k-j))+1:min(Nx,min(Ny-j,Nx+k-j))
s2=s2+g(i)*h(i+j)*g(i+j-k);
end
A(mod(k,Nxc+Nyc-1)+1,mod(j,Nxc+Nyc-1)+1)=(k==j)-s1/(S(mod(k,Nxc+Nyc-1)+1)*Dy)-s2/(S(mod(k,Nxc+Nyc-1)+1)*Dx)+S(mod(j,Nxc+Nyc-1)+1)/(Dx*Dy);
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
tau=circshift((-(Nxc-1):Nyc-1)*dt,[0;Nyc]);
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
xe=[g.*(x-mxe),zeros(1,Ny)];
ge=[g,zeros(1,Ny)];
ye=[h.*(y-mye),zeros(1,Nx)];
he=[h,zeros(1,Nx)];
X=ifft(conj(fft(xe)).*fft(ye));
Y=real(ifft(conj(fft(ge)).*fft(he)));
Cyxp=zeros(1,Nxc+Nyc-1)+0j;
for k=-(Nxc-1):Nyc-1
Cyxp(mod(k,Nxc+Nyc-1)+1)=X(mod(k,Nx+Ny)+1)/Y(mod(k,Nx+Ny)+1);
end
Dx=sum(g);
Dy=sum(h);
A=zeros(Nxc+Nyc-1,Nxc+Nyc-1);
for k=-(Nxc-1):Nyc-1
Gk=real(ifft(conj(fft(ge.*circshift(he,[0;-k]))).*fft(he)));
Hk=real(ifft(conj(fft(ge)).*fft(he.*circshift(ge,[0;k]))));
for j=-(Nxc-1):Nyc-1
A(mod(k,Nxc+Nyc-1)+1,mod(j,Nxc+Nyc-1)+1)=(k==j)+Y(mod(j,Nx+Ny)+1)/(Dx*Dy)-Gk(mod(j,Nx+Ny)+1)/(Y(mod(k,Nx+Ny)+1)*Dy)-Hk(mod(j,Nx+Ny)+1)/(Y(mod(k,Nx+Ny)+1)*Dx);
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)−1gihi+k(xi−x‾)*(yi+k−y‾)]⋅sx2sy2[∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k|xi−x‾|2][∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k|yi+k−y‾|2]
=[∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k(xi−x‾)*(yi+k−y‾)]⋅(∑i=0Nx−1gi|xi−x‾|2)(∑i=0Ny−1hi|yi−y‾|2)(∑i=0Nx−1gi)(∑i=0Ny−1hi)[∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k|xi−x‾|2][∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k|yi+k−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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(h*abs(y)**2)/float(sum(h))-abs(sum(h*y)/float(sum(h)))**2
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  sum2=0
  sum3=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*conj(x[i]-mxe)*(y[i+k]-mye)
    sum2+=g[i]*h[i+k]*abs(x[i]-mxe)**2
    sum3+=g[i]*h[i+k]*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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(h*abs(y)**2)/float(sum(h))-abs(sum(h*y)/float(sum(h)))**2
X=fft(append(g*(x-mxe),zeros(Ny)))
Y=fft(append(h*(y-mye),zeros(Nx)))
QX=fft(append(g*abs(x-mxe)**2,zeros(Ny)))
QY=fft(append(h*abs(y-mye)**2,zeros(Nx)))
OX=fft(append(g,zeros(Ny)))
OY=fft(append(h,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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(h.*abs(y).^2)/sum(h)-abs(sum(h.*y)/sum(h))^2;
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
sum2=0;
sum3=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*conj(x(i)-mxe)*(y(i+k)-mye);
sum2=sum2+g(i)*h(i+k)*abs(x(i)-mxe)^2;
sum3=sum3+g(i)*h(i+k)*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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(h.*abs(y).^2)/sum(h)-abs(sum(h.*y)/sum(h))^2;
X=fft([g.*(x-mxe),zeros(1,Ny)]);
Y=fft([h.*(y-mye),zeros(1,Nx)]);
QX=fft([g.*abs(x-mxe).^2,zeros(1,Ny)]);
QY=fft([h.*abs(y-mye).^2,zeros(1,Nx)]);
OX=fft([g,zeros(1,Ny)]);
OY=fft([h,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)/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)=∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k(xi−x‾)*(yi+k−y‾)[∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k]⋅sx2sy2
=[∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k(xi−x‾)*(yi+k−y‾)]⋅(∑i=0Nx−1gi)(∑i=0Ny−1hi)[∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k]⋅(∑i=0Nx−1gi|xi−x‾|2)(∑i=0Ny−1hi|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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(h*abs(y)**2)/float(sum(h))-abs(sum(h*y)/float(sum(h)))**2
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*conj(x[i]-mxe)*(y[i+k]-mye)
    sum0+=g[i]*h[i+k]
  
  rhoyx[k]=sum1/float(sum0)/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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(h*abs(y)**2)/float(sum(h))-abs(sum(h*y)/float(sum(h)))**2
X1=fft(append(g*(x-mxe),zeros(Ny)))
Y1=fft(append(h*(y-mye),zeros(Nx)))
Eyx1=dt**2*conj(X1)*Y1
REyx1=ifft(Eyx1)/float(dt)
X0=fft(append(g,zeros(Ny)))
Y0=fft(append(h,zeros(Nx)))
Eyx0=dt**2*conj(X0)*Y0
REyx0=real(ifft(Eyx0))/float(dt)
for k in range(-(Nx-1),Ny):
  rhoyx[k]=REyx1[k]/real(REyx0[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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(h.*abs(y).^2)/sum(h)-abs(sum(h.*y)/sum(h))^2;
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*conj(x(i)-mxe)*(y(i+k)-mye);
sum0=sum0+g(i)*h(i+k);
end
rhoyx(mod(k,Nx+Ny-1)+1)=sum1/sum0/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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(h.*abs(y).^2)/sum(h)-abs(sum(h.*y)/sum(h))^2;
X1=fft([g.*(x-mxe),zeros(1,Ny)]);
Y1=fft([h.*(y-mye),zeros(1,Nx)]);
Eyx1=dt^2*conj(X1).*Y1;
REyx1=ifft(Eyx1)/dt;
X0=fft([g,zeros(1,Ny)]);
Y0=fft([h,zeros(1,Nx)]);
Eyx0=dt^2*conj(X0).*Y0;
REyx0=real(ifft(Eyx0))/dt;
for k=-(Nx-1):Ny-1
rhoyx(mod(k,Nx+Ny-1)+1)=REyx1(mod(k,Nx+Ny)+1)/real(REyx0(mod(k,Nx+Ny)+1))/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)−1gihi+k(xi−x‾)*(yi+k−y‾)[∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k|xi−x‾|2][∑i=max(0,−k)min(Nx,Ny−k)−1gihi+k|yi+k−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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0+0j
  sum2=0
  sum3=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*conj(x[i]-mxe)*(y[i+k]-mye)
    sum2+=g[i]*h[i+k]*abs(x[i]-mxe)**2
    sum3+=g[i]*h[i+k]*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=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
X=fft(append(g*(x-mxe),zeros(Ny)))
Y=fft(append(h*(y-mye),zeros(Nx)))
QX=fft(append(g*abs(x-mxe)**2,zeros(Ny)))
QY=fft(append(h*abs(y-mye)**2,zeros(Nx)))
OX=fft(append(g,zeros(Ny)))
OY=fft(append(h,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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0+0j;
sum2=0;
sum3=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*conj(x(i)-mxe)*(y(i+k)-mye);
sum2=sum2+g(i)*h(i+k)*abs(x(i)-mxe)^2;
sum3=sum3+g(i)*h(i+k)*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=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
X=fft([g.*(x-mxe),zeros(1,Ny)]);
Y=fft([h.*(y-mye),zeros(1,Nx)]);
QX=fft([g.*abs(x-mxe).^2,zeros(1,Ny)]);
QY=fft([h.*abs(y-mye).^2,zeros(1,Nx)]);
OX=fft([g,zeros(1,Ny)]);
OY=fft([h,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)/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 der primären Spektren:
from numpy.fft import *
X1=fft(append(g*x,zeros(Ny)))
Y1=fft(append(h*y,zeros(Nx)))
E1=dt**2*conj(X1)*Y1
X0=fft(append(g,zeros(Ny)))
Y0=fft(append(h,zeros(Nx)))
E0=dt**2*conj(X0)*Y0
RA=ifft(E1)/real(ifft(E0))
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]=RA[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
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*conj(x[i])*y[i+k]
    sum0+=g[i]*h[i+k]
  
  Ryx[k]=sum1/float(sum0)

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 der primären Spektren:
X1=fft([g.*x,zeros(1,Ny)]);
Y1=fft([h.*y,zeros(1,Nx)]);
E1=dt^2*conj(X1).*Y1;
X0=fft([g,zeros(1,Ny)]);
Y0=fft([h,zeros(1,Nx)]);
E0=dt^2*conj(X0).*Y0;
RA=ifft(E1)./real(ifft(E0));
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)=RA(mod(k,Nx+Ny)+1);
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;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*conj(x(i))*y(i+k);
sum0=sum0+g(i)*h(i+k);
end
Ryx(mod(k,K)+1)=sum1/sum0;
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')