Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für 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 $N$ komplexen Wertepaaren $(x_{i}, y_{i})$ mit $i = 0 \dots N - 1$ (harmonisches Signal als Beispiel, mit einem Erwartungswert verschieden von null, ohne Rauschen)		`N=100 dt=1.0 t=arange(0,N)dt x=3+1j+1.5exp(0.1jpit) y=2+1.5j+1exp(0.1jpi*t-0.3) plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o') show()`	`N=100; dt=1.0; t=(0:N-1)dt; x=3+1j+1.5exp(0.1jpit); y=2+1.5j+1exp(0.1jpi*t-0.3); plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')`
Überlagerung von Rauschen		`x+=normal(0,0.1,N)+1jnormal(0,0.1,N) y+=normal(0,0.1,N)+1jnormal(0,0.1,N) plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o') show()`	`x=x+0.1randn(1,N)+1j0.1randn(1,N); y=y+0.1randn(1,N)+1j0.1randn(1,N); plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')`
Generieren von $N$ reellen Gewichten $(g_{i}, h_{i})$ mit $i = 0 \dots N - 1$ , passend zu den $(x_{i}, y_{i})$ mit $i = 0 \dots N - 1$ (Kosinusfenster als Beispiel)		`g=sin(pi(arange(0,N)+0.5)/float(N)) h=sin(pi(arange(0,N)+0.5)/float(N)) plot(t,g,'o',t,h,'o') show()`	`g=sin(pi((0:N-1)+0.5)/N); h=sin(pi((0:N-1)+0.5)/N); plot(t,g,'o',t,h,'o')`
Mittelwerte (erwartungstreu)	$\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} h_{i} y_{i}}{\sum_{i = 0}^{N - 1} h_{i}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h);`
Varianzen (ohne Rauschen asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist, mit Rauschen systematische Fehler)	$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} h_{i} {(y_{i} - \overline{y})}^{} (y_{i} - \overline{y})}{\sum_{i = 0}^{N - 1} h_{i}} = \frac{\sum_{i = 0}^{N - 1} h_{i} {\| y_{i} - \overline{y} \|}^{2}}{\sum_{i = 0}^{N - 1} h_{i}}$ $= \frac{\sum_{i = 0}^{N - 1} h_{i} y_{i}^{} y_{i}}{\sum_{i = 0}^{N - 1} h_{i}} - {\overline{y}}^{} \overline{y} = \frac{\sum_{i = 0}^{N - 1} h_{i} {\| y_{i} \|}^{2}}{\sum_{i = 0}^{N - 1} h_{i}} - {\| \overline{y} \|}^{2}$	`vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))**2`	`vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2;`
Kreuzkovarianz (asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist)	$c_{y x} = c_{x y}^{} = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{i} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})}{\sum_{i = 0}^{N - 1} g_{i} h_{i}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) sum(ghconj(x-mxe)(y-mye))/float(sum(gh))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.conj(x-mxe).(y-mye))/sum(g.h)`
(Kreuz-)Korrelationskoeffizient (ohne Rauschen asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist, mit Rauschen systematische Fehler)	$γ_{y x} = γ_{x y}^{} = \frac{c_{y x}}{\sqrt{s_{x}^{2} s_{y}^{2}}} = \frac{[\sum_{i = 0}^{N - 1} g_{i} h_{i} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})] \cdot \sqrt{(\sum_{i = 0}^{N - 1} g_{i}) (\sum_{i = 0}^{N - 1} h_{i})}}{(\sum_{i = 0}^{N - 1} g_{i} h_{i}) \cdot \sqrt{(\sum_{i = 0}^{N - 1} g_{i} {\| x_{i} - \overline{x} \|}^{2}) (\sum_{i = 0}^{N - 1} h_{i} {\| y_{i} - \overline{y} \|}^{2})}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) sum(ghconj(x-mxe)(y-mye))sqrt(sum(g)sum(h))/(sum(gh)sqrt(sum(gabs(x-mxe)*2)sum(habs(y-mye)*2)))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.conj(x-mxe).(y-mye))sqrt(sum(g)sum(h))/(sum(g.h)sqrt(sum(g.abs(x-mxe).^2)sum(h.abs(y-mye).^2)))`
(Kreuz-)Korrelationskoeffizient (mit lokaler Normierung, ohne Rauschen asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist, mit Rauschen systematische Fehler)	$γ_{y x} = γ_{x y}^{} = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{i} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})}{\sqrt{(\sum_{i = 0}^{N - 1} g_{i} h_{i} {\| x_{i} - \overline{x} \|}^{2}) (\sum_{i = 0}^{N - 1} g_{i} h_{i} {\| y_{i} - \overline{y} \|}^{2})}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) sum(ghconj(x-mxe)(y-mye))/sqrt(sum(ghabs(x-mxe)2)sum(ghabs(y-mye)**2))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.conj(x-mxe).(y-mye))/sqrt(sum(g.h.abs(x-mxe).^2)sum(g.h.abs(y-mye).^2))`
Kreuzkorrelationsfunktion (erwartungstreu)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N} x_{i}^{} y_{(i + k) mod N}}{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `Ryx=zeros(N)+0j tau=arange(0,N)dt for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i])y[(i+k)%N] sum0+=g[i]h[(i+k)%N] Ryx[k]=sum1/float(sum0) plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=arange(0,N)dt X1=fft(gx) Y1=fft(hy) P1=conj(X1)Y1/N*2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N**2 Ryx=ifft(P1)/ifft(P0) plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()`	direkte Berechnung: `Ryx=zeros(1,N)+0j; tau=(0:N-1)dt; for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i))y(mod(i+k-1,N)+1); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end Ryx(k+1)=sum1/sum0; end plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; X1=fft(g.x); Y1=fft(h.y); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; Ryx=ifft(P1)./ifft(P0); plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{\underset{j - i \equiv k mod {\tilde{N}}_{x}}{\sum_{i = 0}^{N_{y} - 1} \sum_{j = 0}^{N - 1}} g_{i} h_{j} x_{i}^{} y_{j}}{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) R1=zeros(NT)+0j R0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): R1[(j-i)%NT]+=g[i]h[j]conj(x[i])y[j] R0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt Ryx=zeros(NT)+0j for k in range(0,NT): Ryx[k]=R1[k]/R0[k] plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i]x[i] wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i]y[i] wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) tau=arange(0,NT)dt P1=conj(X1)Y1/NT*2 P0=conj(X0)Y0/NT**2 Ryx=ifft(P1)/real(ifft(P0)) plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()`	direkte Berechnung: `NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); R1=zeros(1,NT)+0j; R0=zeros(1,NT); for i=1:Nx for j=1:Ny R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i))y(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+g(i)h(j); end end tau=(0:NT-1)dt; Ryx=zeros(1,NT)+0j; for k=1:NT Ryx(k)=R1(k)/R0(k); end plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)x(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)y(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); tau=(0:NT-1)dt; P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; Ryx=ifft(P1)./real(ifft(P0)); plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Kreuzkovarianzfunktion (asymptotisch erwartungstreu)	$C_{y x, k} = C_{y x} (τ_{k}) = C_{x y}^{} (- τ_{k}) = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N} {(x_{i} - \overline{x})}^{} [y_{(i + k) mod N} - \overline{y}]}{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `Cyx=zeros(N)+0j tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i]-mxe)(y[(i+k)%N]-mye) sum0+=g[i]h[(i+k)%N] Cyx[k]=sum1/float(sum0) plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) X1=fft(g(x-mxe)) Y1=fft(h(y-mye)) P1=conj(X1)Y1/N*2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N**2 Cyx=ifft(P1)/ifft(P0) plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()`	direkte Berechnung: `Cyx=zeros(1,N)+0j; tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i)-mxe)(y(mod(i+k-1,N)+1)-mye); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end Cyx(k+1)=sum1/sum0; end plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); X1=fft(g.(x-mxe)); Y1=fft(h.(y-mye)); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; Cyx=ifft(P1)./ifft(P0); plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	$C_{y x, k} = C_{y x} (τ_{k}) = C_{x y}^{} (- τ_{k}) = \frac{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j} {(x_{i} - \overline{x})}^{} (y_{j} - \overline{y})}{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) C1=zeros(NT)+0j C0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): C1[(j-i)%NT]+=g[i]h[j]conj(x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt Cyx=zeros(NT)+0j for k in range(0,NT): Cyx[k]=C1[k]/C0[k] plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: from numpy.fft import mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i](x[i]-mxe) wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i](y[i]-mye) wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) tau=arange(0,NT)dt P1=conj(X1)Y1/NT*2 P0=conj(X0)Y0/NT**2 Cyx=ifft(P1)/real(ifft(P0)) plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()	direkte Berechnung: `mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); C1=zeros(1,NT)+0j; C0=zeros(1,NT); for i=1:Nx for j=1:Ny C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i)-mxe)(y(j)-mye); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+g(i)h(j); end end tau=(0:NT-1)dt; Cyx=zeros(1,NT)+0j; for k=1:NT Cyx(k)=C1(k)/C0(k); end plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)(x(i)-mxe); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)(y(i)-mye); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); tau=(0:NT-1)dt; P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; Cyx=ifft(P1)./real(ifft(P0)); plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Kreuzkorrelationskoeffizientenfunktion (ohne Rauschen asymptotisch erwartungstreu, mit Rauschen systematische Fehler)	$ρ_{y x, k} = ρ_{y x} (τ_{k}) = ρ_{x y}^{} (- τ_{k}) = \frac{C_{y x, k}}{\sqrt{s_{x}^{2} s_{y}^{2}}} = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N} {(x_{i} - \overline{x})}^{} [y_{(i + k) mod N} - \overline{y}]}{\sqrt{s_{x}^{2} s_{y}^{2}} \cdot \sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `rhoyx=zeros(N)+0j tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)*2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i]-mxe)(y[(i+k)%N]-mye) sum0+=g[i]h[(i+k)%N] rhoyx[k]=sum1/float(sum0)/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)*2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 X1=fft(g(x-mxe)) Y1=fft(h(y-mye)) P1=conj(X1)Y1/N*2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N*2 rhoyx=ifft(P1)/ifft(P0)/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show()`	direkte Berechnung: `rhoyx=zeros(1,N)+0j; tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i)-mxe)(y(mod(i+k-1,N)+1)-mye); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end rhoyx(k+1)=sum1/sum0/sqrt(vxevye); end plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; X1=fft(g.(x-mxe)); Y1=fft(h.(y-mye)); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; rhoyx=ifft(P1)./ifft(P0)/sqrt(vxevye); plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	$ρ_{y x, k} = ρ_{y x} (τ_{k}) = ρ_{x y}^{} (- τ_{k}) = \frac{C_{y x, k}}{\sqrt{s_{x}^{2} s_{y}^{2}}} = \frac{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j} {(x_{i} - \overline{x})}^{} (y_{j} - \overline{y})}{\sqrt{s_{x}^{2} s_{y}^{2}} \cdot \underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) C1=zeros(NT)+0j C0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): C1[(j-i)%NT]+=g[i]h[j]conj(x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt rhoyx=zeros(NT)+0j for k in range(0,NT): rhoyx[k]=C1[k]/C0[k]/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show() Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: from numpy.fft import mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i](x[i]-mxe) wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i](y[i]-mye) wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) tau=arange(0,NT)dt P1=conj(X1)Y1/NT2 P0=conj(X0)Y0/NT*2 rhoyx=ifft(P1)/real(ifft(P0))/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show()	direkte Berechnung: mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); C1=zeros(1,NT)+0j; C0=zeros(1,NT); for i=1:Nx for j=1:Ny C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i)-mxe)(y(j)-mye); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+g(i)h(j); end end tau=(0:NT-1)dt; rhoyx=zeros(1,NT)+0j; for k=1:NT rhoyx(k)=C1(k)/C0(k)/sqrt(vxevye); end plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)(x(i)-mxe); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)(y(i)-mye); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); tau=(0:NT-1)dt; P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; rhoyx=ifft(P1)./real(ifft(P0))/sqrt(vxevye); plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Kreuzleistungsspektrum (erwartungstreu)	(wegen der nötigen Entfaltung mit dem Leistungsspektrum der Gewichte nur über die Korrelationsfunktion) $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{N} \cdot F F T {R_{y x, k}} = \frac{1}{N} \cdot \sum_{k = 0}^{N - 1} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / N}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{N Δ t}$ $j = - ⌊ N / 2 ⌋ \dots ⌊ (N - 1) / 2 ⌋$	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `from numpy.fft import * X1=fft(gx) Y1=fft(hy) P1=conj(X1)Y1/N2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N*2 Ryx=ifft(P1)/ifft(P0) f=roll(arange(-(N//2),N-(N//2))/float(Ndt),N-(N//2)) Pyx=fft(Ryx)/float(N) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import * Ryx=zeros(N)+0j for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i])y[(i+k)%N] sum0+=g[i]h[(i+k)%N] Ryx[k]=sum1/float(sum0) f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2)) Pyx=fft(Ryx)/float(N) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `X1=fft(g.x); Y1=fft(h.y); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; Ryx=ifft(P1)./ifft(P0); f=circshift(((-N/2:N-N/2-1))/(Ndt),[0;N/2]); Pyx=fft(Ryx)/N; plot(f,real(Pyx),'o',f,imag(Pyx),'o')` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `Ryx=zeros(1,N)+0j; for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i))y(mod(i+k-1,N)+1); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end Ryx(k+1)=sum1/sum0; end f=circshift(((-N/2:N-N/2-1))/(Ndt),[0;N/2]); Pyx=fft(Ryx)/N; plot(f,real(Pyx),'o',f,imag(Pyx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	(wegen der nötigen Entfaltung mit dem Leistungsspektrum der Gewichte nur über die Korrelationsfunktion) $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{\tilde{N}} \cdot F F T {R_{y x, k}} = \frac{1}{\tilde{N}} \cdot \sum_{k = 0}^{\tilde{N} - 1} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / \tilde{N}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{\tilde{N} Δ t}$ $j = - ⌊ \tilde{N} / 2 ⌋ \dots ⌊ (\tilde{N} - 1) / 2 ⌋$	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `from numpy.fft import * NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i]x[i] wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i]y[i] wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) P1=conj(X1)Y1/NT2 P0=conj(X0)Y0/NT*2 Ryx=ifft(P1)/real(ifft(P0)) f=roll(arange(-(NT//2),(NT+1)//2)/float(NTdt),(NT+1)//2) Pyx=fft(Ryx)/float(NT) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import * NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) R1=zeros(NT)+0j R0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): R1[(j-i)%NT]+=g[i]h[j]conj(x[i])y[j] R0[(j-i)%NT]+=g[i]h[j] Ryx=zeros(NT)+0j for k in range(0,NT): Ryx[k]=R1[k]/R0[k] f=roll(arange(-(NT//2),(NT+1)//2)/float(NT*dt),(NT+1)//2) Pyx=fft(Ryx)/float(NT) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)x(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)y(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; Ryx=ifft(P1)./real(ifft(P0)); f=circshift((-fix(NT/2):fix((NT-1)/2))/(NTdt),[0;fix((NT+1)/2)]); Pyx=fft(Ryx)/NT; plot(f,real(Pyx),'o',f,imag(Pyx),'o') Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); R1=zeros(1,NT)+0j; R0=zeros(1,NT); for i=1:Nx for j=1:Ny R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i))y(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+g(i)h(j); end end Ryx=zeros(1,NT)+0j; for k=1:NT Ryx(k)=R1(k)/R0(k); end f=circshift((-fix(NT/2):fix((NT-1)/2))/(NTdt),[0;fix((NT+1)/2)]); Pyx=fft(Ryx)/NT; plot(f,real(Pyx),'o',f,imag(Pyx),'o')`

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von $N$ komplexen Wertepaaren $(x_{i}, y_{i})$ mit $i = 0 \dots N - 1$ (harmonisches Signal als Beispiel, mit einem Erwartungswert verschieden von null, ohne Rauschen)		`N=100 dt=1.0 t=arange(0,N)dt x=3+1j+1.5exp(0.1jpit) y=2+1.5j+1exp(0.1jpi*t-0.3) plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o') show()`	`N=100; dt=1.0; t=(0:N-1)dt; x=3+1j+1.5exp(0.1jpit); y=2+1.5j+1exp(0.1jpi*t-0.3); plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')`
Überlagerung von Rauschen		`x+=normal(0,0.1,N)+1jnormal(0,0.1,N) y+=normal(0,0.1,N)+1jnormal(0,0.1,N) plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o') show()`	`x=x+0.1randn(1,N)+1j0.1randn(1,N); y=y+0.1randn(1,N)+1j0.1randn(1,N); plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')`
Generieren von $N$ reellen Gewichten $(g_{i}, h_{i})$ mit $i = 0 \dots N - 1$ , passend zu den $(x_{i}, y_{i})$ mit $i = 0 \dots N - 1$ (Kosinusfenster als Beispiel)		`g=sin(pi(arange(0,N)+0.5)/float(N)) h=sin(pi(arange(0,N)+0.5)/float(N)) plot(t,g,'o',t,h,'o') show()`	`g=sin(pi((0:N-1)+0.5)/N); h=sin(pi((0:N-1)+0.5)/N); plot(t,g,'o',t,h,'o')`
Mittelwerte (erwartungstreu)	$\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} h_{i} y_{i}}{\sum_{i = 0}^{N - 1} h_{i}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h);`
Varianzen (ohne Rauschen asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist, mit Rauschen systematische Fehler)	$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} h_{i} {(y_{i} - \overline{y})}^{} (y_{i} - \overline{y})}{\sum_{i = 0}^{N - 1} h_{i}} = \frac{\sum_{i = 0}^{N - 1} h_{i} {\| y_{i} - \overline{y} \|}^{2}}{\sum_{i = 0}^{N - 1} h_{i}}$ $= \frac{\sum_{i = 0}^{N - 1} h_{i} y_{i}^{} y_{i}}{\sum_{i = 0}^{N - 1} h_{i}} - {\overline{y}}^{} \overline{y} = \frac{\sum_{i = 0}^{N - 1} h_{i} {\| y_{i} \|}^{2}}{\sum_{i = 0}^{N - 1} h_{i}} - {\| \overline{y} \|}^{2}$	`vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))**2`	`vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2;`
Kreuzkovarianz (asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist)	$c_{y x} = c_{x y}^{} = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{i} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})}{\sum_{i = 0}^{N - 1} g_{i} h_{i}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) sum(ghconj(x-mxe)(y-mye))/float(sum(gh))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.conj(x-mxe).(y-mye))/sum(g.h)`
(Kreuz-)Korrelationskoeffizient (ohne Rauschen asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist, mit Rauschen systematische Fehler)	$γ_{y x} = γ_{x y}^{} = \frac{c_{y x}}{\sqrt{s_{x}^{2} s_{y}^{2}}} = \frac{[\sum_{i = 0}^{N - 1} g_{i} h_{i} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})] \cdot \sqrt{(\sum_{i = 0}^{N - 1} g_{i}) (\sum_{i = 0}^{N - 1} h_{i})}}{(\sum_{i = 0}^{N - 1} g_{i} h_{i}) \cdot \sqrt{(\sum_{i = 0}^{N - 1} g_{i} {\| x_{i} - \overline{x} \|}^{2}) (\sum_{i = 0}^{N - 1} h_{i} {\| y_{i} - \overline{y} \|}^{2})}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) sum(ghconj(x-mxe)(y-mye))sqrt(sum(g)sum(h))/(sum(gh)sqrt(sum(gabs(x-mxe)*2)sum(habs(y-mye)*2)))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.conj(x-mxe).(y-mye))sqrt(sum(g)sum(h))/(sum(g.h)sqrt(sum(g.abs(x-mxe).^2)sum(h.abs(y-mye).^2)))`
(Kreuz-)Korrelationskoeffizient (mit lokaler Normierung, ohne Rauschen asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist, mit Rauschen systematische Fehler)	$γ_{y x} = γ_{x y}^{} = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{i} {(x_{i} - \overline{x})}^{} (y_{i} - \overline{y})}{\sqrt{(\sum_{i = 0}^{N - 1} g_{i} h_{i} {\| x_{i} - \overline{x} \|}^{2}) (\sum_{i = 0}^{N - 1} g_{i} h_{i} {\| y_{i} - \overline{y} \|}^{2})}}$	`mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) sum(ghconj(x-mxe)(y-mye))/sqrt(sum(ghabs(x-mxe)2)sum(ghabs(y-mye)**2))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.conj(x-mxe).(y-mye))/sqrt(sum(g.h.abs(x-mxe).^2)sum(g.h.abs(y-mye).^2))`
Kreuzkorrelationsfunktion (erwartungstreu)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N} x_{i}^{} y_{(i + k) mod N}}{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `Ryx=zeros(N)+0j tau=arange(0,N)dt for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i])y[(i+k)%N] sum0+=g[i]h[(i+k)%N] Ryx[k]=sum1/float(sum0) plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=arange(0,N)dt X1=fft(gx) Y1=fft(hy) P1=conj(X1)Y1/N*2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N**2 Ryx=ifft(P1)/ifft(P0) plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()`	direkte Berechnung: `Ryx=zeros(1,N)+0j; tau=(0:N-1)dt; for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i))y(mod(i+k-1,N)+1); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end Ryx(k+1)=sum1/sum0; end plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; X1=fft(g.x); Y1=fft(h.y); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; Ryx=ifft(P1)./ifft(P0); plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y}^{} (- τ_{k}) = \frac{\underset{j - i \equiv k mod {\tilde{N}}_{x}}{\sum_{i = 0}^{N_{y} - 1} \sum_{j = 0}^{N - 1}} g_{i} h_{j} x_{i}^{} y_{j}}{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) R1=zeros(NT)+0j R0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): R1[(j-i)%NT]+=g[i]h[j]conj(x[i])y[j] R0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt Ryx=zeros(NT)+0j for k in range(0,NT): Ryx[k]=R1[k]/R0[k] plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i]x[i] wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i]y[i] wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) tau=arange(0,NT)dt P1=conj(X1)Y1/NT*2 P0=conj(X0)Y0/NT**2 Ryx=ifft(P1)/real(ifft(P0)) plot(tau,real(Ryx),'o',tau,imag(Ryx),'o') show()`	direkte Berechnung: `NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); R1=zeros(1,NT)+0j; R0=zeros(1,NT); for i=1:Nx for j=1:Ny R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i))y(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+g(i)h(j); end end tau=(0:NT-1)dt; Ryx=zeros(1,NT)+0j; for k=1:NT Ryx(k)=R1(k)/R0(k); end plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)x(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)y(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); tau=(0:NT-1)dt; P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; Ryx=ifft(P1)./real(ifft(P0)); plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Kreuzkovarianzfunktion (asymptotisch erwartungstreu)	$C_{y x, k} = C_{y x} (τ_{k}) = C_{x y}^{} (- τ_{k}) = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N} {(x_{i} - \overline{x})}^{} [y_{(i + k) mod N} - \overline{y}]}{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `Cyx=zeros(N)+0j tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i]-mxe)(y[(i+k)%N]-mye) sum0+=g[i]h[(i+k)%N] Cyx[k]=sum1/float(sum0) plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) X1=fft(g(x-mxe)) Y1=fft(h(y-mye)) P1=conj(X1)Y1/N*2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N**2 Cyx=ifft(P1)/ifft(P0) plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()`	direkte Berechnung: `Cyx=zeros(1,N)+0j; tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i)-mxe)(y(mod(i+k-1,N)+1)-mye); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end Cyx(k+1)=sum1/sum0; end plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); X1=fft(g.(x-mxe)); Y1=fft(h.(y-mye)); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; Cyx=ifft(P1)./ifft(P0); plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	$C_{y x, k} = C_{y x} (τ_{k}) = C_{x y}^{} (- τ_{k}) = \frac{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j} {(x_{i} - \overline{x})}^{} (y_{j} - \overline{y})}{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) C1=zeros(NT)+0j C0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): C1[(j-i)%NT]+=g[i]h[j]conj(x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt Cyx=zeros(NT)+0j for k in range(0,NT): Cyx[k]=C1[k]/C0[k] plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: from numpy.fft import mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i](x[i]-mxe) wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i](y[i]-mye) wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) tau=arange(0,NT)dt P1=conj(X1)Y1/NT*2 P0=conj(X0)Y0/NT**2 Cyx=ifft(P1)/real(ifft(P0)) plot(tau,real(Cyx),'o',tau,imag(Cyx),'o') show()	direkte Berechnung: `mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); C1=zeros(1,NT)+0j; C0=zeros(1,NT); for i=1:Nx for j=1:Ny C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i)-mxe)(y(j)-mye); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+g(i)h(j); end end tau=(0:NT-1)dt; Cyx=zeros(1,NT)+0j; for k=1:NT Cyx(k)=C1(k)/C0(k); end plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)(x(i)-mxe); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)(y(i)-mye); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); tau=(0:NT-1)dt; P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; Cyx=ifft(P1)./real(ifft(P0)); plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Kreuzkorrelationskoeffizientenfunktion (ohne Rauschen asymptotisch erwartungstreu, mit Rauschen systematische Fehler)	$ρ_{y x, k} = ρ_{y x} (τ_{k}) = ρ_{x y}^{} (- τ_{k}) = \frac{C_{y x, k}}{\sqrt{s_{x}^{2} s_{y}^{2}}} = \frac{\sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N} {(x_{i} - \overline{x})}^{} [y_{(i + k) mod N} - \overline{y}]}{\sqrt{s_{x}^{2} s_{y}^{2}} \cdot \sum_{i = 0}^{N - 1} g_{i} h_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `rhoyx=zeros(N)+0j tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)*2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i]-mxe)(y[(i+k)%N]-mye) sum0+=g[i]h[(i+k)%N] rhoyx[k]=sum1/float(sum0)/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)*2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 X1=fft(g(x-mxe)) Y1=fft(h(y-mye)) P1=conj(X1)Y1/N*2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N*2 rhoyx=ifft(P1)/ifft(P0)/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show()`	direkte Berechnung: `rhoyx=zeros(1,N)+0j; tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i)-mxe)(y(mod(i+k-1,N)+1)-mye); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end rhoyx(k+1)=sum1/sum0/sqrt(vxevye); end plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; X1=fft(g.(x-mxe)); Y1=fft(h.(y-mye)); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; rhoyx=ifft(P1)./ifft(P0)/sqrt(vxevye); plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	$ρ_{y x, k} = ρ_{y x} (τ_{k}) = ρ_{x y}^{} (- τ_{k}) = \frac{C_{y x, k}}{\sqrt{s_{x}^{2} s_{y}^{2}}} = \frac{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j} {(x_{i} - \overline{x})}^{} (y_{j} - \overline{y})}{\sqrt{s_{x}^{2} s_{y}^{2}} \cdot \underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} g_{i} h_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) C1=zeros(NT)+0j C0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): C1[(j-i)%NT]+=g[i]h[j]conj(x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt rhoyx=zeros(NT)+0j for k in range(0,NT): rhoyx[k]=C1[k]/C0[k]/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show() Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: from numpy.fft import mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gabs(x)2)/float(sum(g))-abs(sum(gx)/float(sum(g)))*2 vye=sum(habs(y)*2)/float(sum(h))-abs(sum(hy)/float(sum(h)))*2 NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i](x[i]-mxe) wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i](y[i]-mye) wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) tau=arange(0,NT)dt P1=conj(X1)Y1/NT2 P0=conj(X0)Y0/NT*2 rhoyx=ifft(P1)/real(ifft(P0))/sqrt(vxevye) plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') show()	direkte Berechnung: mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); C1=zeros(1,NT)+0j; C0=zeros(1,NT); for i=1:Nx for j=1:Ny C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i)-mxe)(y(j)-mye); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+g(i)h(j); end end tau=(0:NT-1)dt; rhoyx=zeros(1,NT)+0j; for k=1:NT rhoyx(k)=C1(k)/C0(k)/sqrt(vxevye); end plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o') Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.abs(x).^2)/sum(g)-abs(sum(g.x)/sum(g))^2; vye=sum(h.abs(y).^2)/sum(h)-abs(sum(h.y)/sum(h))^2; NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)(x(i)-mxe); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)(y(i)-mye); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); tau=(0:NT-1)dt; P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; rhoyx=ifft(P1)./real(ifft(P0))/sqrt(vxevye); plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Kreuzleistungsspektrum (erwartungstreu)	(wegen der nötigen Entfaltung mit dem Leistungsspektrum der Gewichte nur über die Korrelationsfunktion) $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{N} \cdot F F T {R_{y x, k}} = \frac{1}{N} \cdot \sum_{k = 0}^{N - 1} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / N}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{N Δ t}$ $j = - ⌊ N / 2 ⌋ \dots ⌊ (N - 1) / 2 ⌋$	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `from numpy.fft import * X1=fft(gx) Y1=fft(hy) P1=conj(X1)Y1/N2 X0=fft(g) Y0=fft(h) P0=conj(X0)Y0/N*2 Ryx=ifft(P1)/ifft(P0) f=roll(arange(-(N//2),N-(N//2))/float(Ndt),N-(N//2)) Pyx=fft(Ryx)/float(N) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import * Ryx=zeros(N)+0j for k in range(0,N): sum1=0+0j sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]conj(x[i])y[(i+k)%N] sum0+=g[i]h[(i+k)%N] Ryx[k]=sum1/float(sum0) f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2)) Pyx=fft(Ryx)/float(N) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `X1=fft(g.x); Y1=fft(h.y); P1=conj(X1).Y1/N^2; X0=fft(g); Y0=fft(h); P0=conj(X0).Y0/N^2; Ryx=ifft(P1)./ifft(P0); f=circshift(((-N/2:N-N/2-1))/(Ndt),[0;N/2]); Pyx=fft(Ryx)/N; plot(f,real(Pyx),'o',f,imag(Pyx),'o')` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `Ryx=zeros(1,N)+0j; for k=0:N-1 sum1=0+0j; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)conj(x(i))y(mod(i+k-1,N)+1); sum0=sum0+g(i)h(mod(i+k-1,N)+1); end Ryx(k+1)=sum1/sum0; end f=circshift(((-N/2:N-N/2-1))/(Ndt),[0;N/2]); Pyx=fft(Ryx)/N; plot(f,real(Pyx),'o',f,imag(Pyx),'o')`
Berechnung für den Fall, dass die Signallängen $N_{x}$ bzw. $N_{y}$ von der anzuwendende Periode des Signals $\tilde{N}$ abweichen ( $\tilde{N} \leq N_{x}$ , $\tilde{N} \leq N_{y}$ )	(wegen der nötigen Entfaltung mit dem Leistungsspektrum der Gewichte nur über die Korrelationsfunktion) $P_{y x, j} = P_{y x} (f_{j}) = \frac{1}{\tilde{N}} \cdot F F T {R_{y x, k}} = \frac{1}{\tilde{N}} \cdot \sum_{k = 0}^{\tilde{N} - 1} R_{y x, k} 𝐞^{- 2 π 𝐢 j k / \tilde{N}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{\tilde{N} Δ t}$ $j = - ⌊ \tilde{N} / 2 ⌋ \dots ⌊ (\tilde{N} - 1) / 2 ⌋$	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `from numpy.fft import * NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Nx): xp[i%NT]+=g[i]x[i] wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) yp=zeros(NT)+0j wp=zeros(NT) for i in range(0,Ny): yp[i%NT]+=h[i]y[i] wp[i%NT]+=h[i] Y1=fft(yp) Y0=fft(wp) P1=conj(X1)Y1/NT2 P0=conj(X0)Y0/NT*2 Ryx=ifft(P1)/real(ifft(P0)) f=roll(arange(-(NT//2),(NT+1)//2)/float(NTdt),(NT+1)//2) Pyx=fft(Ryx)/float(NT) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import * NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) R1=zeros(NT)+0j R0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): R1[(j-i)%NT]+=g[i]h[j]conj(x[i])y[j] R0[(j-i)%NT]+=g[i]h[j] Ryx=zeros(NT)+0j for k in range(0,NT): Ryx[k]=R1[k]/R0[k] f=roll(arange(-(NT//2),(NT+1)//2)/float(NT*dt),(NT+1)//2) Pyx=fft(Ryx)/float(NT) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+g(i)x(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+g(i); end X1=fft(xp); X0=fft(wp); yp=zeros(1,NT)+0j; wp=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+h(i)y(i); wp(mod(i-1,NT)+1)=wp(mod(i-1,NT)+1)+h(i); end Y1=fft(yp); Y0=fft(wp); P1=conj(X1).Y1/NT^2; P0=conj(X0).Y0/NT^2; Ryx=ifft(P1)./real(ifft(P0)); f=circshift((-fix(NT/2):fix((NT-1)/2))/(NTdt),[0;fix((NT+1)/2)]); Pyx=fft(Ryx)/NT; plot(f,real(Pyx),'o',f,imag(Pyx),'o') Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); R1=zeros(1,NT)+0j; R0=zeros(1,NT); for i=1:Nx for j=1:Ny R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+g(i)h(j)conj(x(i))y(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+g(i)h(j); end end Ryx=zeros(1,NT)+0j; for k=1:NT Ryx(k)=R1(k)/R0(k); end f=circshift((-fix(NT/2):fix((NT-1)/2))/(NTdt),[0;fix((NT+1)/2)]); Pyx=fft(Ryx)/NT; plot(f,real(Pyx),'o',f,imag(Pyx),'o')`