Signal- und Messdatenverarbeitung

Codes für stochastisch und unabhängig abgetastete, reelle Leistungssignalpaare, 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 korrelierte und gegeneinander verschobene AR1-Prozess als Beispiel, mit Erwartungswerten verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)
Tx=10000.0
Ty=10100.0
Ti=10.0
drx=1.0
dry=1.5
mxs=3.0
mys=2.0
vxs=1.0
vys=2.0
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
tx=[]
ty=[]
x=[]
y=[]
te=0.0
ue=standard_normal()
ve=standard_normal()
while te<maximum(Tx,Ty):
  tp=exponential(1.0/(2*(drx+dry)))
  te+=tp
  phi=exp(-tp/Ti)
  theta=sqrt(1-phi**2)
  ue=ue*phi+normal(0,theta)
  ve=ve*phi+normal(0,theta)
  xe=c11*ue+c12*ve+mxs
  ye=c21*ue+c22*ve+mys
  if random()<drx/float(drx+dry):
    if (te<Tx) and (random()<1/(1+exp(-(xe-mxs)))):
      tx.append(te)
      x.append(xe)
    
  else:
    if (te<Ty) and (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()
Tx=10000;
Ty=10100;
Ti=10;
drx=1;
dry=1.5;
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;
tx=[];
ty=[];
x=[];
y=[];
te=0;
ue=randn();
ve=randn();
while te<max(Tx,Ty)
tp=-log(1-rand())/(2*(drx+dry));
te=te+tp;
phi=exp(-tp/Ti);
theta=sqrt(1-phi^2);
ue=ue*phi+theta*randn();
ve=ve*phi+theta*randn();
xe=c11*ue+c12*ve+mxs;
ye=c21*ue+c22*ve+mys;
if rand<drx/(drx+dry)
if (te<Tx) && (rand<1/(1+exp(-(xe-mxs))))
tx(end+1)=te;
x(end+1)=xe;
end
else
if (te<Ty) && (rand<1/(1+exp(-(ye-mys))))
ty(end+1)=te;
y(end+1)=ye;
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 Kreuzleistungsdichtespektrum über Slotkorrelation 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|<Δt/20sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<min(2Tx,2Ty)/Δt
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⌋
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
K=200
K=minimum(int(ceil(2*minimum(Tx,Ty)/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R1=zeros(K)
R0=zeros(K)
i=0
jb=i
while i<Nx:
  j=jb
  while (jb<Ny) and (ty[jb]<tx[i]+((K-1)//2+0.5)*dt):
    jb+=1
  
  j=jb-1
  while (j>=0) and (ty[j]>tx[i]-(K//2+0.5)*dt):
    k=int(round((ty[j]-tx[i])/float(dt)))
    R1[k]+=g[i]*h[j]*x[i]*y[j]
    R0[k]+=g[i]*h[j]
    j-=1
  
  i+=1

Ryx=R1/R0
plot(tau,Ryx,'o',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R1=zeros(1,K);
R0=zeros(1,K);
i=1;
jb=i;
while i<=Nx
j=jb;
while (jb<=Ny) && (ty(jb)<tx(i)+(fix((K-1)/2)+0.5)*dt)
jb=jb+1;
end
j=jb-1;
while (j>0) && (ty(j)>tx(i)-(fix(K/2)+0.5)*dt)
k=round((ty(j)-tx(i))/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);
j=j-1;
end
i=i+1;
end
Ryx=R1./R0;
plot(tau,Ryx,'o',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Slotkorrelation (mit lokaler Normierung) 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|<Δt/20sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<min(2Tx,2Ty)/Δt
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⌋
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
K=200
K=minimum(int(ceil(2*minimum(Tx,Ty)/float(dt)))-1,K)
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)
i=0
jb=i
while i<Nx:
  j=jb
  while (jb<Ny) and (ty[jb]<tx[i]+((K-1)//2+0.5)*dt):
    jb+=1
  
  j=jb-1
  while (j>=0) and (ty[j]>tx[i]-(K//2+0.5)*dt):
    k=int(round((ty[j]-tx[i])/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
    j-=1
  
  i+=1

Ryx=sqrt(vxe*vye)*R1/sqrt(R2*R3)+mxe*mye
plot(tau,Ryx,'o',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
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);
i=1;
jb=i;
while i<=Nx
j=jb;
while (jb<=Ny) && (ty(jb)<tx(i)+(fix((K-1)/2)+0.5)*dt)
jb=jb+1;
end
j=jb-1;
while (j>0) && (ty(j)>tx(i)-(fix(K/2)+0.5)*dt)
k=round((ty(j)-tx(i))/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;
j=j-1;
end
i=i+1;
end
Ryx=sqrt(vxe*vye)*R1./sqrt(R2.*R3)+mxe*mye;
plot(tau,Ryx,'o',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, ohne Normierung) Eyx,j=Eyx(fj)=TxTyXj*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≥⌈(Tx+Ty)/Δ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)=RE,yx,kmin(Tx,Ty,Tx+τk,Ty−τk)
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<min(2Tx,2Ty)/Δt
Syx,j=Syx(fj)=Δt⋅FFT{Ryx,k}=Δt⋅∑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(ceil((Tx+Ty)/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=Tx*Ty*conj(X)*Y/float(sum(g)*sum(h))
REyx=real(ifft(Eyx))/float(dt)
K=200
K=minimum(int(ceil(2*minimum(Tx,Ty)/float(dt)))-1,K)
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(minimum(Tx,minimum(Ty,minimum(Tx+tau[k],Ty-tau[k]))))

plot(tau,Ryx,'o',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=ceil((Tx+Ty)/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=Tx*Ty*conj(X).*Y/(sum(g)*sum(h));
REyx=real(ifft(Eyx))/dt;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
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)/min(min(Tx,Ty),min(Tx+tau(mod(k+K,K)+1),Ty-tau(mod(k+K,K)+1)));
end
plot(tau,Ryx,'o',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Normierung) Eyx,j=Eyx(fj)=TxTyXj*YjX0′*Y0′
Eyx,j′=Eyx′(fj)=TxTyXj′*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≥⌈(Tx+Ty)/Δ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<min(2Tx,2Ty)/Δt
Syx,j=Syx(fj)=Δt⋅FFT{Ryx,k}=Δt⋅∑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(ceil((Tx+Ty)/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=Tx*Ty*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=Tx*Ty*conj(Xp)**Yp/(conj(Xp[0])*Yp[0])
REyx=real(ifft(Eyx))/float(dt)
REyxp=real(ifft(Eyxp))/float(dt)
K=200
K=minimum(int(ceil(2*minimum(Tx,Ty)/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=ceil((Tx+Ty)/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=Tx*Ty*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=Tx*Ty*conj(Xp).*Yp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
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',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit lokaler Normierung) Eyx,j=Eyx(fj)=TxTyXj*YjX0′*Y0′
Eyx,j′=Eyx′(fj)=TxTyXj′′*Yj′X0′*Y0′
Eyx,j′′=Eyx′′(fj)(fj)=TxTyXj′*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≥⌈(Tx+Ty)/Δ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<min(2Tx,2Ty)/Δt
Syx,j=Syx(fj)=Δt⋅FFT{Ryx,k}=Δt⋅∑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(ceil((Tx+Ty)/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=Tx*Ty*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=Tx*Ty*conj(Xpp)*Yp/(conj(Xp[0])*Yp[0])
Eyxpp=Tx*Ty*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=200
K=minimum(int(ceil(2*minimum(Tx,Ty)/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'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=ceil((Tx+Ty)/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=Tx*Ty*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=Tx*Ty*conj(Xpp).*Yp/(conj(Xp(1))*Yp(1));
Eyxpp=Tx*Ty*conj(Xp).*Ypp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
REyxpp=real(ifft(Eyxpp))/dt;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
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',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
(keine Lomb-Scargle-Methode für Kreuzspektren in Matlab und Python)
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (ohne Normierung)
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
J=int(ceil((Tx+Ty)/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=Tx*Ty*conj(X)*Y/float(sum(g)*sum(h))
REyx=real(ifft(Eyx))/float(dt)
K=200
K=minimum(int(ceil(2*minminimum((Tx,Ty)/float(dt)))-1,K)
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(minimum(Tx,minimum(Ty,minimum(Tx+tau[k],Ty-tau[k]))))

plot(tau,Ryx,'o',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=ceil((Tx+Ty)/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=Tx*Ty*conj(X).*Y/(sum(g)*sum(h));
REyx=real(ifft(Eyx))/dt;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
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)/min(min(Tx,Ty),min(Tx+tau(mod(k+K,K)+1),Ty-tau(mod(k+K,K)+1)));
end
plot(tau,Ryx,'o',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit Normierung)
from numpy.fft import *
Nx=len(tx)
Ny=len(ty)
dt=1.0
J=int(ceil((Tx+Ty)/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=Tx*Ty*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=Tx*Ty*conj(Xp)**Yp/(conj(Xp[0])*Yp[0])
REyx=real(ifft(Eyx))/float(dt)
REyxp=real(ifft(Eyxp))/float(dt)
K=200
K=minimum(int(ceil(2*minminimum((Tx,Ty)/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'o')
show()
Nx=length(tx);
Ny=length(ty);
dt=1.0;
J=ceil((Tx+Ty)/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=Tx*Ty*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=Tx*Ty*conj(Xp).*Yp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
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',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit lokaler Normierung)
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(ceil((Tx+Ty)/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=Tx*Ty*conj(X)*Y/(conj(Xp[0])*Yp[0])
Eyxp=Tx*Ty*conj(Xpp)*Yp/(conj(Xp[0])*Yp[0])
Eyxpp=Tx*Ty*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=200
K=minimum(int(ceil(2*minminimum((Tx,Ty)/float(dt)))-1,K)
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',arange(-5*K,5*K)*dt/10.0,cc*sqrt(vxs*vys)*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs*mys)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
Syx=dt*fft(Ryx)
p1=1-dt/float(Ti)
loglog(f,real(Syx),'o',arange(1,5*K)/float(10*K*dt),cc*sqrt(vxs*vys)*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)),f,imag(Syx),'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=ceil((Tx+Ty)/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=Tx*Ty*conj(X).*Y/(conj(Xp(1))*Yp(1));
Eyxp=Tx*Ty*conj(Xpp).*Yp/(conj(Xp(1))*Yp(1));
Eyxpp=Tx*Ty*conj(Xp).*Ypp/(conj(Xp(1))*Yp(1));
REyx=real(ifft(Eyx))/dt;
REyxp=real(ifft(Eyxp))/dt;
REyxpp=real(ifft(Eyxpp))/dt;
K=200;
K=min(ceil(2*min(Tx,Ty)/dt)-1,K);
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',(-5*K:5*K)*dt/10,cc*sqrt(vxs*vys)*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs*mys)
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
Syx=dt*fft(Ryx);
p1=1-dt/Ti;
loglog(f,real(Syx),'o',(1:5*K)/(10*K*dt),cc*sqrt(vxs*vys)*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)),f,imag(Syx),'o')
(keine Gewichtung für Interpolation)