Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für reelle, periodische Einzelsignale, mit Gewichtung

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von $N$ Werten $x_{i}$ mit $i = 0 \dots N - 1$ (Superposition zweier harmonischer Signale als Beispiel, mit einem Erwartungswert verschieden von null)		`N=100 dt=1.0 t=arange(0,N)dt x=3+1.5sin(0.2pit)+1.5sin(0.1pi*t) plot(t,x,'o') show()`	`N=100; dt=1.0; t=(0:N-1)dt; x=3+1.5sin(0.2pit)+1.5sin(0.1pi*t); plot(t,x,'o')`
Überlagerung von Rauschen		`x+=normal(0,0.1,N) plot(t,x,'o') show()`	`x=x+0.1*randn(1,N); plot(t,x,'o')`
Generieren von $N$ Gewichten $g_{i}$ mit $i = 0 \dots N - 1$ , passend zu den $x_{i}$ mit $i = 0 \dots N - 1$ (Kosinusfenster als Beispiel)		`g=sin(pi*(arange(0,N)+0.5)/float(N)) plot(t,g,'o') show()`	`g=sin(pi*((0:N-1)+0.5)/N); plot(t,g,'o')`
Mittelwert (erwartungstreu)	$\overline{x} = \frac{\sum_{i = 0}^{N - 1} g_{i} x_{i}}{\sum_{i = 0}^{N - 1} g_{i}}$	`sum(g*x)/float(sum(g))`	`sum(g.*x)/sum(g)`
Varianz (ohne Rauschen asymptotisch erwartungstreu, da selbst für streng periodische Signale der Mittelwert eine Schätzgröße ist, mit Rauschen systematische Fehler)	$s^{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}$	`sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))**2`	`sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2`
(Auto-)Korrelationsfunktion (ohne Rauschen erwartungstreu, mit Rauschen systematischer Fehler bei $τ_{k} = 0$ )	$R_{k} = R (τ_{k}) = \frac{\sum_{i = 0}^{N - 1} g_{i} g_{(i + k) mod N} x_{i} x_{(i + k) mod N}}{\sum_{i = 0}^{N - 1} g_{i} g_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `R=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]g[(i+k)%N]x[i]x[(i+k)%N] sum0+=g[i]g[(i+k)%N] R[k]=sum1/float(sum0) plot(tau,R,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=arange(0,N)dt X1=fft(gx) P1=abs(X1)2/N2 X0=fft(g) P0=abs(X0)2/N2 R=real(ifft(P1))/real(ifft(P0)) plot(tau,R,'o') show()`	direkte Berechnung: `R=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)g(mod(i+k-1,N)+1)x(i)x(mod(i+k-1,N)+1); sum0=sum0+g(i)g(mod(i+k-1,N)+1); end R(k+1)=sum1/sum0; end plot(tau,R,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; X1=fft(g.*x); P1=abs(X1).^2/N^2; X0=fft(g); P0=abs(X0).^2/N^2; R=real(ifft(P1))./real(ifft(P0)); plot(tau,R,'o')`
Berechnung für den Fall, dass die anzuwendende Periode des Signals $\tilde{N}$ von der Signallänge $N$ abweicht ( $\tilde{N} \leq N$ )	$R_{k} = R (τ_{k}) = \frac{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1}} g_{i} g_{j} x_{i} x_{j}}{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1}} g_{i} g_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `NT=60 #Periodenlänge NT≤N R1=zeros(NT) R0=zeros(NT) for i in range(0,N): for j in range(0,N): R1[(j-i)%NT]+=g[i]g[j]x[i]x[j] R0[(j-i)%NT]+=g[i]g[j] tau=arange(0,NT)dt R=zeros(NT) for k in range(0,NT): R[k]=R1[k]/R0[k] plot(tau,R,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import NT=60 #Periodenlänge NT≤N xp=zeros(NT) wp=zeros(NT) for i in range(0,N): xp[i%NT]+=g[i]x[i] wp[i%NT]+=g[i] tau=arange(0,NT)dt X1=fft(xp) X0=fft(wp) P1=abs(X1)2/NT2 P0=abs(X0)2/NT2 R=real(ifft(P1))/real(ifft(P0)) plot(tau,R,'o') show()`	direkte Berechnung: `NT=60; %Periodenlänge NT≤N R1=zeros(1,NT); R0=zeros(1,NT); for i=1:N for j=1:N R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+g(i)g(j)x(i)x(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+g(i)g(j); end end tau=(0:NT-1)dt; R=zeros(1,NT); for k=1:NT R(k)=R1(k)/R0(k); end plot(tau,R,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `NT=60; %Periodenlänge NT≤N xp=zeros(1,NT); wp=zeros(1,NT); for i=1:N 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 tau=(0:NT-1)*dt; X1=fft(xp); X0=fft(wp); P1=abs(X1).^2/NT^2; P0=abs(X0).^2/NT^2; R=real(ifft(P1))./real(ifft(P0)); plot(tau,R,'o')`
(Auto-)Kovarianzfunktion (ohne Rauschen asymptotisch erwartungstreu,mit Rauschen systematischer Fehler bei $τ_{k} = 0$ und bei $τ_{k} \neq 0$ asymptotisch erwatungstreu)	$C_{k} = C (τ_{k}) = \frac{\sum_{i = 0}^{N - 1} g_{i} g_{(i + k) mod N} (x_{i} - \overline{x}) [x_{(i + k) mod N} - \overline{x}]}{\sum_{i = 0}^{N - 1} g_{i} g_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `C=zeros(N) tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) for k in range(0,N): sum1=0 sum0=0 for i in range(0,N): sum1+=g[i]g[(i+k)%N](x[i]-mxe)(x[(i+k)%N]-mxe) sum0+=g[i]g[(i+k)%N] C[k]=sum1/float(sum0) plot(tau,C,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) X1=fft(g(x-mxe)) P1=abs(X1)2/N2 X0=fft(g) P0=abs(X0)2/N*2 C=real(ifft(P1))/real(ifft(P0)) plot(tau,C,'o') show()`	direkte Berechnung: `C=zeros(1,N); tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); for k=0:N-1 sum1=0; sum0=0; for i=1:N sum1=sum1+g(i)g(mod(i+k-1,N)+1)(x(i)-mxe)(x(mod(i+k-1,N)+1)-mxe); sum0=sum0+g(i)g(mod(i+k-1,N)+1); end C(k+1)=sum1/sum0; end plot(tau,C,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); X1=fft(g.*(x-mxe)); P1=abs(X1).^2/N^2; X0=fft(g); P0=abs(X0).^2/N^2; C=real(ifft(P1))./real(ifft(P0)); plot(tau,C,'o')`
Berechnung für den Fall, dass die anzuwendende Periode des Signals $\tilde{N}$ von der Signallänge $N$ abweicht ( $\tilde{N} \leq N$ )	$C_{k} = C (τ_{k}) = \frac{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1}} g_{i} g_{j} (x_{i} - \overline{x}) (x_{j} - \overline{x})}{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1}} g_{i} g_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `mxe=sum(gx)/float(sum(g)) NT=60 #Periodenlänge NT≤N C1=zeros(NT) C0=zeros(NT) for i in range(0,N): for j in range(0,N): C1[(j-i)%NT]+=g[i]g[j](x[i]-mxe)(x[j]-mxe) C0[(j-i)%NT]+=g[i]g[j] tau=arange(0,NT)dt C=zeros(NT) for k in range(0,NT): C[k]=C1[k]/C0[k] plot(tau,C,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * mxe=sum(gx)/float(sum(g)) NT=60 #Periodenlänge NT≤N xp=zeros(NT) wp=zeros(NT) for i in range(0,N): xp[i%NT]+=g[i](x[i]-mxe) wp[i%NT]+=g[i] tau=arange(0,NT)dt X1=fft(xp) X0=fft(wp) P1=abs(X1)2/NT2 P0=abs(X0)2/NT*2 C=real(ifft(P1))/real(ifft(P0)) plot(tau,C,'o') show()`	direkte Berechnung: `mxe=sum(g.x)/sum(g); NT=60; %Periodenlänge NT≤N C1=zeros(1,NT); C0=zeros(1,NT); for i=1:N for j=1:N C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+g(i)g(j)(x(i)-mxe)(x(j)-mxe); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+g(i)g(j); end end tau=(0:NT-1)dt; C=zeros(1,NT); for k=1:NT C(k)=C1(k)/C0(k); end plot(tau,C,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `mxe=sum(g.x)/sum(g); NT=60; %Periodenlänge NT≤N xp=zeros(1,NT); wp=zeros(1,NT); for i=1:N 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 tau=(0:NT-1)*dt; X1=fft(xp); X0=fft(wp); P1=abs(X1).^2/NT^2; P0=abs(X0).^2/NT^2; C=real(ifft(P1))./real(ifft(P0)); plot(tau,C,'o')`
(Auto-)Korrelationskoeffizientenfunktion (ohne Rauschen asymptotisch erwartungstreu, mit Rauschen systematischer Fehler bei allen $τ_{k} \neq 0$ )	$ρ_{k} = ρ (τ_{k}) = \frac{C_{k}}{s^{2}} = \frac{\sum_{i = 0}^{N - 1} g_{i} g_{(i + k) mod N} (x_{i} - \overline{x}) [x_{(i + k) mod N} - \overline{x}]}{s^{2} \cdot \sum_{i = 0}^{N - 1} g_{i} g_{(i + k) mod N}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `rho=zeros(N) tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) vxe=sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 for k in range(0,N): sum1=0 sum0=0 for i in range(0,N): sum1+=g[i]g[(i+k)%N](x[i]-mxe)(x[(i+k)%N]-mxe) sum0+=g[i]g[(i+k)%N] rho[k]=sum1/float(sum0)/vxe plot(tau,rho,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=arange(0,N)dt mxe=sum(gx)/float(sum(g)) vxe=sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 X1=fft(g(x-mxe)) P1=abs(X1)2/N2 X0=fft(g) P0=abs(X0)2/N2 rho=real(ifft(P1))/real(ifft(P0))/vxe plot(tau,rho,'o') show()`	direkte Berechnung: `rho=zeros(1,N); tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; for k=0:N-1 sum1=0; sum0=0; for i=1:N sum1=sum1+g(i)g(mod(i+k-1,N)+1)(x(i)-mxe)(x(mod(i+k-1,N)+1)-mxe); sum0=sum0+g(i)g(mod(i+k-1,N)+1); end rho(k+1)=sum1/sum0/vxe; end plot(tau,rho,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `tau=(0:N-1)dt; mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; X1=fft(g.*(x-mxe)); P1=abs(X1).^2/N^2; X0=fft(g); P0=abs(X0).^2/N^2; rho=real(ifft(P1))./real(ifft(P0))/vxe; plot(tau,rho,'o')`
Berechnung für den Fall, dass die anzuwendende Periode des Signals $\tilde{N}$ von der Signallänge $N$ abweicht ( $\tilde{N} \leq N$ )	$ρ_{k} = ρ (τ_{k}) = \frac{C_{k}}{s^{2}} = \frac{\underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1}} g_{i} g_{j} (x_{i} - \overline{x}) (x_{j} - \overline{x})}{s^{2} \cdot \underset{j - i \equiv k mod \tilde{N}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1}} g_{i} g_{j}}$ $τ_{k} = k Δ t$ $k = - \infty \dots \infty$	direkte Berechnung: `mxe=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 NT=60 #Periodenlänge NT≤N C1=zeros(NT) C0=zeros(NT) for i in range(0,N): for j in range(0,N): C1[(j-i)%NT]+=g[i]g[j](x[i]-mxe)(x[j]-mxe) C0[(j-i)%NT]+=g[i]g[j] tau=arange(0,NT)dt rho=zeros(NT) for k in range(0,NT): rho[k]=C1[k]/C0[k]/vxe plot(tau,rho,'o') show()` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * mxe=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 NT=60 #Periodenlänge NT≤N xp=zeros(NT) wp=zeros(NT) for i in range(0,N): xp[i%NT]+=g[i](x[i]-mxe) wp[i%NT]+=g[i] tau=arange(0,NT)dt X1=fft(xp) X0=fft(wp) P1=abs(X1)2/NT2 P0=abs(X0)2/NT*2 rho=real(ifft(P1))/real(ifft(P0))/vxe plot(tau,rho,'o') show()`	direkte Berechnung: `mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; NT=60; %Periodenlänge NT≤N C1=zeros(1,NT); C0=zeros(1,NT); for i=1:N for j=1:N C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+g(i)g(j)(x(i)-mxe)(x(j)-mxe); C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+g(i)g(j); end end tau=(0:NT-1)dt; rho=zeros(1,NT); for k=1:NT rho(k)=C1(k)/C0(k)/vxe; end plot(tau,rho,'o')` Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem: `mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; NT=60; %Periodenlänge NT≤N xp=zeros(1,NT); wp=zeros(1,NT); for i=1:N 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 tau=(0:NT-1)*dt; X1=fft(xp); X0=fft(wp); P1=abs(X1).^2/NT^2; P0=abs(X0).^2/NT^2; rho=real(ifft(P1))./real(ifft(P0))/vxe; plot(tau,rho,'o')`
(Auto-)Leistungsspektrum (ohne Rauschen erwartungstreu, mit Rauschen systematischer Fehler bei allen Frequenzen)	(wegen der nötigen Entfaltung mit dem Leistungsspektrum der Gewichte nur über die Korrelationsfunktion) $P_{j} = P (f_{j}) = \frac{1}{N} \cdot F F T {R_{k}} = \frac{1}{N} \cdot \sum_{k = 0}^{N - 1} R_{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) P1=abs(X1)2/N2 X0=fft(g) P0=abs(X0)2/N2 R=real(ifft(P1))/real(ifft(P0)) f=roll(arange(-(N//2),(N+1)//2)/float(Ndt),(N+1)//2) P=real(fft(R))/float(N) plot(f,P,'o') show()` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import * R=zeros(N) for k in range(0,N): sum1=0 sum0=0 for i in range(0,N): sum1+=g[i]g[(i+k)%N]x[i]x[(i+k)%N] sum0+=g[i]g[(i+k)%N] R[k]=sum1/float(sum0) f=roll(arange(-(N//2),(N+1)//2)/float(N*dt),(N+1)//2) P=real(fft(R))/float(N) plot(f,P,'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `X1=fft(g.x); P1=abs(X1).^2/N^2; X0=fft(g); P0=abs(X0).^2/N^2; R=real(ifft(P1))./real(ifft(P0)); f=circshift((-fix(N/2):fix((N-1)/2))/(Ndt),[0;fix((N+1)/2)]); P=real(fft(R))/N; plot(f,P,'o')` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `R=zeros(1,N); for k=0:N-1 sum1=0; sum0=0; for i=1:N sum1=sum1+g(i)g(mod(i+k-1,N)+1)x(i)x(mod(i+k-1,N)+1); sum0=sum0+g(i)g(mod(i+k-1,N)+1); end R(k+1)=sum1/sum0; end f=circshift((-fix(N/2):fix((N-1)/2))/(N*dt),[0;fix((N+1)/2)]); P=real(fft(R))/N; plot(f,P,'o')`
Berechnung für den Fall, dass die anzuwendende Periode des Signals $\tilde{N}$ von der Signallänge $N$ abweicht ( $\tilde{N} \leq N$ )	(wegen der nötigen Entfaltung mit dem Leistungsspektrum der Gewichte nur über die Korrelationsfunktion) $P_{j} = P (f_{j}) = \frac{1}{\tilde{N}} \cdot F F T {R_{k}} = \frac{1}{\tilde{N}} \cdot \sum_{k = 0}^{\tilde{N} - 1} R_{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 xp=zeros(NT) wp=zeros(NT) for i in range(0,N): xp[i%NT]+=g[i]x[i] wp[i%NT]+=g[i] X1=fft(xp) X0=fft(wp) P1=abs(X1)2/NT2 P0=abs(X0)2/NT2 R=real(ifft(P1))/real(ifft(P0)) f=roll(arange(-(NT//2),(NT+1)//2)/float(NTdt),(NT+1)//2) P=real(fft(R))/float(NT) plot(f,P,'o') show()` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import * NT=60 #Periodenlänge NT≤N R1=zeros(NT) R0=zeros(NT) for i in range(0,N): for j in range(0,N): R1[(j-i)%NT]+=g[i]g[j]x[i]x[j] R0[(j-i)%NT]+=g[i]g[j] R=zeros(NT) for k in range(0,NT): R[k]=R1[k]/R0[k] f=roll(arange(-(NT//2),(NT+1)//2)/float(NT*dt),(NT+1)//2) P=real(fft(R))/float(NT) plot(f,P,'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `NT=60; %Periodenlänge NT≤N xp=zeros(1,NT); wp=zeros(1,NT); for i=1:N 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); P1=abs(X1).^2/NT^2; P0=abs(X0).^2/NT^2; R=real(ifft(P1))./real(ifft(P0)); f=circshift((-fix(NT/2):fix((NT-1)/2))/(NTdt),[0;fix((NT+1)/2)]); P=real(fft(R))/NT; plot(f,P,'o')` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `NT=60; %Periodenlänge NT≤N R1=zeros(1,NT); R0=zeros(1,NT); for i=1:N for j=1:N R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+g(i)g(j)x(i)x(j); R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+g(i)g(j); end end R=zeros(1,NT); for k=1:NT R(k)=R1(k)/R0(k); end f=circshift((-fix(NT/2):fix((NT-1)/2))/(NT*dt),[0;fix((NT+1)/2)]); P=real(fft(R))/NT; plot(f,P,'o')`