Signal- und Messdatenverarbeitung

Codes für reelle, 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 Werten xi mit i=0…Nx−1 und Ny Werten yi mit i=0…Ny−1 (zwei 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
mys=2
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)
v=zeros(Nuv)
x=zeros(Nx)
y=zeros(Ny)
for i in range(0,Nuv):
  good=False
  while (not good):
    good=True
    if (i==0):
      u[i]=standard_normal()
      v[i]=standard_normal()
    
    else:
      u[i]=phi*u[i-1]+normal(0,theta)
      v[i]=phi*v[i-1]+normal(0,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(-(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(-(y[i-maximum(0,-Nd)]-mys))):
        good=False
      
    
  

plot(tx,x,'o',ty,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;
mys=2;
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);
v=zeros(1,Nuv);
x=zeros(1,Nx);
y=zeros(1,Ny);
for i=1:Nuv
good=false;
while (~good)
good=true;
if (i==1)
u(i)=randn();
v(i)=randn();
else
u(i)=phi*u(i-1)+theta*randn();
v(i)=phi*v(i-1)+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(-(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(-(y(i-max(0,-Nd))-mys))))
good=false;
end
end
end
end
plot(tx,x,'o',ty,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(-(x[i]-mxs)))**(1+phi)
  if (i+Nd>=0) and (i+Nd<Ny):
    g[i]*=(1+exp(-(y[i+Nd]-mys)))**(1+phi)
  

h=zeros(Ny)
for i in range(0,Ny):
  h[i]=(1+exp(-(y[i]-mys)))**(1+phi)
  if (i-Nd>=0) and (i-Nd<Nx):
    h[i]*=(1+exp(-(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(-(x(i)-mxs)))^(1+phi);
if ((i+Nd>0) && (i+Nd<=Ny))
g(i)=g(i)*(1+exp(-(y(i+Nd)-mys)))^(1+phi);
end
end
h=zeros(1,Ny);
for i=1:Ny
h(i)=(1+exp(-(y(i)-mys)))^(1+phi);
if ((i-Nd>0) && (i-Nd<=Nx))
h(i)=h(i)*(1+exp(-(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‾)2∑i=0Nx−1gi=∑i=0Nx−1gixi2∑i=0Nx−1gi−x‾2
sy2=∑i=0Ny−1hi(yi−y‾)2∑i=0Ny−1hi=∑i=0Ny−1hiyi2∑i=0Ny−1hi−y‾2
vxe=sum(g*x**2)/float(sum(g))-(sum(g*x)/float(sum(g)))**2
vye=sum(h*y**2)/float(sum(h))-(sum(h*y)/float(sum(h)))**2
vxe=sum(g.*x.^2)/sum(g)-(sum(g.*x)/sum(g))^2;
vye=sum(h.*y.^2)/sum(h)-(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=real(ifft(abs(fft(xe))**2))
Y=real(ifft(abs(fft(ge))**2))
Cp=zeros(2*Nc-1)
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
for k in range(-(Nc-1),Nc):
  sC+=(Nx-abs(k))*Ce[k]

vxe=sum(g*x**2)/float(sum(g))-(sum(g*x)/float(sum(g)))**2+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=real(ifft(abs(fft(ye))**2))
Y=real(ifft(abs(fft(he))**2))
Cp=zeros(2*Nc-1)
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
for k in range(-(Nc-1),Nc):
  sC+=(Ny-abs(k))*Ce[k]

vye=sum(h*y**2)/float(sum(h))-(sum(h*y)/float(sum(h)))**2+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=real(ifft(abs(fft(xe)).^2));
Y=real(ifft(abs(fft(ge)).^2));
Cp=zeros(1,2*Nc-1);
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;
for k=-(Nc-1):Nc-1
sC=sC+(Nx-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vxe=sum(g.*x.^2)/sum(g)-(sum(g.*x)/sum(g))^2+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=real(ifft(abs(fft(ye)).^2));
Y=real(ifft(abs(fft(he)).^2));
Cp=zeros(1,2*Nc-1);
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;
for k=-(Nc-1):Nc-1
sC=sC+(Ny-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vye=sum(h.*y.^2)/sum(h)-(sum(h.*y)/sum(h))^2+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)]*(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)).*(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=real(ifft(abs(fft(xe))**2))
Y=real(ifft(abs(fft(ge))**2))
Cp=zeros(2*Nc-1)
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
for k in range(-(Nc-1),Nc):
  sC+=(Nx-abs(k))*Ce[k]

vxe=sum(g*x**2)/float(sum(g))-(sum(g*x)/float(sum(g)))**2+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=real(ifft(abs(fft(ye))**2))
Y=real(ifft(abs(fft(he))**2))
Cp=zeros(2*Nc-1)
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
for k in range(-(Nc-1),Nc):
  sC+=(Ny-abs(k))*Ce[k]

vye=sum(h*y**2)/float(sum(h))-(sum(h*y)/float(sum(h)))**2+sC/float(Ny**2)
sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*(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*(x-mxe)**2)*sum(h*(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=real(ifft(abs(fft(xe)).^2));
Y=real(ifft(abs(fft(ge)).^2));
Cp=zeros(1,2*Nc-1);
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;
for k=-(Nc-1):Nc-1
sC=sC+(Nx-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vxe=sum(g.*x.^2)/sum(g)-(sum(g.*x)/sum(g))^2+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=real(ifft(abs(fft(ye)).^2));
Y=real(ifft(abs(fft(he)).^2));
Cp=zeros(1,2*Nc-1);
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;
for k=-(Nc-1):Nc-1
sC=sC+(Ny-abs(k))*Ce(mod(k,2*Nc-1)+1);
end
vye=sum(h.*y.^2)/sum(h)-(sum(h.*y)/sum(h))^2+sC/(Ny^2);
sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*(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.*(x-mxe).^2)*sum(h.*(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)]*(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*(x-mxe)**2)*sum(h*(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)).*(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.*(x-mxe).^2)*sum(h.*(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)]*(x[0:minimum(Nx,Ny)]-mxe)*(y[0:minimum(Nx,Ny)]-mye))/sqrt(sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*(x[0:minimum(Nx,Ny)]-mxe)**2)*sum(g[0:minimum(Nx,Ny)]*h[0:minimum(Nx,Ny)]*(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)).*(x(1:min(Nx,Ny))-mxe).*(y(1:min(Nx,Ny))-mye))/sqrt(sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*(x(1:min(Nx,Ny))-mxe).^2)*sum(g(1:min(Nx,Ny)).*h(1:min(Nx,Ny)).*(y(1:min(Nx,Ny))-mye).^2))
Kreuzkorrelationsfunktion (erwartungstreu) Ryx,k=Ryx(τk)=Rxy(−τk)=∑i=max(0,−k)min(Nx,Ny−k)−1gihi+kxiyi+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)
tau=zeros(Nx+Ny-1)
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*x[i]*y[i+k]
    sum0+=g[i]*h[i+k]
  
  Ryx[k]=sum1/float(sum0)

plot(tau,Ryx,'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
Ryx=zeros(Nx+Ny-1)
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=real(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,Ryx,'o')
show()
 direkte Berechnung:
Ryx=zeros(1,Nx+Ny-1);
tau=zeros(1,Nx+Ny-1);
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*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,Ryx,'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Ryx=zeros(1,Nx+Ny-1);
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=real(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,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)
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
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*(x[i]-mxe)*(y[i+k]-mye)
    sum0+=g[i]*h[i+k]
  
  Cyx[k]=sum1/float(sum0)

plot(tau,Cyx,'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
Cyx=zeros(Nx+Ny-1)
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
CEyx1=real(ifft(Eyx1))/float(dt)
X0=fft(append(g,zeros(Ny)))
Y0=fft(append(h,zeros(Nx)))
Eyx0=dt**2*conj(X0)*Y0
CEyx0=real(ifft(Eyx0))/float(dt)
for k in range(-(Nx-1),Ny):
  Cyx[k]=CEyx1[k]/CEyx0[k]

plot(tau,Cyx,'o')
show()
 direkte Berechnung:
Cyx=zeros(1,Nx+Ny-1);
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;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*(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,Cyx,'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Cyx=zeros(1,Nx+Ny-1);
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;
CEyx1=real(ifft(Eyx1))/dt;
X0=fft([g,zeros(1,Ny)]);
Y0=fft([h,zeros(1,Nx)]);
Eyx0=dt^2*conj(X0).*Y0;
CEyx0=real(ifft(Eyx0))/dt;
for k=-(Nx-1):Ny-1
Cyx(mod(k,Nx+Ny-1)+1)=CEyx1(mod(k,Nx+Ny)+1)/CEyx0(mod(k,Nx+Ny)+1);
end
plot(tau,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)
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
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*(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,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=real(ifft(conj(fft(xe))*fft(ye)))
Y=real(ifft(conj(fft(ge))*fft(he)))
Cyxp=zeros(Nxc+Nyc-1)
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,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);
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;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*(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,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=real(ifft(conj(fft(xe)).*fft(ye)));
Y=real(ifft(conj(fft(ge)).*fft(he)));
Cyxp=zeros(1,Nxc+Nyc-1);
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,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)
tau=zeros(Nx+Ny-1)
mxe=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
vxe=sum(g*x**2)/float(sum(g))-(sum(g*x)/float(sum(g)))**2
vye=sum(h*y**2)/float(sum(h))-(sum(h*y)/float(sum(h)))**2
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0
  sum2=0
  sum3=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*(x[i]-mxe)*(y[i+k]-mye)
    sum2+=g[i]*h[i+k]*(x[i]-mxe)**2
    sum3+=g[i]*h[i+k]*(y[i+k]-mye)**2
  
  Cyx[k]=sqrt(vxe*vye)*sum1/sqrt(sum2*sum3)

plot(tau,Cyx,'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
Cyx=zeros(Nx+Ny-1)
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*x**2)/float(sum(g))-(sum(g*x)/float(sum(g)))**2
vye=sum(h*y**2)/float(sum(h))-(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*(x-mxe)**2,zeros(Ny)))
QY=fft(append(h*(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=real(ifft(Eyx))/float(dt)
CEA=real(ifft(EA))/float(dt)
CEB=real(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,Cyx,'o')
show()
 direkte Berechnung:
Cyx=zeros(1,Nx+Ny-1);
tau=zeros(1,Nx+Ny-1);
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*x.^2)/sum(g)-(sum(g.*x)/sum(g))^2;
vye=sum(h.*y.^2)/sum(h)-(sum(h.*y)/sum(h))^2;
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0;
sum2=0;
sum3=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*(x(i)-mxe)*(y(i+k)-mye);
sum2=sum2+g(i)*h(i+k)*(x(i)-mxe)^2;
sum3=sum3+g(i)*h(i+k)*(y(i+k)-mye)^2;
end
Cyx(mod(k,Nx+Ny-1)+1)=sqrt(vxe*vye)*sum1/sqrt(sum2*sum3);
end
plot(tau,Cyx,'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
Cyx=zeros(1,Nx+Ny-1);
tau=circshift((-(Nx-1):Ny-1)*dt,[0;Ny]);
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*x.^2)/sum(g)-(sum(g.*x)/sum(g))^2;
vye=sum(h.*y.^2)/sum(h)-(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.*(x-mxe).^2,zeros(1,Ny)]);
QY=fft([h.*(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=real(ifft(Eyx))/dt;
CEA=real(ifft(EA))/dt;
CEB=real(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,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)
tau=zeros(Nx+Ny-1)
mxe=sum(g*x)/float(sum(g))
mye=sum(h*y)/float(sum(h))
vxe=sum(g*x**2)/float(sum(g))-(sum(g*x)/float(sum(g)))**2
vye=sum(h*y**2)/float(sum(h))-(sum(h*y)/float(sum(h)))**2
for k in range(-(Nx-1),Ny):
  tau[k]=k*dt
  sum1=0
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*(x[i]-mxe)*(y[i+k]-mye)
    sum0+=g[i]*h[i+k]
  
  rhoyx[k]=sum1/float(sum0)/sqrt(vxe*vye)

plot(tau,rhoyx,'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
rhoyx=zeros(Nx+Ny-1)
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*x**2)/float(sum(g))-(sum(g*x)/float(sum(g)))**2
vye=sum(h*y**2)/float(sum(h))-(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
CEyx1=real(ifft(Eyx1))/float(dt)
X0=fft(append(g,zeros(Ny)))
Y0=fft(append(h,zeros(Nx)))
Eyx0=dt**2*conj(X0)*Y0
CEyx0=real(ifft(Eyx0))/float(dt)
for k in range(-(Nx-1),Ny):
  rhoyx[k]=CEyx1[k]/CEyx0[k]/sqrt(vxe*vye)

plot(tau,rhoyx,'o')
show()
 direkte Berechnung:
rhoyx=zeros(1,Nx+Ny-1);
tau=zeros(1,Nx+Ny-1);
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*x.^2)/sum(g)-(sum(g.*x)/sum(g))^2;
vye=sum(h.*y.^2)/sum(h)-(sum(h.*y)/sum(h))^2;
for k=-(Nx-1):Ny-1
tau(mod(k,Nx+Ny-1)+1)=k*dt;
sum1=0;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*(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,rhoyx,'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
rhoyx=zeros(1,Nx+Ny-1);
tau=circshift((-(Nx-1):Ny-1)*dt,[0;Ny]);
mxe=sum(g.*x)/sum(g);
mye=sum(h.*y)/sum(h);
vxe=sum(g.*x.^2)/sum(g)-(sum(g.*x)/sum(g))^2;
vye=sum(h.*y.^2)/sum(h)-(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;
CEyx1=real(ifft(Eyx1))/dt;
X0=fft([g,zeros(1,Ny)]);
Y0=fft([h,zeros(1,Nx)]);
Eyx0=dt^2*conj(X0).*Y0;
CEyx0=real(ifft(Eyx0))/dt;
for k=-(Nx-1):Ny-1
rhoyx(mod(k,Nx+Ny-1)+1)=CEyx1(mod(k,Nx+Ny)+1)/CEyx0(mod(k,Nx+Ny)+1)/sqrt(vxe*vye);
end
plot(tau,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)
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
  sum2=0
  sum3=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*(x[i]-mxe)*(y[i+k]-mye)
    sum2+=g[i]*h[i+k]*(x[i]-mxe)**2
    sum3+=g[i]*h[i+k]*(y[i+k]-mye)**2
  
  rhoyx[k]=sum1/sqrt(sum2*sum3)

plot(tau,rhoyx,'o')
show()
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
rhoyx=zeros(Nx+Ny-1)
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*(x-mxe)**2,zeros(Ny)))
QY=fft(append(h*(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=real(ifft(Eyx))/float(dt)
CEA=real(ifft(EA))/float(dt)
CEB=real(ifft(EB))/float(dt)
for k in range(-(Nx-1),Ny):
  rhoyx[k]=CEyx[k]/sqrt(CEA[k]*CEB[k])

plot(tau,rhoyx,'o')
show()
 direkte Berechnung:
rhoyx=zeros(1,Nx+Ny-1);
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;
sum2=0;
sum3=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*(x(i)-mxe)*(y(i+k)-mye);
sum2=sum2+g(i)*h(i+k)*(x(i)-mxe)^2;
sum3=sum3+g(i)*h(i+k)*(y(i+k)-mye)^2;
end
rhoyx(mod(k,Nx+Ny-1)+1)=sum1/sqrt(sum2*sum3);
end
plot(tau,rhoyx,'o')
Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem:
rhoyx=zeros(1,Nx+Ny-1);
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.*(x-mxe).^2,zeros(1,Ny)]);
QY=fft([h.*(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=real(ifft(Eyx))/dt;
CEA=real(ifft(EA))/dt;
CEB=real(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,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=real(ifft(E1))/real(ifft(E0))
K=200
K=minimum(2*minimum(Nx,Ny)-1,K)
Ryx=zeros(K)
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)
for k in range(-(K//2),(K-1)//2+1):
  sum1=0
  sum0=0
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=g[i]*h[i+k]*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=real(ifft(E1))./real(ifft(E0));
K=200;
K=min(2*min(Nx,Ny)-1,K);
Ryx=zeros(1,K);
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);
for k=-floor(K/2):floor((K-1)/2)
sum1=0;
sum0=0;
for i=max(1,1-k):min(Nx,Ny-k)
sum1=sum1+g(i)*h(i+k)*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')