Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

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

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von Wertepaaren $(t_{x, i}, x_{i})$ mit den Messzeitpunkten $t_{x, i}$ und den Messwerten $x_{i}$ mit $i = 0 \dots N_{x} - 1$ sowie den Wertepaaren $(t_{y, i}, y_{i})$ mit den Messzeitpunkten $t_{y, i}$ und den Messwerten $y_{i}$ mit $i = 0 \dots N_{y} - 1$ (zwei harmonische Signale als Beispiel, mit einem Erwartungswert verschieden von null)		`T=1000.0 drx=4.0 dry=6.0 mxs=3.0 mys=2.0 axs=1.5 ays=2.0 tx=[] ty=[] x=[] y=[] te=0.0 while te<T: tp=exponential(1.0/(drx+dry)) te+=tp if te<T: if random()<drx/float(drx+dry): tx.append(te) x.append(axssin(0.05pite)+mxs) else: ty.append(te) y.append(ayssin(0.05pite+0.5)+mys) tx=array(tx) ty=array(ty) x=array(x) y=array(y) plot(tx,x,'o',ty,y,'o') show()`	`T=1000; drx=4; dry=6; mxs=3; mys=2; axs=1.5; ays=2; tx=[]; ty=[]; x=[]; y=[]; te=0; while te<T tp=-log(1-rand())/(drx+dry); te=te+tp; if te<T if rand<drx/(drx+dry) tx(end+1)=te; x(end+1)=axssin(0.05pite)+mxs; else ty(end+1)=te; y(end+1)=ayssin(0.05pite+0.5)+mys; end end end plot(tx,x,'o',ty,y,'o')`
Mittelwerte	$\overline{x} = \frac{1}{N_{x}} \cdot \sum_{i = 0}^{N_{x} - 1} x_{i}$ $\overline{y} = \frac{1}{N_{y}} \cdot \sum_{i = 0}^{N_{y} - 1} y_{i}$	`mxe=mean(x) mye=mean(y)`	`mxe=mean(x) mye=mean(y)`
Varianzen (ohne Bessel-Korrektur, asymptotisch erwartungstreu)	$s_{x}^{2} = \frac{1}{N_{x}} \cdot \sum_{i = 0}^{N_{x} - 1} {(x_{i} - \overline{x})}^{2} = (\frac{1}{N_{x}} \cdot \sum_{i = 0}^{N_{x} - 1} x_{i}^{2}) - {\overline{x}}^{2}$ $s_{y}^{2} = \frac{1}{N_{y}} \cdot \sum_{i = 0}^{N_{y} - 1} {(y_{i} - \overline{y})}^{2} = (\frac{1}{N_{y}} \cdot \sum_{i = 0}^{N_{y} - 1} y_{i}^{2}) - {\overline{y}}^{2}$	`vxe=var(x) vye=var(y)`	`vxe=var(x,1) vye=var(y,1)`
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Slotkorrelation (auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y} (- τ_{k}) = \frac{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1} b_{k} (t_{j} - t_{i}) x_{i} y_{j}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1} b_{k} (t_{j} - t_{i})}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t + c T \| < Δ t / 2 für irgendein ganzes c \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [T / Δ t]$ $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{K} \cdot F F T {R_{y x, k}} = \frac{1}{K} \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 * Nx=len(tx) Ny=len(ty) dt=1.0 K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K) R0=zeros(K) for i in range(0,Nx): for j in range(0,Ny): k=int(round((ty[j]-tx[i]-Tround((ty[j]-tx[i])/float(T)))/float(dt))) R1[k]+=x[i]y[j] R0[k]+=1.0 Ryx=R1/R0 plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()`	`Nx=length(tx); Ny=length(ty); dt=1.0; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K); R0=zeros(1,K); for i=1:Nx for j=1:Ny k=round((ty(j)-tx(i)-Tround((ty(j)-tx(i))/T))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+x(i)y(j); R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+1; end end Ryx=R1./R0; plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')`
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Slotkorrelation (mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y} (- τ_{k}) = \frac{[\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1} b_{k} (t_{j} - t_{i}) (x_{i} - \overline{x}) (y_{j} - \overline{y})] \cdot \sqrt{s_{x}^{2} s_{y}^{2}}}{\sqrt{[\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1} b_{k} (t_{j} - t_{i}) {(x_{i} - \overline{x})}^{2}] [\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1} b_{k} (t_{j} - t_{i}) {(y_{j} - \overline{y})}^{2}]}} + \overline{x} \cdot \overline{y}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t + c T \| < Δ t / 2 für irgendein ganzes c \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [T / Δ t]$ $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{K} \cdot F F T {R_{y x, k}} = \frac{1}{K} \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 * Nx=len(tx) Ny=len(ty) dt=1.0 K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) mxe=mean(x) mye=mean(y) vxe=var(x) vye=var(y) R1=zeros(K) R2=zeros(K) R3=zeros(K) for i in range(0,Nx): for j in range(0,Ny): k=int(round((ty[j]-tx[i]-Tround((ty[j]-tx[i])/float(T)))/float(dt))) R1[k]+=(x[i]-mxe)(y[j]-mye) R2[k]+=(x[i]-mxe)2 R3[k]+=(y[j]-mye)2 Ryx=sqrt(vxevye)R1/sqrt(R2R3)+mxemye plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); mxe=mean(x); mye=mean(y); vxe=var(x,1); vye=var(y,1); R1=zeros(1,K); R2=zeros(1,K); R3=zeros(1,K); for i=1:Nx for j=1:Ny k=round((ty(j)-tx(i)-Tround((ty(j)-tx(i))/T))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+(x(i)-mxe)(y(j)-mye); R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+(x(i)-mxe)^2; R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+(y(j)-mye)^2; end end Ryx=sqrt(vxevye)R1./sqrt(R2.R3)+mxemye; plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, ohne Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = \frac{T^{2} X_{j}^{*} Y_{j}}{N_{x} N_{y}}$ $X_{j} = \sum_{i = 0}^{N_{x} - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j} = \sum_{i = 0}^{N_{y} - 1} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J = [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{1}{T} R_{E, y x, k}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = J = [T / Δ t]$ $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{K} \cdot F F T {R_{y x, k}} = \frac{1}{K} \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 * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(round(T/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,Nx): X+=x[i]exp(-2jpifptx[i]) for i in range(0,Ny): Y+=y[i]exp(-2jpifpty[i]) Eyx=T2conj(X)Y/float(NxNy) REyx=real(ifft(Eyx))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) Ryx=zeros(K) for k in range(-(K//2),(K+1)//2): Ryx[k]=REyx[k]/float(T) plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=round(T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Y=zeros(size(fp))+0j; for i=1:Nx X=X+x(i)exp(-2jpifptx(i)); end for i=1:Ny Y=Y+y(i)exp(-2jpifpty(i)); end Eyx=T^2conj(X).Y/(NxNy); REyx=real(ifft(Eyx))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); Ryx=zeros([1,K]); for k=-fix(K/2):fix((K-1)/2) Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/T; end plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$E_{y x, j} = E_{y x} (f_{j}) = \frac{T^{2} X_{j}^{} Y_{j}}{N_{x} N_{y}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = \frac{T^{2} X_{j}^{' } Y_{j}^{'}}{N_{x} N_{y}}$ $X_{j} = \sum_{i = 0}^{N_{x} - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j} = \sum_{i = 0}^{N_{y} - 1} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ $X_{j}^{'} = \sum_{i = 0}^{N_{x} - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j}^{'} = \sum_{i = 0}^{N_{y} - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J = [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 = J = [T / Δ t]$ $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{K} \cdot F F T {R_{y x, k}} = \frac{1}{K} \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 * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(round(T/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 Yp=zeros(size(fp))+0j for i in range(0,Nx): X+=x[i]exp(-2jpifptx[i]) Xp+=exp(-2jpifptx[i]) for i in range(0,Ny): Y+=y[i]exp(-2jpifpty[i]) Yp+=exp(-2jpifpty[i]) Eyx=T2conj(X)Y/float(NxNy) Eyxp=T*2conj(Xp)Yp/float(NxNy) REyx=real(ifft(Eyx))/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) Ryx=zeros(K) for k in range(-(K//2),(K+1)//2): Ryx[k]=REyx[k]/float(REyxp[k]) plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=round(T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Y=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; Yp=zeros(size(fp))+0j; for i=1:Nx X=X+x(i)exp(-2jpifptx(i)); Xp=Xp+exp(-2jpifptx(i)); end for i=1:Ny Y=Y+y(i)exp(-2jpifpty(i)); Yp=Yp+exp(-2jpifpty(i)); end Eyx=T^2conj(X).Y/(NxNy); Eyxp=T^2conj(Xp).Yp/(NxNy); REyx=real(ifft(Eyx))/dt; REyxp=real(ifft(Eyxp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); Ryx=zeros([1,K]); for k=-fix(K/2):fix((K-1)/2) Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/REyxp(mod(k+length(REyxp),length(REyxp))+1); end plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$E_{y x, j} = E_{y x} (f_{j}) = \frac{T^{2} X_{j}^{} Y_{j}}{N_{x} N_{y}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = \frac{T^{2} X_{j}^{'' } Y_{j}^{'}}{N_{x} N_{y}}$ $E_{y x, j}^{''} = E_{y x}^{''} (f_{j}) = \frac{T^{2} X_{j}^{' *} Y_{j}^{''}}{N_{x} N_{y}}$ $X_{j} = \sum_{i = 0}^{N_{x} - 1} (x_{i} - \overline{x}) 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j} = \sum_{i = 0}^{N_{y} - 1} (y_{i} - \overline{y}) 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ $X_{j}^{'} = \sum_{i = 0}^{N_{x} - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j}^{'} = \sum_{i = 0}^{N_{y} - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ $X_{j}^{''} = \sum_{i = 0}^{N_{x} - 1} {(x_{i} - \overline{x})}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j}^{''} = \sum_{i = 0}^{N_{y} - 1} {(y_{i} - \overline{y})}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J = [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} \cdot \overline{y}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = J = [T / Δ t]$ $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{K} \cdot F F T {R_{y x, k}} = \frac{1}{K} \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 * Nx=len(tx) Ny=len(ty) dt=1.0 mxe=mean(x) mye=mean(y) vxe=var(x) vye=var(y) J=int(round(T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j Y=zeros(size(fp))+0j Xp=zeros(size(fp))+0j Yp=zeros(size(fp))+0j Xpp=zeros(size(fp))+0j Ypp=zeros(size(fp))+0j for i in range(0,Nx): X+=(x[i]-mxe)exp(-2jpifptx[i]) Xp+=exp(-2jpifptx[i]) Xpp+=(x[i]-mxe)*2exp(-2jpifptx[i]) for i in range(0,Ny): Y+=(y[i]-mye)exp(-2jpifpty[i]) Yp+=exp(-2jpifpty[i]) Ypp+=(y[i]-mye)*2exp(-2jpifpty[i]) Eyx=T2conj(X)Y/float(NxNy) Eyxp=T*2conj(Xpp)Yp/float(NxNy) Eyxpp=T*2conj(Xp)Ypp/float(NxNy) REyx=real(ifft(Eyx))/float(dt) REyxp=real(ifft(Eyxp))/float(dt) REyxpp=real(ifft(Eyxpp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) Ryx=zeros(K) for k in range(-(K//2),(K+1)//2): Ryx[k]=sqrt(vxevye)REyx[k]/sqrt(REyxp[k]REyxpp[k])+mxemye plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; mxe=mean(x); mye=mean(y); vxe=var(x,1); vye=var(y,1); J=round(T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Y=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; Yp=zeros(size(fp))+0j; Xpp=zeros(size(fp))+0j; Ypp=zeros(size(fp))+0j; for i=1:Nx X=X+(x(i)-mxe)exp(-2jpifptx(i)); Xp=Xp+exp(-2jpifptx(i)); Xpp=Xpp+(x(i)-mxe)^2exp(-2jpifptx(i)); end for i=1:Ny Y=Y+(y(i)-mye)exp(-2jpifpty(i)); Yp=Yp+exp(-2jpifpty(i)); Ypp=Ypp+(y(i)-mye)^2exp(-2jpifpty(i)); end Eyx=T^2conj(X).Y/(NxNy); Eyxp=T^2conj(Xpp).Yp/(NxNy); Eyxpp=T^2conj(Xp).Ypp/(NxNy); REyx=real(ifft(Eyx))/dt; REyxp=real(ifft(Eyxp))/dt; REyxpp=real(ifft(Eyxpp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); Ryx=zeros([1,K]); for k=-fix(K/2):fix((K-1)/2) Ryx(mod(k+K,K)+1)=sqrt(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))+mxemye; end plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')
(keine Lomb-Scargle-Methode für Kreuzspektren in Matlab und Python)
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Zeitquantisierung (ohne Normierung)		from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(round(T/float(dt))) x1=zeros(J) y1=zeros(J) for i in range(0,Nx): j=int(floor((tx[i]-Tfloor(tx[i]/float(T)))/float(dt))) x1[j]+=x[i]; for i in range(0,Ny): j=int(floor((ty[i]-Tfloor(ty[i]/float(T)))/float(dt))) y1[j]+=y[i]; X=fft(x1) Y=fft(y1) Eyx=T*2conj(X)Y/float(NxNy) REyx=real(ifft(Eyx))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) Ryx=zeros(K) for k in range(-(K//2),(K+1)//2): Ryx[k]=REyx[k]/float(T) plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=round(T/dt); x1=zeros([1,J]); y1=zeros([1,J]); for i=1:Nx j=fix((tx(i)-Tfix(tx(i)/T))/dt)+1; x1(j)=x1(j)+x(i); end for i=1:Ny j=fix((ty(i)-Tfix(ty(i)/T))/dt)+1; y1(j)=y1(j)+y(i); end X=fft(x1); Y=fft(y1); Eyx=T^2conj(X).Y/(NxNy); REyx=real(ifft(Eyx))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); Ryx=zeros([1,K]); for k=-fix(K/2):fix((K-1)/2) Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/T; end plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Zeitquantisierung (mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)		from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(round(T/float(dt))) x0=zeros(J) y0=zeros(J) x1=zeros(J) y1=zeros(J) for i in range(0,Nx): j=int(floor((tx[i]-Tfloor(tx[i]/float(T)))/float(dt))) x0[j]+=1.0; x1[j]+=x[i]; for i in range(0,Ny): j=int(floor((ty[i]-Tfloor(ty[i]/float(T)))/float(dt))) y0[j]+=1.0; y1[j]+=y[i]; X=fft(x1) Y=fft(y1) Xp=fft(x0) Yp=fft(y0) Eyx=T*2conj(X)Y/float(NxNy) Eyxp=T*2conj(Xp)Yp/float(NxNy) REyx=real(ifft(Eyx))/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) Ryx=zeros(K) for k in range(-(K//2),(K+1)//2): Ryx[k]=REyx[k]/float(REyxp[k]) plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=round(T/dt); x0=zeros([1,J]); y0=zeros([1,J]); x1=zeros([1,J]); y1=zeros([1,J]); for i=1:Nx j=fix((tx(i)-Tfix(tx(i)/T))/dt)+1; x0(j)=x0(j)+1; x1(j)=x1(j)+x(i); end for i=1:Ny j=fix((ty(i)-Tfix(ty(i)/T))/dt)+1; y0(j)=y0(j)+1; y1(j)=y1(j)+y(i); end X=fft(x1); Y=fft(y1); Xp=fft(x0); Yp=fft(y0); Eyx=T^2conj(X).Y/(NxNy); Eyxp=T^2conj(Xp).Yp/(NxNy); REyx=real(ifft(Eyx))/dt; REyxp=real(ifft(Eyxp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); Ryx=zeros([1,K]); for k=-fix(K/2):fix((K-1)/2) Ryx(mod(k+K,K)+1)=REyx(mod(k+length(REyx),length(REyx))+1)/REyxp(mod(k+length(REyxp),length(REyxp))+1); end plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Zeitquantisierung (mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)		from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 mxe=mean(x) mye=mean(y) vxe=var(x) vye=var(y) J=int(round(T/float(dt))) x0=zeros(J) y0=zeros(J) x1=zeros(J) y1=zeros(J) x2=zeros(J) y2=zeros(J) for i in range(0,Nx): j=int(floor((tx[i]-Tfloor(tx[i]/float(T)))/float(dt))) x0[j]+=1.0; x1[j]+=x[i]-mxe; x2[j]+=(x[i]-mxe)2; for i in range(0,Ny): j=int(floor((ty[i]-Tfloor(ty[i]/float(T)))/float(dt))) y0[j]+=1.0; y1[j]+=y[i]-mye; y2[j]+=(y[i]-mye)2; X=fft(x1) Y=fft(y1) Xp=fft(x0) Yp=fft(y0) Xpp=fft(x2) Ypp=fft(y2) Eyx=T2conj(X)Y/float(NxNy) Eyxp=T2conj(Xpp)Yp/float(NxNy) Eyxpp=T*2conj(Xp)Ypp/float(NxNy) REyx=real(ifft(Eyx))/float(dt) REyxp=real(ifft(Eyxp))/float(dt) REyxpp=real(ifft(Eyxpp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) Ryx=zeros(K) for k in range(-(K//2),(K+1)//2): Ryx[k]=sqrt(vxevye)REyx[k]/sqrt(REyxp[k]REyxpp[k])+mxemye plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; mxe=mean(x); mye=mean(y); vxe=var(x,1); vye=var(y,1); J=round(T/dt); x0=zeros([1,J]); y0=zeros([1,J]); x1=zeros([1,J]); y1=zeros([1,J]); x2=zeros([1,J]); y2=zeros([1,J]); for i=1:Nx j=fix((tx(i)-Tfix(tx(i)/T))/dt)+1; x0(j)=x0(j)+1; x1(j)=x1(j)+x(i)-mxe; x2(j)=x2(j)+(x(i)-mxe)^2; end for i=1:Ny j=fix((ty(i)-Tfix(ty(i)/T))/dt)+1; y0(j)=y0(j)+1; y1(j)=y1(j)+y(i)-mye; y2(j)=y2(j)+(y(i)-mye)^2; end X=fft(x1); Y=fft(y1); Xp=fft(x0); Yp=fft(y0); Xpp=fft(x2); Ypp=fft(y2); Eyx=T^2conj(X).Y/(NxNy); Eyxp=T^2conj(Xpp).Yp/(NxNy); Eyxpp=T^2conj(Xp).Ypp/(NxNy); REyx=real(ifft(Eyx))/dt; REyxp=real(ifft(Eyxp))/dt; REyxpp=real(ifft(Eyxpp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); Ryx=zeros([1,K]); for k=-fix(K/2):fix((K-1)/2) Ryx(mod(k+K,K)+1)=sqrt(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))+mxemye; end plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Kreuzkorrelationsfunktion und Kreuzleistungsspektrum über Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)		from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(round(T/float(dt))) xp=x[Nx-1]ones(J) yp=y[Ny-1]ones(J) for i in range(0,Nx-1): for j in range(int(floor(tx[i]/float(dt)))+1,int(floor(tx[i+1]/float(dt)))+1): if j<J: xp[j]=x[i] for i in range(0,Ny-1): for j in range(int(floor(ty[i]/float(dt)))+1,int(floor(ty[i+1]/float(dt)))+1): if j<J: yp[j]=y[i] Xp=fft(xp) Yp=fft(yp) Eyxp=dt*2conj(Xp)Yp 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) cx=exp(-ddx/float(dt))/((1-exp(-ddx/float(dt)))(1-exp(-ddy/float(dt)))) cy=exp(-ddy/float(dt))/((1-exp(-ddx/float(dt)))(1-exp(-ddy/float(dt)))) for k in range(-(K//2),(K+1)//2): Ryx[k]=((cx+cy+1)REyxp[k]-cyREyxp[k-1]-cxREyxp[k+1])/float(T) plot(tau,Ryx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Pyx=fft(Ryx)/float(K) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=round(T/dt); xp=x(Nx)ones([1,J]); yp=y(Ny)ones([1,J]); for i=1:Nx-1 for j=fix(tx(i)/dt)+2:fix(tx(i+1)/dt)+1 if j<=J xp(j)=x(i); end end end for i=1:Ny-1 for j=fix(ty(i)/dt)+2:fix(ty(i+1)/dt)+1 if j<=J yp(j)=y(i); end end end Xp=fft(xp); Yp=fft(yp); Eyxp=dt^2conj(Xp).Yp; 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]); cx=exp(-ddx/dt)/((1-exp(-ddx/dt))(1-exp(-ddy/dt))); cy=exp(-ddy/dt)/((1-exp(-ddx/dt))(1-exp(-ddy/dt))); for k=-fix(K/2):fix((K-1)/2) Ryx(mod(k+K,K)+1)=((cx+cy+1)REyxp(mod(k+length(REyxp),length(REyxp))+1)-cyREyxp(mod(k-1+length(REyxp),length(REyxp))+1)-cxREyxp(mod(k+1+length(REyxp),length(REyxp))+1))/T; end plot(tau,Ryx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); Pyx=fft(Ryx)/K; plot(f,real(Pyx),'o',f,imag(Pyx),'o')