Signal- und Messdatenverarbeitung

Codes für stochastisch und unabhängig abgetastete, komplexe, periodische Signalpaare, ohne 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 komplexen Messwerten xi mit i=0…Nx−1 sowie den Wertepaaren (ty,i,yi) mit den Messzeitpunkten ty,i und den komplexen Messwerten yi mit i=0…Ny−1 (zwei harmonische Signale als Beispiel, mit einem Erwartungswert verschieden von null)
T=1000.0
drx=4.0
dry=6.0
mxs=3.0+1.5j
mys=2.0+0.5j
axs=1.5
ays=2.0
tx=[]
ty=[]
x=[]
y=[]
te=0.0
while te<T:
  tp=exponential(1.0/(drx+dry))
  te+=tp
  if te<T:
    if random()<drx/float(drx+dry):
      tx.append(te)
      x.append(axs*exp(0.05j*pi*te)+mxs)
    else:
      ty.append(te)
      y.append(ays*exp(0.05j*pi*te+0.5)+mys)
    
  

tx=array(tx)
ty=array(ty)
x=array(x)
y=array(y)
plot(tx,real(x),'o',tx,imag(x),'o',ty,real(y),'o',ty,imag(y),'o')
show()
T=1000;
drx=4;
dry=6;
mxs=3+1.5j;
mys=2+0.5j;
axs=1.5;
ays=2;
tx=[];
ty=[];
x=[];
y=[];
te=0;
while te<T
tp=-log(1-rand())/(drx+dry);
te=te+tp;
if te<T
if rand<drx/(drx+dry)
tx(end+1)=te;
x(end+1)=axs*exp(0.05j*pi*te)+mxs;
else
ty(end+1)=te;
y(end+1)=ays*exp(0.05j*pi*te+0.5)+mys;
end
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)
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)xi*yj∑i=0Nx−1∑j=0Ny−1bk(tj−ti)
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)+0j
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]+=conj(x[i])*y[j]
    R0[k]+=1.0
  

Ryx=R1/R0
plot(tau,real(Ryx),'o',tau,imag(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)+0j;
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)+conj(x(i))*y(j);
R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+1;
end
end
Ryx=R1./R0;
plot(tau,real(Ryx),'o',tau,imag(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)(xi−x‾)*(yj−y‾)]⋅sx2sy2[∑i=0Nx−1∑j=0Ny−1bk(tj−ti)|xi−x‾|2][∑i=0Nx−1∑j=0Ny−1bk(tj−ti)|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=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
R1=zeros(K)+0j
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]+=conj(x[i]-mxe)*(y[j]-mye)
    R2[k]+=abs(x[i]-mxe)**2
    R3[k]+=abs(y[j]-mye)**2
  

