Signal- und Messdatenverarbeitung

Codes für stochastisch und gemeinsam abgetastete, komplexe, 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 Wertetripeln (ti,xi,yi) mit den Messzeitpunkten ti und den komplexen Messwerten xi und yi (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
dr=5.0
mxs=3.0+1.5j
mys=2.0+0.5j
axs=1.5
ays=2.0
t=[]
x=[]
y=[]
te=0.0
while te<T:
  tp=exponential(1.0/(2*dr))
  te+=tp
  xe=axs*exp(0.05j*pi*te)+mxs
  ye=ays*exp(0.05j*pi*te+0.5)+mys
  if (te<T) and (random()<1/(1+exp(-real(xe-mxs)-real(ye-mys)))):
    t.append(te)
    x.append(xe)
    y.append(ye)
  

t=array(t)
x=array(x)
y=array(y)
plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')
show()
T=1000;
dr=5;
mxs=3+1.5j;
mys=2+0.5j;
axs=1.5;
ays=2;
t=[];
x=[];
y=[];
te=0;
while te<T
tp=-log(1-rand())/(2*dr);
te=te+tp;
xe=axs*exp(0.05j*pi*te)+mxs;
ye=ays*exp(0.05j*pi*te+0.5)+mys;
if (te<T) && (rand<1/(1+exp(-real(xe-mxs)-real(ye-mys))))
t(end+1)=te;
x(end+1)=xe;
y(end+1)=ye;
end
end
plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')
Generieren von reellen Gewichten gi passend zu den (ti,xi,yi) (Inverse der verschobenen Sigmoidfunktion)
g=1+exp(-real(x-mxs)-real(y-mys))
plot(t,g,'o')
show()
g=1+exp(-real(x-mxs)-real(y-mys));
plot(t,g,'o')
Mittelwerte x‾=∑i=0N−1gixi∑i=0N−1gi
y‾=∑i=0N−1giyi∑i=0N−1gi
mxe=sum(g*x)/float(sum(g))
mye=sum(g*y)/float(sum(g))
mxe=sum(g.*x)/sum(g)
mye=sum(g.*y)/sum(g)
Varianzen (ohne Bessel-Korrektur, asymptotisch erwartungstreu) sx2=∑i=0N−1gi(xi−x‾)*(xi−x‾)∑i=0N−1gi=∑i=0N−1gi|xi−x‾|2∑i=0N−1gi
=∑i=0N−1gixi*xi∑i=0N−1gi−x‾*x‾=∑i=0N−1gi|xi|2∑i=0N−1gi−|x‾|2
sy2=∑i=0N−1gi(yi−y‾)*(yi−y‾)∑i=0N−1gi=∑i=0N−1gi|yi−y‾|2∑i=0N−1gi
=∑i=0N−1giyi*yi∑i=0N−1gi−y‾*y‾=∑i=0N−1gi|yi|2∑i=0N−1gi−|y‾|2
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(g*abs(y)**2)/float(sum(g))-abs(sum(g*y)/float(sum(g)))**2
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2
vye=sum(g.*abs(y).^2)/sum(g)-abs(sum(g.*y)/sum(g))^2
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Slotkorrelation (ohne koinzidente Produkte, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Ryx,k=Ryx(τk)=Rxy*(−τk)=∑i=0N−1∑j=0N−1bk(tj−ti)gigjxi*yj∑i=0N−1∑j=0N−1bk(tj−ti)gigj
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 *
N=len(t)
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,N):
  for j in range(0,N):
    if j!=i:
      k=int(round((t[j]-t[i]-T*round((t[j]-t[i])/float(T)))/float(dt)))
      R1[k]+=g[i]*g[j]*conj(x[i])*y[j]
      R0[k]+=g[i]*g[j]
    
  

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()
N=length(t);
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:N
for j=1:N
if i~=j
k=round((t(j)-t(i)-T*round((t(j)-t(i))/T))/dt);
R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)*g(j)*conj(x(i))*y(j);
R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+g(i)*g(j);
end
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 (ohne koinzidente Produkte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Ryx,k=Ryx(τk)=Rxy*(−τk)=[∑i=0N−1∑j=0N−1bk(tj−ti)gigj(xi−x‾)*(yj−y‾)]⋅sx2sy2[∑i=0N−1∑j=0N−1bk(tj−ti)gigj|xi−x‾|2][∑i=0N−1∑j=0N−1bk(tj−ti)gigj|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 *
N=len(t)
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(g*y)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(g*abs(y)**2)/float(sum(g))-abs(sum(g*y)/float(sum(g)))**2
R1=zeros(K)+0j
R2=zeros(K)
R3=zeros(K)
for i in range(0,N):
  for j in range(0,N):
    if j!=i:
      k=int(round((t[j]-t[i]-T*round((t[j]-t[i])/float(T)))/float(dt)))
      R1[k]+=g[i]*g[j]*conj(x[i]-mxe)*(y[j]-mye)
      R2[k]+=g[i]*g[j]*abs(x[i]-mxe)**2
      R3[k]+=g[i]*g[j]*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()
N=length(t);
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(g.*y)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(g.*abs(y).^2)/sum(g)-abs(sum(g.*y)/sum(g))^2;
R1=zeros(1,K)+0j;
R2=zeros(1,K);
R3=zeros(1,K);
for i=1:N
for j=1:N
if i~=j
k=round((t(j)-t(i)-T*round((t(j)-t(i))/T))/dt);
R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)*g(j)*conj(x(i)-mxe)*(y(j)-mye);
R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+g(i)*g(j)*abs(x(i)-mxe)^2;
R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+g(i)*g(j)*abs(y(j)-mye)^2;
end
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, mit Abzug der koinzidenten Produkte, ohne Normierung) Eyx,j=Eyx(fj)=T2Xj*Yj−∑i=0N−1gi2xi*yi(∑i=0N−1gi)2−∑i=0N−1gi2
Xj=∑i=0N−1gixi𝐞−2π𝐢fjti
Yj=∑i=0N−1giyi𝐞−2π𝐢fjti
(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 *
N=len(t)
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,N):
  X+=g[i]*x[i]*exp(-2j*pi*fp*t[i])
  Y+=g[i]*y[i]*exp(-2j*pi*fp*t[i])

