Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für stochastisch und gemeinsam abgetastete, komplexe Leistungssignalpaare, mit Gewichtung

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von Wertetripeln $(t_{i}, x_{i}, y_{i})$ mit den Messzeitpunkten $t_{i}$ und den komplexen Messwerten $x_{i}$ und $y_{i}$ (zwei korrelierte AR1-Prozess als Beispiel, mit Erwartungswerten verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)		T=10000.0 Ti=10.0 dr=1.0 mxs=3.0+1.5j mys=2.0+0.5j vxs=1.0 vys=2.0 cc=0.5 c1=sqrt(0.5+0.5sqrt(1-cc2/(vxsvys))) c2=cc/sqrt(vxsvys)/(2c1) c11=sqrt(vxs)c1 c12=sqrt(vxs)c2 c21=sqrt(vys)c2 c22=sqrt(vys)c1 t=[] x=[] y=[] te=0.0 ue=normal(0,sqrt(0.5))+1jnormal(0,sqrt(0.5)) ve=normal(0,sqrt(0.5))+1jnormal(0,sqrt(0.5)) while te<T: tp=exponential(1.0/(2dr)) te+=tp phi=exp(-tp/Ti) theta=sqrt(0.5(1-phi*2)) ue=uephi+normal(0,theta)+1jnormal(0,theta) ve=vephi+normal(0,theta)+1jnormal(0,theta) xe=c11ue+c12ve+mxs ye=c21ue+c22*ve+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=10000; Ti=10; dr=1; mxs=3+1.5j; mys=2+0.5j; vxs=1; vys=2; cc=0.5; c1=sqrt(0.5+0.5sqrt(1-cc^2/(vxsvys))); c2=cc/sqrt(vxsvys)/(2c1); c11=sqrt(vxs)c1; c12=sqrt(vxs)c2; c21=sqrt(vys)c2; c22=sqrt(vys)c1; t=[]; x=[]; y=[]; te=0; ue=sqrt(0.5)randn()+1jsqrt(0.5)randn(); ve=sqrt(0.5)randn()+1jsqrt(0.5)randn(); while te<T tp=-log(1-rand())/(2dr); te=te+tp; phi=exp(-tp/Ti); theta=sqrt(0.5(1-phi^2)); ue=uephi+thetarandn()+1jthetarandn(); ve=vephi+thetarandn()+1jthetarandn(); xe=c11ue+c12ve+mxs; ye=c21ue+c22ve+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 $g_{i}$ passend zu den $(t_{i}, x_{i}, y_{i})$ (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	$\overline{x} = \frac{\sum_{i = 0}^{N - 1} g_{i} x_{i}}{\sum_{i = 0}^{N - 1} g_{i}}$ $\overline{y} = \frac{\sum_{i = 0}^{N - 1} g_{i} y_{i}}{\sum_{i = 0}^{N - 1} g_{i}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g))`	`mxe=sum(g.x)/sum(g) mye=sum(g.y)/sum(g)`
Varianzen (ohne Bessel-Korrektur, asymptotisch erwartungstreu)	$s_{x}^{2} = \frac{\sum_{i = 0}^{N - 1} g_{i} {(x_{i} - \overline{x})}^{} (x_{i} - \overline{x})}{\sum_{i = 0}^{N - 1} g_{i}} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| x_{i} - \overline{x} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}}$ $= \frac{\sum_{i = 0}^{N - 1} g_{i} x_{i}^{} x_{i}}{\sum_{i = 0}^{N - 1} g_{i}} - {\overline{x}}^{} \overline{x} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| x_{i} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}} - {\| \overline{x} \|}^{2}$ $s_{y}^{2} = \frac{\sum_{i = 0}^{N - 1} g_{i} {(y_{i} - \overline{y})}^{} (y_{i} - \overline{y})}{\sum_{i = 0}^{N - 1} g_{i}} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| y_{i} - \overline{y} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}}$ $= \frac{\sum_{i = 0}^{N - 1} g_{i} y_{i}^{} y_{i}}{\sum_{i = 0}^{N - 1} g_{i}} - {\overline{y}}^{} \overline{y} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| y_{i} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}} - {\| \overline{y} \|}^{2}$	`vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/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 Kreuzleistungsdichtespektrum über Slotkorrelation (ohne koinzidente Produkte)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} x_{i}^{} y_{j}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j}}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K)+0j R0=zeros(K) i=0 jb=i while i<N: j=jb while (jb<N) and (t[jb]<t[i]+((K-1)//2+0.5)dt): jb+=1 j=jb-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/float(dt))) R1[k]+=g[i]g[j]conj(x[i])y[j] R0[k]+=g[i]g[j] j-=1 i+=1 Ryx=R1/R0 plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p1*2)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K)+0j; R0=zeros(1,K); i=2; jb=i; while i<=N j=jb; while (jb<=N) && (t(jb)<t(i)+(fix((K-1)/2)+0.5)dt) jb=jb+1; end j=jb-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/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); j=j-1; end i=i+1; end Ryx=R1./R0; plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Slotkorrelation (ohne koinzidente Produkte, mit lokaler Normierung)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{[\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {(x_{i} - \overline{x})}^{} (y_{j} - \overline{y})] \cdot \sqrt{s_{x}^{2} s_{y}^{2}}}{\sqrt{[\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {\| x_{i} - \overline{x} \|}^{2}] [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {\| y_{j} - \overline{y} \|}^{2}]}} + {\overline{x}}^{*} \overline{y}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/float(sum(g)))*2 R1=zeros(K)+0j R2=zeros(K) R3=zeros(K) i=0 jb=i while i<N: j=jb while (jb<N) and (t[jb]<t[i]+((K-1)//2+0.5)dt): jb+=1 j=jb-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/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 j-=1 i+=1 Ryx=sqrt(vxevye)R1/sqrt(R2R3)+conj(mxe)mye plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/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(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); i=2; jb=i; while i<=N j=jb; while (jb<=N) && (t(jb)<t(i)+(fix((K-1)/2)+0.5)dt) jb=jb+1; end j=jb-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/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; j=j-1; end i=i+1; end Ryx=sqrt(vxevye)R1./sqrt(R2.R3)+conj(mxe)mye; plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der koinzidenten Produkte, ohne Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{} y_{i}}{{(\sum_{i = 0}^{N - 1} g_{i})}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{*} (- τ_{k}) = \frac{R_{E, y x, k}}{T - \| τ_{k} \|}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(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(-2jpifpt[i]) Y+=g[i]y[i]exp(-2jpifpt[i]) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/float(sum(g)2-sum(g2)) REyx=ifft(Eyx)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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-abs(tau[k])) plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[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(-2jpifpt(i)); Y=Y+g(i)y(i)exp(-2jpifpt(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=200; K=min(ceil(2T/dt)-1,K); 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-abs(tau(mod(k+K,K)+1))); end plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der koinzidenten Produkte, mit Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{} y_{i}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = T^{2} \frac{{\| X_{j}^{'} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k}^{'} = R_{E, y x}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{*} (- τ_{k}) = \frac{R_{E, y x, k}}{R_{E, y x, k}^{'}}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(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(-2jpifpt[i]) Y+=g[i]y[i]exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p1*2)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[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(-2jpifpt(i)); Y=Y+g(i)y(i)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(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=200; K=min(ceil(2T/dt)-1,K); 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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der koinzidenten Produkte, mit lokaler Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = T^{2} \frac{X_{j}^{'' } X_{j}^{'} - \sum_{i = 0}^{N - 1} g_{i}^{2} {\| x_{i} - \overline{x} \|}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{y x, j}^{''} = E_{y x}^{''} (f_{j}) = T^{2} \frac{X_{j}^{' } Y_{j}^{''} - \sum_{i = 0}^{N - 1} g_{i}^{2} {\| y_{i} - \overline{y} \|}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} (x_{i} - \overline{x}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} (y_{i} - \overline{y}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{''} = \sum_{i = 0}^{N - 1} g_{i} {\| x_{i} - \overline{x} \|}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j}^{''} = \sum_{i = 0}^{N - 1} g_{i} {\| y_{i} - \overline{y} \|}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k}^{'} = R_{E, y x}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k}^{''} = R_{E}^{''} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}^{''}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j}^{''} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{R_{E, y x, k} \sqrt{s_{x}^{2} s_{y}^{2}}}{\sqrt{R_{E, y x, k}^{'} R_{E, y x, k}^{''}}} + {\overline{x}}^{} \overline{y}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/float(sum(g)))*2 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(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(-2jpifpt[i]) Y+=g[i](y[i]-mye)exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) Xpp+=g[i]abs(x[i]-mxe)*2exp(-2jpifpt[i]) Ypp+=g[i]abs(y[i]-mxe)*2exp(-2jpifpt[i]) Eyx=T2(conj(X)Y-sum((gconj(x-mxe)(y-mye)))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(conj(Xpp)Xp-sum((gabs(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Eyxpp=T2(conj(Xp)Ypp-sum((gabs(y-mye))2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) REyxpp=real(ifft(Eyxpp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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(vxevye)REyx[k]/sqrt(REyxp[k]REyxpp[k])+conj(mxe)mye plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'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=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[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(-2jpifpt(i)); Y=Y+g(i)(y(i)-mye)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(i)); Xpp=Xpp+g(i)abs(x(i)-mxe)^2exp(-2jpifpt(i)); Ypp=Ypp+g(i)abs(y(i)-mye)^2exp(-2jpifpt(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=200; K=min(ceil(2T/dt)-1,K); 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(vxevye)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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
(keine Lomb-Scargle-Methode für Kreuzspektren in Matlab und Python)
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit Abzug der koinzidenten Produkte, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x1=zeros(J)+0j y1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(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(g2conj(x)y))/float(sum(g)2-sum(g2)) REyx=ifft(Eyx)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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-abs(tau[k])) plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); x1=zeros([1,J])+0j; y1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(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=200; K=min(ceil(2T/dt)-1,K); 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-abs(tau(mod(k+K,K)+1))); end plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit Abzug der koinzidenten Produkte, mit Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x0=zeros(J) x1=zeros(J)+0j y1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(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(g2conj(x)y))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p1*2)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); x0=zeros([1,J]); x1=zeros([1,J])+0j; y1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(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=200; K=min(ceil(2T/dt)-1,K); 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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit Abzug der koinzidenten Produkte, mit lokaler Normierung)		from numpy.fft import * N=len(t) dt=1.0 mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/float(sum(g)))*2 J=int(ceil(2T/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]-Tfloor(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=T2(conj(X)Y-sum(g*2conj(x-mxe)(y-mye)))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(conj(Xpp)Xp-sum((gabs(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Eyxpp=T2(conj(Xp)Ypp-sum((gabs(y-mye))2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) REyxpp=real(ifft(Eyxpp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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(vxevye)REyx[k]/sqrt(REyxp[k]REyxpp[k])+conj(mxe)mye plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'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=ceil(2T/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)-Tfix(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=200; K=min(ceil(2T/dt)-1,K); 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(vxevye)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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
(keine Gewichtung für Interpolation)

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von Wertetripeln $(t_{i}, x_{i}, y_{i})$ mit den Messzeitpunkten $t_{i}$ und den komplexen Messwerten $x_{i}$ und $y_{i}$ (zwei korrelierte AR1-Prozess als Beispiel, mit Erwartungswerten verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)		T=10000.0 Ti=10.0 dr=1.0 mxs=3.0+1.5j mys=2.0+0.5j vxs=1.0 vys=2.0 cc=0.5 c1=sqrt(0.5+0.5sqrt(1-cc2/(vxsvys))) c2=cc/sqrt(vxsvys)/(2c1) c11=sqrt(vxs)c1 c12=sqrt(vxs)c2 c21=sqrt(vys)c2 c22=sqrt(vys)c1 t=[] x=[] y=[] te=0.0 ue=normal(0,sqrt(0.5))+1jnormal(0,sqrt(0.5)) ve=normal(0,sqrt(0.5))+1jnormal(0,sqrt(0.5)) while te<T: tp=exponential(1.0/(2dr)) te+=tp phi=exp(-tp/Ti) theta=sqrt(0.5(1-phi*2)) ue=uephi+normal(0,theta)+1jnormal(0,theta) ve=vephi+normal(0,theta)+1jnormal(0,theta) xe=c11ue+c12ve+mxs ye=c21ue+c22*ve+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=10000; Ti=10; dr=1; mxs=3+1.5j; mys=2+0.5j; vxs=1; vys=2; cc=0.5; c1=sqrt(0.5+0.5sqrt(1-cc^2/(vxsvys))); c2=cc/sqrt(vxsvys)/(2c1); c11=sqrt(vxs)c1; c12=sqrt(vxs)c2; c21=sqrt(vys)c2; c22=sqrt(vys)c1; t=[]; x=[]; y=[]; te=0; ue=sqrt(0.5)randn()+1jsqrt(0.5)randn(); ve=sqrt(0.5)randn()+1jsqrt(0.5)randn(); while te<T tp=-log(1-rand())/(2dr); te=te+tp; phi=exp(-tp/Ti); theta=sqrt(0.5(1-phi^2)); ue=uephi+thetarandn()+1jthetarandn(); ve=vephi+thetarandn()+1jthetarandn(); xe=c11ue+c12ve+mxs; ye=c21ue+c22ve+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 $g_{i}$ passend zu den $(t_{i}, x_{i}, y_{i})$ (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	$\overline{x} = \frac{\sum_{i = 0}^{N - 1} g_{i} x_{i}}{\sum_{i = 0}^{N - 1} g_{i}}$ $\overline{y} = \frac{\sum_{i = 0}^{N - 1} g_{i} y_{i}}{\sum_{i = 0}^{N - 1} g_{i}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g))`	`mxe=sum(g.x)/sum(g) mye=sum(g.y)/sum(g)`
Varianzen (ohne Bessel-Korrektur, asymptotisch erwartungstreu)	$s_{x}^{2} = \frac{\sum_{i = 0}^{N - 1} g_{i} {(x_{i} - \overline{x})}^{} (x_{i} - \overline{x})}{\sum_{i = 0}^{N - 1} g_{i}} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| x_{i} - \overline{x} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}}$ $= \frac{\sum_{i = 0}^{N - 1} g_{i} x_{i}^{} x_{i}}{\sum_{i = 0}^{N - 1} g_{i}} - {\overline{x}}^{} \overline{x} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| x_{i} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}} - {\| \overline{x} \|}^{2}$ $s_{y}^{2} = \frac{\sum_{i = 0}^{N - 1} g_{i} {(y_{i} - \overline{y})}^{} (y_{i} - \overline{y})}{\sum_{i = 0}^{N - 1} g_{i}} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| y_{i} - \overline{y} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}}$ $= \frac{\sum_{i = 0}^{N - 1} g_{i} y_{i}^{} y_{i}}{\sum_{i = 0}^{N - 1} g_{i}} - {\overline{y}}^{} \overline{y} = \frac{\sum_{i = 0}^{N - 1} g_{i} {\| y_{i} \|}^{2}}{\sum_{i = 0}^{N - 1} g_{i}} - {\| \overline{y} \|}^{2}$	`vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/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 Kreuzleistungsdichtespektrum über Slotkorrelation (ohne koinzidente Produkte)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} x_{i}^{} y_{j}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j}}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K)+0j R0=zeros(K) i=0 jb=i while i<N: j=jb while (jb<N) and (t[jb]<t[i]+((K-1)//2+0.5)dt): jb+=1 j=jb-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/float(dt))) R1[k]+=g[i]g[j]conj(x[i])y[j] R0[k]+=g[i]g[j] j-=1 i+=1 Ryx=R1/R0 plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p1*2)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K)+0j; R0=zeros(1,K); i=2; jb=i; while i<=N j=jb; while (jb<=N) && (t(jb)<t(i)+(fix((K-1)/2)+0.5)dt) jb=jb+1; end j=jb-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/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); j=j-1; end i=i+1; end Ryx=R1./R0; plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Slotkorrelation (ohne koinzidente Produkte, mit lokaler Normierung)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{[\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {(x_{i} - \overline{x})}^{} (y_{j} - \overline{y})] \cdot \sqrt{s_{x}^{2} s_{y}^{2}}}{\sqrt{[\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {\| x_{i} - \overline{x} \|}^{2}] [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {\| y_{j} - \overline{y} \|}^{2}]}} + {\overline{x}}^{*} \overline{y}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/float(sum(g)))*2 R1=zeros(K)+0j R2=zeros(K) R3=zeros(K) i=0 jb=i while i<N: j=jb while (jb<N) and (t[jb]<t[i]+((K-1)//2+0.5)dt): jb+=1 j=jb-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/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 j-=1 i+=1 Ryx=sqrt(vxevye)R1/sqrt(R2R3)+conj(mxe)mye plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/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(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); i=2; jb=i; while i<=N j=jb; while (jb<=N) && (t(jb)<t(i)+(fix((K-1)/2)+0.5)dt) jb=jb+1; end j=jb-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/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; j=j-1; end i=i+1; end Ryx=sqrt(vxevye)R1./sqrt(R2.R3)+conj(mxe)mye; plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der koinzidenten Produkte, ohne Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{} y_{i}}{{(\sum_{i = 0}^{N - 1} g_{i})}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{*} (- τ_{k}) = \frac{R_{E, y x, k}}{T - \| τ_{k} \|}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(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(-2jpifpt[i]) Y+=g[i]y[i]exp(-2jpifpt[i]) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/float(sum(g)2-sum(g2)) REyx=ifft(Eyx)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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-abs(tau[k])) plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[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(-2jpifpt(i)); Y=Y+g(i)y(i)exp(-2jpifpt(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=200; K=min(ceil(2T/dt)-1,K); 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-abs(tau(mod(k+K,K)+1))); end plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der koinzidenten Produkte, mit Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{} y_{i}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = T^{2} \frac{{\| X_{j}^{'} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k}^{'} = R_{E, y x}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{*} (- τ_{k}) = \frac{R_{E, y x, k}}{R_{E, y x, k}^{'}}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(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(-2jpifpt[i]) Y+=g[i]y[i]exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p1*2)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[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(-2jpifpt(i)); Y=Y+g(i)y(i)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(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=200; K=min(ceil(2T/dt)-1,K); 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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der koinzidenten Produkte, mit lokaler Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = T^{2} \frac{X_{j}^{'' } X_{j}^{'} - \sum_{i = 0}^{N - 1} g_{i}^{2} {\| x_{i} - \overline{x} \|}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{y x, j}^{''} = E_{y x}^{''} (f_{j}) = T^{2} \frac{X_{j}^{' } Y_{j}^{''} - \sum_{i = 0}^{N - 1} g_{i}^{2} {\| y_{i} - \overline{y} \|}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} (x_{i} - \overline{x}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} (y_{i} - \overline{y}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{''} = \sum_{i = 0}^{N - 1} g_{i} {\| x_{i} - \overline{x} \|}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j}^{''} = \sum_{i = 0}^{N - 1} g_{i} {\| y_{i} - \overline{y} \|}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k}^{'} = R_{E, y x}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k}^{''} = R_{E}^{''} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}^{''}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j}^{''} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{R_{E, y x, k} \sqrt{s_{x}^{2} s_{y}^{2}}}{\sqrt{R_{E, y x, k}^{'} R_{E, y x, k}^{''}}} + {\overline{x}}^{} \overline{y}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $S_{y x, j} = S_{y x} (f_{j}) = Δ t \cdot F F T {R_{y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/float(sum(g)))*2 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(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(-2jpifpt[i]) Y+=g[i](y[i]-mye)exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) Xpp+=g[i]abs(x[i]-mxe)*2exp(-2jpifpt[i]) Ypp+=g[i]abs(y[i]-mxe)*2exp(-2jpifpt[i]) Eyx=T2(conj(X)Y-sum((gconj(x-mxe)(y-mye)))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(conj(Xpp)Xp-sum((gabs(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Eyxpp=T2(conj(Xp)Ypp-sum((gabs(y-mye))2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) REyxpp=real(ifft(Eyxpp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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(vxevye)REyx[k]/sqrt(REyxp[k]REyxpp[k])+conj(mxe)mye plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'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=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[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(-2jpifpt(i)); Y=Y+g(i)(y(i)-mye)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(i)); Xpp=Xpp+g(i)abs(x(i)-mxe)^2exp(-2jpifpt(i)); Ypp=Ypp+g(i)abs(y(i)-mye)^2exp(-2jpifpt(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=200; K=min(ceil(2T/dt)-1,K); 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(vxevye)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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
(keine Lomb-Scargle-Methode für Kreuzspektren in Matlab und Python)
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit Abzug der koinzidenten Produkte, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x1=zeros(J)+0j y1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(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(g2conj(x)y))/float(sum(g)2-sum(g2)) REyx=ifft(Eyx)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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-abs(tau[k])) plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); x1=zeros([1,J])+0j; y1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(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=200; K=min(ceil(2T/dt)-1,K); 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-abs(tau(mod(k+K,K)+1))); end plot(tau,real(Ryx),'o',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit Abzug der koinzidenten Produkte, mit Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x0=zeros(J) x1=zeros(J)+0j y1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(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(g2conj(x)y))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p1*2)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); x0=zeros([1,J]); x1=zeros([1,J])+0j; y1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(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=200; K=min(ceil(2T/dt)-1,K); 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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)),f,imag(Syx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsdichtespektrum über Zeitquantisierung (mit Abzug der koinzidenten Produkte, mit lokaler Normierung)		from numpy.fft import * N=len(t) dt=1.0 mxe=sum(gx)/float(sum(g)) mye=sum(gy)/float(sum(g)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(gabs(y)*2)/float(sum(g))-abs(sum(gy)/float(sum(g)))*2 J=int(ceil(2T/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]-Tfloor(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=T2(conj(X)Y-sum(g*2conj(x-mxe)(y-mye)))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(conj(Xpp)Xp-sum((gabs(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Eyxpp=T2(conj(Xp)Ypp-sum((gabs(y-mye))2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) REyxpp=real(ifft(Eyxpp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) 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(vxevye)REyx[k]/sqrt(REyxp[k]REyxpp[k])+conj(mxe)mye plot(tau,real(Ryx),'o',arange(-5K,5K)dt/10.0,ccsqrt(vxsvys)exp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+real(conj(mxs)mys),tau,imag(Ryx),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Syx=dtfft(Ryx) p1=1-dt/float(Ti) loglog(f,real(Syx),'o',arange(1,5K)/float(10Kdt),ccsqrt(vxsvys)(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt)),f,imag(Syx),'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=ceil(2T/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)-Tfix(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=200; K=min(ceil(2T/dt)-1,K); 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(vxevye)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',(-5K:5K)dt/10,ccsqrt(vxsvys)exp(-abs((-5K:5K)dt/10/Ti))+real(conj(mxs)mys),tau,imag(Ryx),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Syx=dtfft(Ryx); p1=1-dt/Ti; loglog(f,real(Syx),'o',(1:5K)/(10Kdt),ccsqrt(vxsvys)(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)),f,imag(Syx),'o')
(keine Gewichtung für Interpolation)