Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für reelle, 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$ reellen Wertepaaren $(x_{i}, y_{i})$ mit $i = 0 \dots N - 1$ (zwei harmonische Signal als Beispiel, mit einem Erwartungswert verschieden von null, ohne Rauschen)		`N=100 dt=1.0 t=arange(0,N)dt x=3+1.5sin(0.1pit) y=2+1sin(0.1pi*t-0.3) plot(t,x,'o',t,y,'o') show()`	`N=100; dt=1.0; t=(0:N-1)dt; x=3+1.5sin(0.1pit); y=2+1sin(0.1pi*t-0.3); plot(t,x,'o',t,y,'o')`
Überlagerung von Rauschen		`x+=normal(0,0.1,N) y+=normal(0,0.1,N) plot(t,x,'o',t,y,'o') show()`	`x=x+0.1randn(1,N); y=y+0.1randn(1,N); plot(t,x,'o',t,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})}^{2}}{\sum_{i = 0}^{N - 1} g_{i}} = \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})}^{2}}{\sum_{i = 0}^{N - 1} h_{i}} = \frac{\sum_{i = 0}^{N - 1} h_{i} y_{i}^{2}}{\sum_{i = 0}^{N - 1} h_{i}} - {\overline{y}}^{2}$	`vxe=sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 vye=sum(hy*2)/float(sum(h))-(sum(hy)/float(sum(h)))**2`	`vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; vye=sum(h.y.^2)/sum(h)-(sum(h.y)/sum(h))^2;`
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(gh(x-mxe)(y-mye))/float(sum(gh))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.(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(gh(x-mxe)(y-mye))sqrt(sum(g)sum(h))/(sum(gh)sqrt(sum(g(x-mxe)*2)sum(h(y-mye)*2)))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.(x-mxe).(y-mye))sqrt(sum(g)sum(h))/(sum(g.h)sqrt(sum(g.(x-mxe).^2)sum(h.(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(gh(x-mxe)(y-mye))/sqrt(sum(gh(x-mxe)2)sum(gh(y-mye)**2))`	`mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); sum(g.h.(x-mxe).(y-mye))/sqrt(sum(g.h.(x-mxe).^2)sum(g.h.(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) tau=arange(0,N)dt for k in range(0,N): sum1=0 sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]x[i]y[(i+k)%N] sum0+=g[i]h[(i+k)%N] Ryx[k]=sum1/float(sum0) plot(tau,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=real(ifft(P1))/real(ifft(P0)) plot(tau,Ryx,'o') show()`	direkte Berechnung: `Ryx=zeros(1,N); tau=(0:N-1)dt; for k=0:N-1 sum1=0; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)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,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=real(ifft(P1))./real(ifft(P0)); plot(tau,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) R0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): R1[(j-i)%NT]+=g[i]h[j]x[i]y[j] R0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt Ryx=zeros(NT) for k in range(0,NT): Ryx[k]=R1[k]/R0[k] plot(tau,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) 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) 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=real(ifft(P1))/real(ifft(P0)) plot(tau,Ryx,'o') show()`	direkte Berechnung: `NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); R1=zeros(1,NT); 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)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); for k=1:NT Ryx(k)=R1(k)/R0(k); end plot(tau,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); 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); 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=real(ifft(P1))./real(ifft(P0)); plot(tau,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) 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 sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N](x[i]-mxe)(y[(i+k)%N]-mye) sum0+=g[i]h[(i+k)%N] Cyx[k]=sum1/float(sum0) plot(tau,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=real(ifft(P1))/real(ifft(P0)) plot(tau,Cyx,'o') show()`	direkte Berechnung: `Cyx=zeros(1,N); tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); for k=0:N-1 sum1=0; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)(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,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=real(ifft(P1))./real(ifft(P0)); plot(tau,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) C0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): C1[(j-i)%NT]+=g[i]h[j](x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt Cyx=zeros(NT) for k in range(0,NT): Cyx[k]=C1[k]/C0[k] plot(tau,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) 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) 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=real(ifft(P1))/real(ifft(P0)) plot(tau,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); 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)(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); for k=1:NT Cyx(k)=C1(k)/C0(k); end plot(tau,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); 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); 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=real(ifft(P1))./real(ifft(P0)); plot(tau,Cyx,'o')
Kreuzkorrelationskoeffizientenfunktion (ohne Rauschen asymptotisch erwartungstreu, mit Rauschen systematische Fehle)	$ρ_{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) tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) mye=sum(hy)/float(sum(h)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 vye=sum(hy*2)/float(sum(h))-(sum(hy)/float(sum(h)))*2 for k in range(0,N): sum1=0 sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N](x[i]-mxe)(y[(i+k)%N]-mye) sum0+=g[i]h[(i+k)%N] rhoyx[k]=sum1/float(sum0)/sqrt(vxevye) plot(tau,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(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 vye=sum(hy*2)/float(sum(h))-(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=real(ifft(P1)/ifft(P0))/sqrt(vxevye) plot(tau,rhoyx,'o') show()`	direkte Berechnung: `rhoyx=zeros(1,N); tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; vye=sum(h.y.^2)/sum(h)-(sum(h.y)/sum(h))^2; for k=0:N-1 sum1=0; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)(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,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.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; vye=sum(h.y.^2)/sum(h)-(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=real(ifft(P1)./ifft(P0))/sqrt(vxevye); plot(tau,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(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 vye=sum(hy*2)/float(sum(h))-(sum(hy)/float(sum(h)))*2 NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) C1=zeros(NT) C0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): C1[(j-i)%NT]+=g[i]h[j](x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=g[i]h[j] tau=arange(0,NT)dt rhoyx=zeros(NT) for k in range(0,NT): rhoyx[k]=C1[k]/C0[k]/sqrt(vxevye) plot(tau,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(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 vye=sum(hy*2)/float(sum(h))-(sum(hy)/float(sum(h)))*2 NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT) 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) 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=real(ifft(P1))/real(ifft(P0))/sqrt(vxevye) plot(tau,rhoyx,'o') show()	direkte Berechnung: `mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; vye=sum(h.y.^2)/sum(h)-(sum(h.y)/sum(h))^2; NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); C1=zeros(1,NT); 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)(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); for k=1:NT rhoyx(k)=C1(k)/C0(k)/sqrt(vxevye); end plot(tau,rhoyx,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: mxe=sum(g.x)/sum(g); mye=sum(h.y)/sum(h); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; vye=sum(h.y.^2)/sum(h)-(sum(h.y)/sum(h))^2; NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT); 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); 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=real(ifft(P1))./real(ifft(P0))/sqrt(vxevye); plot(tau,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=real(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) for k in range(0,N): sum1=0 sum0=0 for i in range(0,N): sum1+=g[i]h[(i+k)%N]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=real(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); for k=0:N-1 sum1=0; sum0=0; for i=1:N sum1=sum1+g(i)h(mod(i+k-1,N)+1)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) 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) 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=real(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) R0=zeros(NT) for i in range(0,Nx): for j in range(0,Ny): R1[(j-i)%NT]+=g[i]h[j]x[i]y[j] R0[(j-i)%NT]+=g[i]h[j] Ryx=zeros(NT) 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); 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); 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=real(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); 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)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); 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')`