Ryx=sqrt(vxe*vye)*R1/sqrt(R2*R3)+conj(mxe)*mye
plot(tau,real(Ryx),'o',tau,imag(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=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
R1=zeros(1,K)+0j;
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)+conj(x(i)-mxe)*(y(j)-mye);
R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+abs(x(i)-mxe)^2;
R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+abs(y(j)-mye)^2;
end
end
Ryx=sqrt(vxe*vye)*R1./sqrt(R2.*R3)+conj(mxe)*mye;
plot(tau,real(Ryx),'o',tau,imag(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*YjNxNy
Xj=∑i=0Nx−1xi𝐞−2π𝐢fjtx,i
Yj=∑i=0Ny−1yi𝐞−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+=x[i]*exp(-2j*pi*fp*tx[i])

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

Eyx=T**2*conj(X)*Y/float(Nx*Ny)
REyx=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)+0j
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(T)

plot(tau,real(Ryx),'o',tau,imag(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+x(i)*exp(-2j*pi*fp*tx(i));
end
for i=1:Ny
Y=Y+y(i)*exp(-2j*pi*fp*ty(i));
end
Eyx=T^2*conj(X).*Y/(Nx*Ny);
REyx=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])+0j;
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,real(Ryx),'o',tau,imag(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*YjNxNy
Eyx,j′=Eyx′(fj)=T2Xj′*Yj′NxNy
Xj=∑i=0Nx−1xi𝐞−2π𝐢fjtx,i
Yj=∑i=0Ny−1yi𝐞−2π𝐢fjty,i
Xj′=∑i=0Nx−1𝐞−2π𝐢fjtx,i
Yj′=∑i=0Ny−1𝐞−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+=x[i]*exp(-2j*pi*fp*tx[i])
  Xp+=exp(-2j*pi*fp*tx[i])

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

Eyx=T**2*conj(X)*Y/float(Nx*Ny)
Eyxp=T**2*conj(Xp)*Yp/float(Nx*Ny)
REyx=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)+0j
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(REyxp[k])

plot(tau,real(Ryx),'o',tau,imag(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+x(i)*exp(-2j*pi*fp*tx(i));
Xp=Xp+exp(-2j*pi*fp*tx(i));
end
for i=1:Ny
Y=Y+y(i)*exp(-2j*pi*fp*ty(i));
Yp=Yp+exp(-2j*pi*fp*ty(i));
end
Eyx=T^2*conj(X).*Y/(Nx*Ny);
Eyxp=T^2*conj(Xp).*Yp/(Nx*Ny);
REyx=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])+0j;
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,real(Ryx),'o',tau,imag(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*YjNxNy
Eyx,j′=Eyx′(fj)=T2Xj′′*Yj′NxNy
Eyx,j′′=Eyx′′(fj)=T2Xj′*Yj′′NxNy
Xj=∑i=0Nx−1(xi−x‾)𝐞−2π𝐢fjtx,i
Yj=∑i=0Ny−1(yi−y‾)𝐞−2π𝐢fjty,i
Xj′=∑i=0Nx−1𝐞−2π𝐢fjtx,i
Yj′=∑i=0Ny−1𝐞−2π𝐢fjty,i
Xj′′=∑i=0Nx−1|xi−x‾|2𝐞−2π𝐢fjtx,i
Yj′′=∑i=0Ny−1|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=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
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+=(x[i]-mxe)*exp(-2j*pi*fp*tx[i])
  Xp+=exp(-2j*pi*fp*tx[i])
  Xpp+=abs(x[i]-mxe)**2*exp(-2j*pi*fp*tx[i])

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

Eyx=T**2*conj(X)*Y/float(Nx*Ny)
Eyxp=T**2*conj(Xpp)*Yp/float(Nx*Ny)
Eyxpp=T**2*conj(Xp)*Ypp/float(Nx*Ny)
REyx=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)+0j
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=sqrt(vxe*vye)*REyx[k]/sqrt(REyxp[k]*REyxpp[k])+conj(mxe)*mye

plot(tau,real(Ryx),'o',tau,imag(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=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
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+(x(i)-mxe)*exp(-2j*pi*fp*tx(i));
Xp=Xp+exp(-2j*pi*fp*tx(i));
Xpp=Xpp+abs(x(i)-mxe)^2*exp(-2j*pi*fp*tx(i));
end
for i=1:Ny
Y=Y+(y(i)-mye)*exp(-2j*pi*fp*ty(i));
Yp=Yp+exp(-2j*pi*fp*ty(i));
Ypp=Ypp+abs(y(i)-mye)^2*exp(-2j*pi*fp*ty(i));
end
Eyx=T^2*conj(X).*Y/(Nx*Ny);
Eyxp=T^2*conj(Xpp).*Yp/(Nx*Ny);
Eyxpp=T^2*conj(Xp).*Ypp/(Nx*Ny);
REyx=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])+0j;
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))+conj(mxe)*mye;
end
plot(tau,real(Ryx),'o',tau,imag(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)+0j
y1=zeros(J)+0j
for i in range(0,Nx):
  j=int(floor((tx[i]-T*floor(tx[i]/float(T)))/float(dt)))
  x1[j]+=x[i];

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

X=fft(x1)
Y=fft(y1)
Eyx=T**2*conj(X)*Y/float(Nx*Ny)
REyx=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)+0j
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(T)

plot(tau,real(Ryx),'o',tau,imag(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])+0j;
y1=zeros([1,J])+0j;
for i=1:Nx
j=fix((tx(i)-T*fix(tx(i)/T))/dt)+1;
x1(j)=x1(j)+x(i);
end
for i=1:Ny
j=fix((ty(i)-T*fix(ty(i)/T))/dt)+1;
y1(j)=y1(j)+y(i);
end
X=fft(x1);
Y=fft(y1);
Eyx=T^2*conj(X).*Y/(Nx*Ny);
REyx=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])+0j;
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,real(Ryx),'o',tau,imag(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)+0j
y1=zeros(J)+0j
for i in range(0,Nx):
  j=int(floor((tx[i]-T*floor(tx[i]/float(T)))/float(dt)))
  x0[j]+=1.0;
  x1[j]+=x[i];

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

X=fft(x1)
Y=fft(y1)
Xp=fft(x0)
Yp=fft(y0)
Eyx=T**2*conj(X)*Y/float(Nx*Ny)
Eyxp=T**2*conj(Xp)*Yp/float(Nx*Ny)
REyx=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)+0j
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=REyx[k]/float(REyxp[k])

