Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für 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 $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')`
Mittelwerte (ohne Rauschen exakt, mit Rauschen erwartungstreu)	$\overline{x} = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} x_{i}$ $\overline{y} = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} y_{i}$	`mxe=mean(x) mye=mean(y)`	`mxe=mean(x); mye=mean(y);`
Varianzen (für streng periodische Signale keine Bessel-Korrektur, da Mittelwert keine Schätzgröße, ohne Rauschen exakt, mit Rauschen systematischer Fehler)	$s_{x}^{2} = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} {(x_{i} - \overline{x})}^{2} = (\frac{1}{N} \cdot \sum_{i = 0}^{N - 1} x_{i}^{2}) - {\overline{x}}^{2}$ $s_{y}^{2} = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} {(y_{i} - \overline{y})}^{2} = (\frac{1}{N} \cdot \sum_{i = 0}^{N - 1} y_{i}^{2}) - {\overline{y}}^{2}$	`vxe=var(x) vye=var(y)`	`vxe=var(x,1); vye=var(y,1);`
Kreuzkovarianz (für streng periodische Signale keine Bessel-Korrektur, da Mittelwert keine Schätzgröße, ohne Rauschen exakt, mit Rauschen erwartungstreu)	$c_{y x} = c_{x y} = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} (x_{i} - \overline{x}) (y_{i} - \overline{y}) = (\frac{1}{N} \cdot \sum_{i = 0}^{N - 1} x_{i} y_{i}) - \overline{x} \cdot \overline{y}$	`cov(x,y,ddof=0)[0,1]`	`sum(x.y)/N-mean(x)mean(y)`
(Kreuz-)Korrelationskoeffizient (für streng periodische Signale keine Bessel-Korrektur, da Mittelwert keine Schätzgröße, ohne Rauschen exakt, 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} (x_{i} - \overline{x}) (y_{i} - \overline{y})}{\sqrt{[\sum_{i = 0}^{N - 1} {(x_{i} - \overline{x})}^{2}] [\sum_{i = 0}^{N - 1} {(y_{i} - \overline{y})}^{2}]}}$ $= \frac{(\sum_{i = 0}^{N - 1} x_{i} y_{i}) - N \overline{x} \cdot \overline{y}}{\sqrt{[(\sum_{i = 0}^{N - 1} x_{i}^{2}) - N {\overline{x}}^{2}] [(\sum_{i = 0}^{N - 1} y_{i}^{2}) - N {\overline{y}}^{2}]}}$	`cov(x,y,ddof=0)[0,1]/sqrt(var(x)*var(y))`	`(sum(x.y)/N-mean(x)mean(y))/sqrt(var(x,1)*var(y,1))`
Kreuzkorrelationsfunktion (ohne Rauschen exakt, mit Rauschen erwartungstreu)	$R_{y x, k} = R_{y x} (τ_{k}) = R_{x y} (- τ_{k}) = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} x_{i} y_{(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 for i in range(0,N): sum1+=x[i]y[(i+k)%N] Ryx[k]=sum1/float(N) plot(tau,Ryx,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * tau=arange(0,N)dt X=fft(x) Y=fft(y) Pyx=conj(X)Y/N*2 Ryx=Nreal(ifft(Pyx)) plot(tau,Ryx,'o') show()`	direkte Berechnung: `Ryx=zeros(1,N); tau=(0:N-1)dt; for k=0:N-1 sum1=0; for i=1:N sum1=sum1+x(i)y(mod(i+k-1,N)+1); end Ryx(k+1)=sum1/N; end plot(tau,Ryx,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; X=fft(x); Y=fft(y); Pyx=conj(X).Y/N^2; Ryx=N*real(ifft(Pyx)); 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}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} x_{i} y_{j}}{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1}} 1}$ $τ_{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]+=x[i]y[j] R0[(j-i)%NT]+=1 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) xp=zeros(NT) np=zeros(NT) for i in range(0,Nx): xp[i%NT]+=x[i] np[i%NT]+=1 X1=fft(xp) X0=fft(np) Ny=len(x) yp=zeros(NT) np=zeros(NT) for i in range(0,Ny): yp[i%NT]+=y[i] np[i%NT]+=1 Y1=fft(yp) Y0=fft(np) 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)+x(i)y(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+1; 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); np=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i); np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end tau=(0:NT-1)dt; X1=fft(xp); X0=fft(np); yp=zeros(1,NT); np=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i); np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end Y1=fft(yp); Y0=fft(np); 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 (ohne Rauschen exakt, mit Rauschen erwartungstreu)	$C_{y x, k} = C_{y x} (τ_{k}) = C_{x y} (- τ_{k}) = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} (x_{i} - \overline{x}) [y_{(i + k) mod N} - \overline{y}]$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `Cyx=zeros(N) tau=arange(0,N)dt mxe=mean(x) mye=mean(y) for k in range(0,N): sum1=0 for i in range(0,N): sum1+=(x[i]-mxe)(y[(i+k)%N]-mye) Cyx[k]=sum1/float(N) plot(tau,Cyx,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * tau=arange(0,N)dt mxe=mean(x) mye=mean(y) X=fft(x-mxe) Y=fft(y-mye) Pyx=conj(X)Y/N*2 Cyx=Nreal(ifft(Pyx)) plot(tau,Cyx,'o') show()`	direkte Berechnung: `Cyx=zeros(1,N); tau=(0:N-1)dt; mxe=mean(x); mye=mean(y); for k=0:N-1 sum1=0; for i=1:N sum1=sum1+(x(i)-mxe)(y(mod(i+k-1,N)+1)-mye); end Cyx(k+1)=sum1/N; end plot(tau,Cyx,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=mean(x); mye=mean(y); X=fft(x-mxe); Y=fft(y-mye); Pyx=conj(X).Y/N^2; Cyx=N*real(ifft(Pyx)); 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}} (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}} 1}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `mxe=mean(x) mye=mean(y) 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]+=(x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=1 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=mean(x) mye=mean(y) NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT) np=zeros(NT) for i in range(0,Nx): xp[i%NT]+=x[i]-mxe np[i%NT]+=1 X1=fft(xp) X0=fft(np) yp=zeros(NT) np=zeros(NT) for i in range(0,Ny): yp[i%NT]+=y[i]-mye np[i%NT]+=1 Y1=fft(yp) Y0=fft(np) 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=mean(x); mye=mean(y); 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)+(x(i)-mxe)(y(j)-mye); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+1; 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=mean(x); mye=mean(y); NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT); np=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i)-mxe; np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end X1=fft(xp); X0=fft(np); yp=zeros(1,NT); np=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i)-mye; np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end Y1=fft(yp); Y0=fft(np); 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 exakt, 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{1}{N \sqrt{s_{x}^{2} s_{y}^{2}}} \cdot \sum_{i = 0}^{N - 1} (x_{i} - \overline{x}) [y_{(i + k) mod N} - \overline{y}]$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `rhoyx=zeros(N) tau=arange(0,N)dt mxe=mean(x) mye=mean(y) vxe=var(x) vye=var(y) for k in range(0,N): sum1=0 for i in range(0,N): sum1+=(x[i]-mxe)(y[(i+k)%N]-mye) rhoyx[k]=sum1/float(N)/sqrt(vxevye) plot(tau,rhoyx,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=arange(0,N)dt mxe=mean(x) mye=mean(y) vxe=var(x) vye=var(y) X=fft(x-mxe) Y=fft(y-mye) Pyx=conj(X)Y/N*2 rhoyx=Nreal(ifft(Pyx))/sqrt(vxe*vye) plot(tau,rhoyx,'o') show()`	direkte Berechnung: `rhoyx=zeros(1,N); tau=(0:N-1)dt; mxe=mean(x); mye=mean(y); vxe=var(x,1); vye=var(y,1); for k=0:N-1 sum1=0; for i=1:N sum1=sum1+(x(i)-mxe)(y(mod(i+k-1,N)+1)-mye); end rhoyx(k+1)=sum1/N/sqrt(vxevye); end plot(tau,rhoyx,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=mean(x); mye=mean(y); vxe=var(x,1); vye=var(y,1); X=fft(x-mxe); Y=fft(y-mye); Pyx=conj(X).Y/N^2; rhoyx=Nreal(ifft(Pyx))/sqrt(vxe*vye); 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}} (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}} 1}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `mxe=mean(x) mye=mean(y) vxe=var(x) vye=var(y) 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]+=(x[i]-mxe)(y[j]-mye) C0[(j-i)%NT]+=1 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=mean(x) mye=mean(y) vxe=var(x) vye=var(y) NT=60 #Periodenlänge NT≤N Nx=len(x) Ny=len(y) xp=zeros(NT) np=zeros(NT) for i in range(0,Nx): xp[i%NT]+=x[i]-mxe np[i%NT]+=1 X1=fft(xp) X0=fft(np) yp=zeros(NT) np=zeros(NT) for i in range(0,Ny): yp[i%NT]+=y[i]-mye np[i%NT]+=1 Y1=fft(yp) Y0=fft(np) tau=arange(0,NT)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') show()`	direkte Berechnung: `mxe=mean(x); mye=mean(y); vxe=var(x,1); vye=var(y,1); 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)+(x(i)-mxe)(y(j)-mye); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+1; 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=mean(x); mye=mean(y); vxe=var(x,1); vye=var(y,1); NT=60; %Periodenlänge NT≤N Nx=length(x); Ny=length(y); xp=zeros(1,NT); np=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i)-mxe; np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end X1=fft(xp); X0=fft(np); yp=zeros(1,NT); np=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i)-mye; np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end Y1=fft(yp); Y0=fft(np); 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(vxe*vye); plot(tau,rhoyx,'o')
Kreuzleistungsspektrum (ohne Rauschen exakt, mit Rauschen erwartungstreu)	$X_{j} = F F T {x_{i}} = \sum_{i = 0}^{N - 1} x_{i} 𝐞^{- 2 π 𝐢 i j / N}$ $Y_{j} = F F T {y_{i}} = \sum_{i = 0}^{N - 1} y_{i} 𝐞^{- 2 π 𝐢 i j / N}$ (imaginäre Einheit $𝐢$ ) $P_{y x, j} = P_{y x} (f_{j}) = \frac{X_{j}^{*} Y_{j}}{N^{2}}$ $f_{j} = \frac{j}{N Δ t}$ $j = - ⌊ N / 2 ⌋ \dots ⌊ (N - 1) / 2 ⌋$	direkte Berechnung: `from numpy.fft import * f=roll(arange(-(N//2),(N+1)//2)/float(Ndt),(N+1)//2) X=fft(x) Y=fft(y) Pyx=conj(X)Y/N*2 plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()` Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import Ryx=zeros(N) for k in range(0,N): sum1=0 for i in range(0,N): sum1+=x[i]y[(i+k)%N] Ryx[k]=sum1/N f=roll(arange(-(N//2),(N+1)//2)/float(Ndt),(N+1)//2) Pyx=fft(Ryx)/float(N) plot(f,real(Pyx),'o',f,imag(Pyx),'o') show()`	direkte Berechnung: `f=circshift((-fix(N/2):fix((N-1)/2))/(Ndt),[0;fix((N+1)/2)]); X=fft(x); Y=fft(y); Pyx=conj(X).Y/N^2; plot(f,real(Pyx),'o',f,imag(Pyx),'o')` Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `Ryx=zeros(1,N); for k=0:N-1 sum1=0; for i=1:N sum1=sum1+x(i)y(mod(i+k-1,N)+1); end Ryx(k+1)=sum1/N; end f=circshift((-fix(N/2):fix((N-1)/2))/(Ndt),[0;fix((N+1)/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 Überlappungen 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) np=zeros(NT) for i in range(0,Nx): xp[i%NT]+=x[i] np[i%NT]+=1 X1=fft(xp) X0=fft(np) yp=zeros(NT) np=zeros(NT) for i in range(0,Ny): yp[i%NT]+=y[i] np[i%NT]+=1 Y1=fft(yp) Y0=fft(np) 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]+=x[i]y[j] R0[(j-i)%NT]+=1 Ryx=zeros(NT) for k in range(0,NT): Ryx[k]=R1[k]/R0[k] 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 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); np=zeros(1,NT); for i=1:Nx xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i); np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end X1=fft(xp); X0=fft(np); yp=zeros(1,NT); np=zeros(1,NT); for i=1:Ny yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i); np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1; end Y1=fft(yp); Y0=fft(np); 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)+x(i)y(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+1; 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))/(NT*dt),[0;fix((NT+1)/2)]); Pyx=fft(Ryx)/NT; plot(f,real(Pyx),'o',f,imag(Pyx),'o')`