Eyx=T**2*(conj(X)*Y-sum(g**2*conj(x)*y))/float(sum(g)**2-sum(g**2))
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()
N=length(t);
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:N
X=X+g(i)*x(i)*exp(-2j*pi*fp*t(i));
Y=Y+g(i)*y(i)*exp(-2j*pi*fp*t(i));
end
Eyx=T^2*(conj(X).*Y-sum(g.^2.*conj(x).*y))/(sum(g)^2-sum(g.^2));
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 Abzug der koinzidenten Produkte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Eyx,j=Eyx(fj)=T2Xj*Yj−∑i=0N−1gi2xi*yiX0′2−∑i=0N−1gi2
Eyx,j′=Eyx′(fj)=T2|Xj′|2−∑i=0N−1gi2X0′2−∑i=0N−1gi2
Xj=∑i=0N−1gixi𝐞−2π𝐢fjti
Yj=∑i=0N−1giyi𝐞−2π𝐢fjti
Xj′=∑i=0N−1gi𝐞−2π𝐢fjti
(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 *
N=len(t)
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
for i in range(0,N):
  X+=g[i]*x[i]*exp(-2j*pi*fp*t[i])
  Y+=g[i]*y[i]*exp(-2j*pi*fp*t[i])
  Xp+=g[i]*exp(-2j*pi*fp*t[i])

Eyx=T**2*(conj(X)*Y-sum(g**2*conj(x)*y))/(abs(Xp[0])**2-sum(g**2))
Eyxp=T**2*(abs(Xp)**2-sum(g**2))/(abs(Xp[0])**2-sum(g**2))
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()
N=length(t);
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;
for i=1:N
X=X+g(i)*x(i)*exp(-2j*pi*fp*t(i));
Y=Y+g(i)*y(i)*exp(-2j*pi*fp*t(i));
Xp=Xp+g(i)*exp(-2j*pi*fp*t(i));
end
Eyx=T^2*(conj(X).*Y-sum(g.^2.*conj(x).*y))/(abs(Xp(1))^2-sum(g.^2));
Eyxp=T^2*(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2));
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 Abzug der koinzidenten Produkte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist) Eyx,j=Eyx(fj)=T2Xj*Yj−∑i=0N−1gi2(xi−x‾)*(yi−y‾)X0′2−∑i=0N−1gi2
Eyx,j′=Eyx′(fj)=T2Xj′′*Xj′−∑i=0N−1gi2|xi−x‾|2X0′2−∑i=0N−1gi2
Eyx,j′′=Eyx′′(fj)Eyx′(fj)=T2Xj′*Yj′′−∑i=0N−1gi2|yi−y‾|2X0′2−∑i=0N−1gi2
Xj=∑i=0N−1gi(xi−x‾)𝐞−2π𝐢fjti
Yj=∑i=0N−1gi(yi−y‾)𝐞−2π𝐢fjti
Xj′=∑i=0N−1gi𝐞−2π𝐢fjti
Xj′′=∑i=0N−1gi|xi−x‾|2𝐞−2π𝐢fjti
Yj′′=∑i=0N−1gi|yi−y‾|2𝐞−2π𝐢fjti
(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 *
N=len(t)
dt=1.0
mxe=sum(g*x)/float(sum(g))
mye=sum(g*y)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(g*abs(y)**2)/float(sum(g))-abs(sum(g*y)/float(sum(g)))**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
Xpp=zeros(size(fp))+0j
Ypp=zeros(size(fp))+0j
for i in range(0,N):
  X+=g[i]*(x[i]-mxe)*exp(-2j*pi*fp*t[i])
  Y+=g[i]*(y[i]-mye)*exp(-2j*pi*fp*t[i])
  Xp+=g[i]*exp(-2j*pi*fp*t[i])
  Xpp+=g[i]*abs(x[i]-mxe)**2*exp(-2j*pi*fp*t[i])
  Ypp+=g[i]*abs(y[i]-mxe)**2*exp(-2j*pi*fp*t[i])

Eyx=T**2*(conj(X)*Y-sum((g*conj(x-mxe)*(y-mye)))/(abs(Xp[0])**2-sum(g**2))
Eyxp=T**2*(conj(Xpp)*Xp-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
Eyxpp=T**2*(conj(Xp)*Ypp-sum((g*abs(y-mye))**2))/(abs(Xp[0])**2-sum(g**2))
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()
N=length(t);
dt=1.0;
mxe=sum(g.*x)/sum(g);
mye=sum(g.*y)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(g.*abs(y).^2)/sum(g)-abs(sum(g.*y)/sum(g))^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;
Xpp=zeros(size(fp))+0j;
Ypp=zeros(size(fp))+0j;
for i=1:N
X=X+g(i)*(x(i)-mxe)*exp(-2j*pi*fp*t(i));
Y=Y+g(i)*(y(i)-mye)*exp(-2j*pi*fp*t(i));
Xp=Xp+g(i)*exp(-2j*pi*fp*t(i));
Xpp=Xpp+g(i)*abs(x(i)-mxe)^2*exp(-2j*pi*fp*t(i));
Ypp=Ypp+g(i)*abs(y(i)-mye)^2*exp(-2j*pi*fp*t(i));
end
Eyx=T^2*(conj(X).*Y-sum(g.^2.*conj(x-mxe).*(y-mye)))/(abs(Xp(1))^2-sum(g.^2));
Eyxp=T^2*(conj(Xpp).*Xp-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
Eyxpp=T^2*(conj(Xp).*Ypp-sum((g.*abs(y-mye)).^2))/(abs(Xp(1))^2-sum(g.^2));
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 (mit Abzug der koinzidenten Produkte, ohne Normierung)
from numpy.fft import *
N=len(t)
dt=1.0
J=int(round(T/float(dt)))
x1=zeros(J)+0j
y1=zeros(J)+0j
for i in range(0,N):
  j=int(floor((t[i]-T*floor(t[i]/float(T)))/float(dt)))
  x1[j]+=g[i]*x[i];
  y1[j]+=g[i]*y[i];

X=fft(x1)
Y=fft(y1)
Eyx=T**2*(conj(X)*Y-sum(g**2*conj(x)*y))/float(sum(g)**2-sum(g**2))
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()
N=length(t);
dt=1.0;
J=round(T/dt);
x1=zeros([1,J])+0j;
y1=zeros([1,J])+0j;
for i=1:N
j=fix((t(i)-T*fix(t(i)/T))/dt)+1;
x1(j)=x1(j)+g(i)*x(i);
y1(j)=y1(j)+g(i)*y(i);
end
X=fft(x1);
Y=fft(y1);
Eyx=T^2*(conj(X).*Y-sum(g.^2.*conj(x).*y))/(sum(g)^2-sum(g.^2));
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 Abzug der koinzidenten Produkte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)
from numpy.fft import *
N=len(t)
dt=1.0
J=int(round(T/float(dt)))
x0=zeros(J)
x1=zeros(J)+0j
y1=zeros(J)+0j
for i in range(0,N):
  j=int(floor((t[i]-T*floor(t[i]/float(T)))/float(dt)))
  x0[j]+=g[i];
  x1[j]+=g[i]*x[i];
  y1[j]+=g[i]*y[i];

X=fft(x1)
Y=fft(y1)
Xp=fft(x0)
Eyx=T**2*(conj(X)*Y-sum(g**2*conj(x)*y))/(abs(Xp[0])**2-sum(g**2))
Eyxp=T**2*(abs(Xp)**2-sum(g**2))/(abs(Xp[0])**2-sum(g**2))
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()
N=length(t);
dt=1.0;
J=round(T/dt);
x0=zeros([1,J]);
x1=zeros([1,J])+0j;
y1=zeros([1,J])+0j;
for i=1:N
j=fix((t(i)-T*fix(t(i)/T))/dt)+1;
x0(j)=x0(j)+g(i);
x1(j)=x1(j)+g(i)*x(i);
y1(j)=y1(j)+g(i)*y(i);
end
X=fft(x1);
Y=fft(y1);
Xp=fft(x0);
Eyx=T^2*(conj(X).*Y-sum(g.^2.*conj(x).*y))/(abs(Xp(1))^2-sum(g.^2));
Eyxp=T^2*(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2));
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 Abzug der koinzidenten Produkte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)
from numpy.fft import *
N=len(t)
dt=1.0
mxe=sum(g*x)/float(sum(g))
mye=sum(g*y)/float(sum(g))
vxe=sum(g*abs(x)**2)/float(sum(g))-abs(sum(g*x)/float(sum(g)))**2
vye=sum(g*abs(y)**2)/float(sum(g))-abs(sum(g*y)/float(sum(g)))**2
J=int(round(T/float(dt)))
x0=zeros(J)
x1=zeros(J)+0j
y1=zeros(J)+0j
x2=zeros(J)
y2=zeros(J)
for i in range(0,N):
  j=int(floor((t[i]-T*floor(t[i]/float(T)))/float(dt)))
  x0[j]+=g[i];
  x1[j]+=g[i]*(x[i]-mxe);
  y1[j]+=g[i]*(y[i]-mye);
  x2[j]+=g[i]*abs(x[i]-mxe)**2;
  y2[j]+=g[i]*abs(y[i]-mye)**2;

X=fft(x1)
Y=fft(y1)
Xp=fft(x0)
Xpp=fft(x2)
Ypp=fft(y2)
Eyx=T**2*(conj(X)*Y-sum(g**2*conj(x-mxe)*(y-mye)))/(abs(Xp[0])**2-sum(g**2))
Eyxp=T**2*(conj(Xpp)*Xp-sum((g*abs(x-mxe))**2))/(abs(Xp[0])**2-sum(g**2))
Eyxpp=T**2*(conj(Xp)*Ypp-sum((g*abs(y-mye))**2))/(abs(Xp[0])**2-sum(g**2))
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()
N=length(t);
dt=1.0;
mxe=sum(g.*x)/sum(g);
mye=sum(g.*y)/sum(g);
vxe=sum(g.*abs(x).^2)/sum(g)-abs(sum(g.*x)/sum(g))^2;
vye=sum(g.*abs(y).^2)/sum(g)-abs(sum(g.*y)/sum(g))^2;
J=round(T/dt);
x0=zeros([1,J]);
x1=zeros([1,J])+0j;
y1=zeros([1,J])+0j;
x2=zeros([1,J]);
y2=zeros([1,J]);
for i=1:N
j=fix((t(i)-T*fix(t(i)/T))/dt)+1;
x0(j)=x0(j)+g(i);
x1(j)=x1(j)+g(i)*(x(i)-mxe);
y1(j)=y1(j)+g(i)*(y(i)-mye);
x2(j)=x2(j)+g(i)*abs(x(i)-mxe)^2;
y2(j)=y2(j)+g(i)*abs(y(i)-mye)^2;
end
X=fft(x1);
Y=fft(y1);
Xp=fft(x0);
Xpp=fft(x2);
Ypp=fft(y2);
Eyx=T^2*(conj(X).*Y-sum(g.^2.*conj(x-mxe).*(y-mye)))/(abs(Xp(1))^2-sum(g.^2));
Eyxp=T^2*(conj(Xpp).*Xp-sum((g.*abs(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2));
Eyxpp=T^2*(conj(Xp).*Ypp-sum((g.*abs(y-mye)).^2))/(abs(Xp(1))^2-sum(g.^2));
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 Gewichtung für Interpolation)