plot(tau,real(Ryx),'o',tau,imag(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])+0j;
y1=zeros([1,J])+0j;
for i=1:Nx
j=fix((tx(i)-T*fix(tx(i)/T))/dt)+1;
x0(j)=x0(j)+1;
x1(j)=x1(j)+x(i);
end
for i=1:Ny
j=fix((ty(i)-T*fix(ty(i)/T))/dt)+1;
y0(j)=y0(j)+1;
y1(j)=y1(j)+y(i);
end
X=fft(x1);
Y=fft(y1);
Xp=fft(x0);
Yp=fft(y0);
Eyx=T^2*conj(X).*Y/(Nx*Ny);
Eyxp=T^2*conj(Xp).*Yp/(Nx*Ny);
REyx=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])+0j;
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,real(Ryx),'o',tau,imag(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 Normierun, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall istg)
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
J=int(round(T/float(dt)))
x0=zeros(J)
y0=zeros(J)
x1=zeros(J)+0j
y1=zeros(J)+0j
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]+=1.0;
  x1[j]+=x[i]-mxe;
  x2[j]+=abs(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]+=1.0;
  y1[j]+=y[i]-mye;
  y2[j]+=abs(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/float(Nx*Ny)
Eyxp=T**2*conj(Xpp)*Yp/float(Nx*Ny)
Eyxpp=T**2*conj(Xp)*Ypp/float(Nx*Ny)
REyx=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)+0j
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=sqrt(vxe*vye)*REyx[k]/sqrt(REyxp[k]*REyxpp[k])+conj(mxe)*mye

plot(tau,real(Ryx),'o',tau,imag(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=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
J=round(T/dt);
x0=zeros([1,J]);
y0=zeros([1,J]);
x1=zeros([1,J])+0j;
y1=zeros([1,J])+0j;
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)+1;
x1(j)=x1(j)+x(i)-mxe;
x2(j)=x2(j)+abs(x(i)-mxe)^2;
end
for i=1:Ny
j=fix((ty(i)-T*fix(ty(i)/T))/dt)+1;
y0(j)=y0(j)+1;
y1(j)=y1(j)+y(i)-mye;
y2(j)=y2(j)+abs(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/(Nx*Ny);
Eyxp=T^2*conj(Xpp).*Yp/(Nx*Ny);
Eyxpp=T^2*conj(Xp).*Ypp/(Nx*Ny);
REyx=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])+0j;
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))+conj(mxe)*mye;
end
plot(tau,real(Ryx),'o',tau,imag(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 Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
J=int(round(T/float(dt)))
xp=x[Nx-1]*ones(J)
yp=y[Ny-1]*ones(J)
for i in range(0,Nx-1):
  for j in range(int(floor(tx[i]/float(dt)))+1,int(floor(tx[i+1]/float(dt)))+1):
    if j<J:
      xp[j]=x[i]
    
  

for i in range(0,Ny-1):
  for j in range(int(floor(ty[i]/float(dt)))+1,int(floor(ty[i+1]/float(dt)))+1):
    if j<J:
      yp[j]=y[i]
    
  

Xp=fft(xp)
Yp=fft(yp)
Eyxp=dt**2*conj(Xp)*Yp
REyxp=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)+0j
cx=exp(-ddx/float(dt))/((1-exp(-ddx/float(dt)))*(1-exp(-ddy/float(dt))))
cy=exp(-ddy/float(dt))/((1-exp(-ddx/float(dt)))*(1-exp(-ddy/float(dt))))
for k in range(-(K//2),(K+1)//2):
  Ryx[k]=((cx+cy+1)*REyxp[k]-cy*REyxp[k-1]-cx*REyxp[k+1])/float(T)

plot(tau,real(Ryx),'o',tau,imag(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);
xp=x(Nx)*ones([1,J]);
yp=y(Ny)*ones([1,J]);
for i=1:Nx-1
for j=fix(tx(i)/dt)+2:fix(tx(i+1)/dt)+1
if j<=J
xp(j)=x(i);
end
end
end
for i=1:Ny-1
for j=fix(ty(i)/dt)+2:fix(ty(i+1)/dt)+1
if j<=J
yp(j)=y(i);
end
end
end
Xp=fft(xp);
Yp=fft(yp);
Eyxp=dt^2*conj(Xp).*Yp;
REyxp=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])+0j;
cx=exp(-ddx/dt)/((1-exp(-ddx/dt))*(1-exp(-ddy/dt)));
cy=exp(-ddy/dt)/((1-exp(-ddx/dt))*(1-exp(-ddy/dt)));
for k=-fix(K/2):fix((K-1)/2)
Ryx(mod(k+K,K)+1)=((cx+cy+1)*REyxp(mod(k+length(REyxp),length(REyxp))+1)-cy*REyxp(mod(k-1+length(REyxp),length(REyxp))+1)-cx*REyxp(mod(k+1+length(REyxp),length(REyxp))+1))/T;
end
plot(tau,real(Ryx),'o',tau,imag(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')