Signal- und Messdatenverarbeitung

Codes für stochastisch und unabhängig abgetastete, reelle, periodische Signalpaare, mit Gewichtung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren von Wertepaaren (tx,i,xi) mit den Messzeitpunkten tx,i und den Messwerten xi mit i=0…Nx−1 sowie den Wertepaaren (ty,i,yi) mit den Messzeitpunkten ty,i und den Messwerten yi mit i=0…Ny−1 (zwei harmonische Signale als Beispiel, mit einem Erwartungswert verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)
T=1000.0
drx=4.0
dry=6.0
mxs=3.0
mys=2.0
axs=1.5
ays=2.0
tx=[]
ty=[]
x=[]
y=[]
te=0.0
while te<T:
  tp=exponential(1.0/(2*(drx+dry)))
  te+=tp
  if te<T:
    if random()<drx/float(drx+dry):
      xe=axs*sin(0.05*pi*te)+mxs
      if random()<1/(1+exp(-(xe-mxs))):
          tx.append(te)
          x.append(xe)
    else:
      ye=ays*sin(0.05*pi*te+0.5)+mys
      if random()<1/(1+exp(-(ye-mys))):
          ty.append(te)
          y.append(ye)
        
      
    
  

tx=array(tx)
ty=array(ty)
x=array(x)
y=array(y)
plot(tx,x,'o',ty,y,'o')
show()
T=1000;
drx=4;
dry=6;
mxs=3;
mys=2;
axs=1.5;
ays=2;
tx=[];
ty=[];
x=[];
y=[];
te=0;
while te<T
tp=-log(1-rand())/(2*(drx+dry));
te=te+tp;
if te<T
if rand<drx/(drx+dry)
xe=axs*sin(0.05*pi*te)+mxs;
if rand<1/(1+exp(-(xe-mxs)))
tx(end+1)=te;
x(end+1)=xe;
end
else
ye=ays*sin(0.05*pi*te+0.5)+mys;
if rand<1/(1+exp(-(ye-mys)))
ty(end+1)=te;
y(end+1)=ye;
end
end
end
end
plot(tx,x,'o',ty,y,'o')
Generieren von Gewichten gi passend zu den (tx,i,xi) und hi passend zu den (ty,i,yi) (Inverse der verschobenen Sigmoidfunktion)
g=1+exp(-(x-mxs))
h=1+exp(-(y-mys))
plot(tx,g,'o',ty,h,'o')
show()
g=1+exp(-(x-mxs));
h=1+exp(-(y-mys));
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
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Slotkorrelation (auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Ryx,k=Ryx(τk)=Rxy(−τk)=∑i=0Nx−1∑j=0Ny−1bk(tj−ti)gihjxiyj∑i=0Nx−1∑j=0Ny−1bk(tj−ti)gihj
bk(t)={1falls 0<|t+cT|<Δt/2 für irgendein ganzes c0sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=[T/Δt]
Pyx,j=Pyx(fj)=1K⋅FFT{Ryx,k}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Ryx,k𝐞−2π𝐢jk/K
(imaginäre Einheit 𝐢)
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R1=zeros(K)
R0=zeros(K)
for i in range(0,Nx):
  for j in range(0,Ny):
    k=int(round((ty[j]-tx[i]-T*round((ty[j]-tx[i])/float(T)))/float(dt)))
    R1[k]+=g[i]*h[j]*x[i]*y[j]
    R0[k]+=g[i]*h[j]
  

Ryx=R1/R0
plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R1=zeros(1,K);
R0=zeros(1,K);
for i=1:Nx
for j=1:Ny
k=round((ty(j)-tx(i)-T*round((ty(j)-tx(i))/T))/dt);
R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)*h(j)*x(i)*y(j);
R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+g(i)*h(j);
end
end
Ryx=R1./R0;
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Slotkorrelation (mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Ryx,k=Ryx(τk)=Rxy(−τk)=[∑i=0Nx−1∑j=0Ny−1bk(tj−ti)gihj(xi−x‾)(yj−y‾)]⋅sx2sy2[∑i=0Nx−1∑j=0Ny−1bk(tj−ti)gihj(xi−x‾)2][∑i=0Nx−1∑j=0Ny−1bk(tj−ti)gihj(yj−y‾)2]+x‾⋅y‾
bk(t)={1falls 0<|t+cT|<Δt/2 für irgendein ganzes c0sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=[T/Δt]
Pyx,j=Pyx(fj)=1K⋅FFT{Ryx,k}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Ryx,k𝐞−2π𝐢jk/K
(imaginäre Einheit 𝐢)
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
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
R1=zeros(K)
R2=zeros(K)
R3=zeros(K)
for i in range(0,Nx):
  for j in range(0,Ny):
    k=int(round((ty[j]-tx[i]-T*round((ty[j]-tx[i])/float(T)))/float(dt)))
    R1[k]+=g[i]*h[j]*(x[i]-mxe)*(y[j]-mye)
    R2[k]+=g[i]*h[j]*(x[i]-mxe)**2
    R3[k]+=g[i]*h[j]*(y[j]-mye)**2
  

Ryx=sqrt(vxe*vye)*R1/sqrt(R2*R3)+mxe*mye
plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
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;
R1=zeros(1,K);
R2=zeros(1,K);
R3=zeros(1,K);
for i=1:Nx
for j=1:Ny
k=round((ty(j)-tx(i)-T*round((ty(j)-tx(i))/T))/dt);
R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)*h(j)*(x(i)-mxe)*(y(j)-mye);
R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+g(i)*h(j)*(x(i)-mxe)^2;
R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+g(i)*h(j)*(y(j)-mye)^2;
end
end
Ryx=sqrt(vxe*vye)*R1./sqrt(R2.*R3)+mxe*mye;
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, ohne Normierung) Eyx,j=Eyx(fj)=T2Xj*Yj(∑i=0Nx−1gi)(∑i=0Ny−1hi)
Xj=∑i=0Nx−1gixi𝐞−2π𝐢fjtx,i
Yj=∑i=0Ny−1hiyi𝐞−2π𝐢fjty,i
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J=[T/Δt]
RE,yx,k=RE,yx(τk)=1Δt⋅IFFT{Eyx,j}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Eyx,j𝐞2π𝐢fjτk/J
Ryx,k=Ryx(τk)=Rxy(−τk)=1TRE,yx,k
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=J=[T/Δt]
Pyx,j=Pyx(fj)=1K⋅FFT{Ryx,k}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Ryx,k𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
J=int(round(T/float(dt)))
fp=roll(arange(-(J//2),(J+1)//2)/float(J*dt),(J+1)//2)
X=zeros(size(fp))+0j
Y=zeros(size(fp))+0j
for i in range(0,Nx):
  X+=g[i]*x[i]*exp(-2j*pi*fp*tx[i])

for i in range(0,Ny):
  Y+=h[i]*y[i]*exp(-2j*pi*fp*ty[i])

Eyx=T**2*conj(X)*Y/float(sum(g)*sum(h))
REyx=real(ifft(Eyx))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
Ryx=zeros(K)
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(T)

plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=round(T/dt);
fp=circshift((-fix(J/2):fix((J-1)/2))/(J*dt),[0;fix((J+1)/2)]);
X=zeros(size(fp))+0j;
Y=zeros(size(fp))+0j;
for i=1:Nx
X=X+g(i)*x(i)*exp(-2j*pi*fp*tx(i));
end
for i=1:Ny
Y=Y+h(i)*y(i)*exp(-2j*pi*fp*ty(i));
end
Eyx=T^2*conj(X).*Y/(sum(g)*sum(h));
REyx=real(ifft(Eyx))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
Ryx=zeros([1,K]);
for k=-fix(K/2):fix((K-1)/2)
Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/T;
end
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Eyx,j=Eyx(fj)=T2Xj*YjX0′*Y0′
Eyx,j′=Eyx′(fj)=T2Xj′*Yj′X0′*Y0′
Xj=∑i=0Nx−1gixi𝐞−2π𝐢fjtx,i
Yj=∑i=0Ny−1hiyi𝐞−2π𝐢fjty,i
Xj′=∑i=0Nx−1gi𝐞−2π𝐢fjtx,i
Yj′=∑i=0Ny−1hi𝐞−2π𝐢fjty,i
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J=[T/Δt]
RE,yx,k=RE,yx(τk)=1Δt⋅IFFT{Eyx,j}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Eyx,j𝐞2π𝐢fjτk/J
RE,yx,k′=RE,yx′(τk)=1Δt⋅IFFT{Eyx,j′}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Eyx,j′𝐞2π𝐢fjτk/J
Ryx,k=Ryx(τk)=Rxy(−τk)=RE,yx,kRE,yx,k′
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=J=[T/Δt]
Pyx,j=Pyx(fj)=1K⋅FFT{Ryx,k}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Ryx,k𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
J=int(round(T/float(dt)))
fp=roll(arange(-(J//2),(J+1)//2)/float(J*dt),(J+1)//2)
X=zeros(size(fp))+0j
Y=zeros(size(fp))+0j
Xp=zeros(size(fp))+0j
Yp=zeros(size(fp))+0j
for i in range(0,Nx):
  X+=g[i]*x[i]*exp(-2j*pi*fp*tx[i])
  Xp+=g[i]*exp(-2j*pi*fp*tx[i])

for i in range(0,Ny):
  Y+=h[i]*y[i]*exp(-2j*pi*fp*ty[i])
  Yp+=h[i]*exp(-2j*pi*fp*ty[i])

Eyx=T**2*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=T**2*conj(Xp)**Yp/(conj(Xp[0])*Yp[0])
REyx=real(ifft(Eyx))/float(dt)
REyxp=real(ifft(Eyxp))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
Ryx=zeros(K)
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(REyxp[k])

plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=round(T/dt);
fp=circshift((-fix(J/2):fix((J-1)/2))/(J*dt),[0;fix((J+1)/2)]);
X=zeros(size(fp))+0j;
Y=zeros(size(fp))+0j;
Xp=zeros(size(fp))+0j;
Yp=zeros(size(fp))+0j;
for i=1:Nx
X=X+g(i)*x(i)*exp(-2j*pi*fp*tx(i));
Xp=Xp+g(i)*exp(-2j*pi*fp*tx(i));
end
for i=1:Ny
Y=Y+h(i)*y(i)*exp(-2j*pi*fp*ty(i));
Yp=Yp+h(i)*exp(-2j*pi*fp*ty(i));
end
Eyx=T^2*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=T^2*conj(Xp).*Yp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
Ryx=zeros([1,K]);
for k=-fix(K/2):fix((K-1)/2)
Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/REyxp(mod(k+length(REyxp),length(REyxp))+1);
end
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Eyx,j=Eyx(fj)=T2Xj*YjX0′*Y0′
Eyx,j′=Eyx′(fj)=T2Xj′′*Yj′X0′*Y0′
Eyx,j′′=Eyx′′(fj)=T2Xj′*Yj′′X0′*Y0′
Xj=∑i=0Nx−1gi(xi−x‾)𝐞−2π𝐢fjtx,i
Yj=∑i=0Ny−1hi(yi−y‾)𝐞−2π𝐢fjty,i
Xj′=∑i=0Nx−1gi𝐞−2π𝐢fjtx,i
Yj′=∑i=0Ny−1hi𝐞−2π𝐢fjty,i
Xj′′=∑i=0Nx−1gi(xi−x‾)2𝐞−2π𝐢fjtx,i
Yj′′=∑i=0Ny−1hi(yi−y‾)2𝐞−2π𝐢fjty,i
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J=[T/Δt]
RE,yx,k=RE,yx(τk)=1Δt⋅IFFT{Eyx,j}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Eyx,j𝐞2π𝐢fjτk/J
RE,yx,k′=RE,yx′(τk)=1Δt⋅IFFT{Eyx,j′}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Eyx,j′𝐞2π𝐢fjτk/J
RE,yx,k′′=RE′′(τk)=1Δt⋅IFFT{Eyx,j′′}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Eyx,j′′𝐞2π𝐢fjτk/J
Ryx,k=Ryx(τk)=Rxy(−τk)=RE,yx,ksx2sy2RE,yx,k′RE,yx,k′′+x‾⋅y‾
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K=J=[T/Δt]
Pyx,j=Pyx(fj)=1K⋅FFT{Ryx,k}=1K⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Ryx,k𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
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
J=int(round(T/float(dt)))
fp=roll(arange(-(J//2),(J+1)//2)/float(J*dt),(J+1)//2)
X=zeros(size(fp))+0j
Y=zeros(size(fp))+0j
Xp=zeros(size(fp))+0j
Yp=zeros(size(fp))+0j
Xpp=zeros(size(fp))+0j
Ypp=zeros(size(fp))+0j
for i in range(0,Nx):
  X+=g[i]*(x[i]-mxe)*exp(-2j*pi*fp*tx[i])
  Xp+=g[i]*exp(-2j*pi*fp*tx[i])
  Xpp+=g[i]*(x[i]-mxe)**2*exp(-2j*pi*fp*tx[i])

for i in range(0,Ny):
  Y+=h[i]*(y[i]-mye)*exp(-2j*pi*fp*ty[i])
  Yp+=h[i]*exp(-2j*pi*fp*ty[i])
  Ypp+=h[i]*(y[i]-mye)**2*exp(-2j*pi*fp*ty[i])

Eyx=T**2*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=T**2*conj(Xpp)*Yp/(conj(Xp[0])*Yp[0])
Eyxpp=T**2*conj(Xp)*Ypp/(conj(Xp[0])*Yp[0])
REyx=real(ifft(Eyx))/float(dt)
REyxp=real(ifft(Eyxp))/float(dt)
REyxpp=real(ifft(Eyxpp))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
Ryx=zeros(K)
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=sqrt(vxe*vye)*REyx[k]/sqrt(REyxp[k]*REyxpp[k])+mxe*mye

plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
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;
J=round(T/dt);
fp=circshift((-fix(J/2):fix((J-1)/2))/(J*dt),[0;fix((J+1)/2)]);
X=zeros(size(fp))+0j;
Y=zeros(size(fp))+0j;
Xp=zeros(size(fp))+0j;
Yp=zeros(size(fp))+0j;
Xpp=zeros(size(fp))+0j;
Ypp=zeros(size(fp))+0j;
for i=1:Nx
X=X+g(i)*(x(i)-mxe)*exp(-2j*pi*fp*tx(i));
Xp=Xp+g(i)*exp(-2j*pi*fp*tx(i));
Xpp=Xpp+g(i)*(x(i)-mxe)^2*exp(-2j*pi*fp*tx(i));
end
for i=1:Ny
Y=Y+h(i)*(y(i)-mye)*exp(-2j*pi*fp*ty(i));
Yp=Yp+h(i)*exp(-2j*pi*fp*ty(i));
Ypp=Ypp+h(i)*(y(i)-mye)^2*exp(-2j*pi*fp*ty(i));
end
Eyx=T^2*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=T^2*conj(Xpp).*Yp/(conj(Xp(1))*Yp(1));
Eyxpp=T^2*conj(Xp).*Ypp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
REyxpp=real(ifft(Eyxpp))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
Ryx=zeros([1,K]);
for k=-fix(K/2):fix((K-1)/2)
Ryx(mod(k+K,K)+1)=sqrt(vxe*vye)*REyx(mod(k+length(REyx),length(REyx))+1)/sqrt(REyxp(mod(k+length(REyxp),length(REyxp))+1)*REyxpp(mod(k+length(REyxpp),length(REyxpp))+1))+mxe*mye;
end
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
(keine Lomb-Scargle-Methode für Kreuzspektren in Matlab und Python)
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Zeitquantisierung (ohne Normierung)
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
J=int(round(T/float(dt)))
x1=zeros(J)
y1=zeros(J)
for i in range(0,Nx):
  j=int(floor((tx[i]-T*floor(tx[i]/float(T)))/float(dt)))
  x1[j]+=g[i]*x[i];

for i in range(0,Ny):
  j=int(floor((ty[i]-T*floor(ty[i]/float(T)))/float(dt)))
  y1[j]+=h[i]*y[i];

X=fft(x1)
Y=fft(y1)
Eyx=T**2*conj(X)*Y/float(sum(g)*sum(h))
REyx=real(ifft(Eyx))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
Ryx=zeros(K)
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(T)

plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=round(T/dt);
x1=zeros([1,J]);
y1=zeros([1,J]);
for i=1:Nx
j=fix((tx(i)-T*fix(tx(i)/T))/dt)+1;
x1(j)=x1(j)+g(i)*x(i);
end
for i=1:Ny
j=fix((ty(i)-T*fix(ty(i)/T))/dt)+1;
y1(j)=y1(j)+h(i)*y(i);
end
X=fft(x1);
Y=fft(y1);
Eyx=T^2*conj(X).*Y/(sum(g)*sum(h));
REyx=real(ifft(Eyx))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
Ryx=zeros([1,K]);
for k=-fix(K/2):fix((K-1)/2)
Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/T;
end
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Zeitquantisierung (mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
J=int(round(T/float(dt)))
x0=zeros(J)
y0=zeros(J)
x1=zeros(J)
y1=zeros(J)
for i in range(0,Nx):
  j=int(floor((tx[i]-T*floor(tx[i]/float(T)))/float(dt)))
  x0[j]+=g[i];
  x1[j]+=g[i]*x[i];

for i in range(0,Ny):
  j=int(floor((ty[i]-T*floor(ty[i]/float(T)))/float(dt)))
  y0[j]+=h[i];
  y1[j]+=h[i]*y[i];

X=fft(x1)
Y=fft(y1)
Xp=fft(x0)
Yp=fft(y0)
Eyx=T**2*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=T**2*conj(Xp)**Yp/(conj(Xp[0])*Yp[0])
REyx=real(ifft(Eyx))/float(dt)
REyxp=real(ifft(Eyxp))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
Ryx=zeros(K)
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(REyxp[k])

plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=round(T/dt);
x0=zeros([1,J]);
y0=zeros([1,J]);
x1=zeros([1,J]);
y1=zeros([1,J]);
for i=1:Nx
j=fix((tx(i)-T*fix(tx(i)/T))/dt)+1;
x0(j)=x0(j)+g(i);
x1(j)=x1(j)+g(i)*x(i);
end
for i=1:Ny
j=fix((ty(i)-T*fix(ty(i)/T))/dt)+1;
y0(j)=y0(j)+h(i);
y1(j)=y1(j)+h(i)*y(i);
end
X=fft(x1);
Y=fft(y1);
Xp=fft(x0);
Yp=fft(y0);
Eyx=T^2*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=T^2*conj(Xp).*Yp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
Ryx=zeros([1,K]);
for k=-fix(K/2):fix((K-1)/2)
Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/REyxp(mod(k+length(REyxp),length(REyxp))+1);
end
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Zeitquantisierung (mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
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
J=int(round(T/float(dt)))
x0=zeros(J)
y0=zeros(J)
x1=zeros(J)
y1=zeros(J)
x2=zeros(J)
y2=zeros(J)
for i in range(0,Nx):
  j=int(floor((tx[i]-T*floor(tx[i]/float(T)))/float(dt)))
  x0[j]+=g[i];
  x1[j]+=g[i]*(x[i]-mxe);
  x2[j]+=g[i]*(x[i]-mxe)**2;

for i in range(0,Ny):
  j=int(floor((ty[i]-T*floor(ty[i]/float(T)))/float(dt)))
  y0[j]+=h[i];
  y1[j]+=h[i]*(y[i]-mye);
  y2[j]+=h[i]*(y[i]-mye)**2;

X=fft(x1)
Y=fft(y1)
Xp=fft(x0)
Yp=fft(y0)
Xpp=fft(x2)
Ypp=fft(y2)
Eyx=T**2*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=T**2*conj(Xpp)*Yp/(conj(Xp[0])*Yp[0])
Eyxpp=T**2*conj(Xp)*Ypp/(conj(Xp[0])*Yp[0])
REyx=real(ifft(Eyx))/float(dt)
REyxp=real(ifft(Eyxp))/float(dt)
REyxpp=real(ifft(Eyxpp))/float(dt)
K=int(round(T/float(dt)))
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
Ryx=zeros(K)
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=sqrt(vxe*vye)*REyx[k]/sqrt(REyxp[k]*REyxpp[k])+mxe*mye

plot(tau,Ryx,'o')
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Pyx=fft(Ryx)/float(K)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
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;
J=round(T/dt);
x0=zeros([1,J]);
y0=zeros([1,J]);
x1=zeros([1,J]);
y1=zeros([1,J]);
x2=zeros([1,J]);
y2=zeros([1,J]);
for i=1:Nx
j=fix((tx(i)-T*fix(tx(i)/T))/dt)+1;
x0(j)=x0(j)+g(i);
x1(j)=x1(j)+g(i)*(x(i)-mxe);
x2(j)=x2(j)+g(i)*(x(i)-mxe)^2;
end
for i=1:Ny
j=fix((ty(i)-T*fix(ty(i)/T))/dt)+1;
y0(j)=y0(j)+h(i);
y1(j)=y1(j)+h(i)*(y(i)-mye);
y2(j)=y2(j)+h(i)*(y(i)-mye)^2;
end
X=fft(x1);
Y=fft(y1);
Xp=fft(x0);
Yp=fft(y0);
Xpp=fft(x2);
Ypp=fft(y2);
Eyx=T^2*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=T^2*conj(Xpp).*Yp/(conj(Xp(1))*Yp(1));
Eyxpp=T^2*conj(Xp).*Ypp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
REyxpp=real(ifft(Eyxpp))/dt;
K=round(T/dt);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
Ryx=zeros([1,K]);
for k=-fix(K/2):fix((K-1)/2)
Ryx(mod(k+K,K)+1)=sqrt(vxe*vye)*REyx(mod(k+length(REyx),length(REyx))+1)/sqrt(REyxp(mod(k+length(REyxp),length(REyxp))+1)*REyxpp(mod(k+length(REyxpp),length(REyxpp))+1))+mxe*mye;
end
plot(tau,Ryx,'o')
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Pyx=fft(Ryx)/K;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
(keine Gewichtung für Interpolation)