Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für reelle, stochastische Einzelleistungssignale, ohne 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$ (AR1-Prozess als Beispiel, mit einem Erwartungswert verschieden von null)		`N=1000 dt=1.0 t=arange(0,N)dt phi=0.95 mxs=3 vxs=1 theta=sqrt((1-phi2)vxs) x=zeros(N) x[0]=normal(mxs,sqrt(vxs)) for i in range(1,N): x[i]=phi*(x[i-1]-mxs)+normal(mxs,theta) plot(t,x,'o') show()`	`N=1000; dt=1.0; t=(0:N-1)dt; phi=0.95; mxs=3; vxs=1; theta=sqrt((1-phi^2)vxs); x=zeros(1,N); x(1)=mxs+sqrt(vxs)randn(); for i=2:N x(i)=phi(x(i-1)-mxs)+mxs+theta*randn(); end plot(t,x,'o')`
Mittelwert	$\overline{x} = \frac{1}{N} \cdot \sum_{i = 0}^{N - 1} x_{i}$	`mean(x)`	`mean(x)`
Varianz (ohne Bessel-Korrektur, asymptotisch erwartungstreu)	$s^{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}$	`var(x)`	`var(x,1)`
Varianz (mit Bessel-Korrektur, nur für unabhängige $x_{i}$ erwartungstreu)	$s^{2} = \frac{1}{N - 1} \cdot \sum_{i = 0}^{N - 1} {(x_{i} - \overline{x})}^{2}$	`var(x,ddof=1)`	`var(x)`
lineare Trendbereinigung		`import scipy.signal as scisig x=scisig.detrend(x)`	`x=detrend(x);`
Trendbereinigung mit Polynom $p$ -ten Grades (Ergebnis der Regression $\tilde{x}$ )		`p=5 a=arange(0,N) c=polyfit(a,x,p) xtilde=zeros(N) for j in range(0,p+1): xtilde=(a*xtilde+c[j]) x-=xtilde`	`p=5; a=(1:N); c=polyfit(a,x,p); xtilde=zeros(1,N); for j=1:p+1 xtilde=(a.*xtilde+c(j)); end x=x-xtilde;`
Varianz nach Trendbereinigung mit Polynom $p$ -ten Grades (Ergebnis der Regression $\tilde{x}$ , mit Bessel-Korrektur, nur für unabhängige $x_{i}$ erwartungstreu)	$s^{2} = \frac{1}{N - p - 1} \cdot \sum_{i = 0}^{N - 1} {(x_{i} - \tilde{x})}^{2}$	`var(x,ddof=p+1)`	`N*var(x,1)/(N-p-1);`
Varianz (mit Bessel-Korrektur für korrelierte Daten, erwartungstreu)	$s^{2} = C (0)$ mit der erwartungstreuen Schätzung $C (τ_{k})$ der (Auto-)Kovarianzfunktion mit Bessel-Korrektur für korrelierte Daten (s. unten)	from numpy.fft import * from numpy.linalg import * Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N Cp=zeros(2Nc-1) mxe=mean(x) X=fft(append(x-mxe,zeros(N))) E=dt2abs(X)*2 CE=real(ifft(E))/float(dt) for k in range(-(Nc-1),Nc): Cp[k]=CE[k]/float(dt)/float(N-abs(k)) A=zeros([2Nc-1,2Nc-1]) for i in range(-(Nc-1),Nc): for j in range(-(Nc-1),Nc): A[i,j]=int(i==j)+(N-abs(j))/float(N2)-2(N-maximum(abs(i),maximum(abs(j),minimum(N,abs(i-j)))))/float(N*(N-abs(i))) Ainv=inv(A) Ce=dot(Ainv,Cp) Ce[0]	Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N Cp=zeros(1,2Nc-1); mxe=mean(x); X=fft([x-mxe,zeros(1,N)]); E=dt^2abs(X).^2; CE=real(ifft(E))/dt; for k=-(Nc-1):Nc-1 Cp(mod(k,2Nc-1)+1)=CE(mod(k,2N)+1)/dt/(N-abs(k)); end A=zeros(2Nc-1,2Nc-1); for i=-(Nc-1):Nc-1 for j=-(Nc-1):Nc-1 A(mod(i,2Nc-1)+1,mod(j,2Nc-1)+1)=(i==j)+(N-abs(j))/N^2-2(N-max(max(abs(i),abs(j)),min(N,abs(i-j))))/(N(N-abs(i))); end end Ainv=inv(A); Ce=Ainv*Cp'; Ce(1)
(Auto-)Korrelationsfunktion (erwartungstreu)	$R_{k} = R (τ_{k}) = \frac{1}{N - \| k \|} \cdot \sum_{i = 0}^{N - 1 - \| k \|} x_{i} x_{i + \| k \|}$ $τ_{k} = k Δ t$ $k = - (N - 1) \dots N - 1$ (außerhalb dieses Bereichs nicht bekannt)	direkte Berechnung: `R=zeros(2N-1) tau=zeros(2N-1) for k in range(-(N-1),N): tau[k]=kdt sum1=0 for i in range(0,N-abs(k)): sum1+=x[i]x[i+abs(k)] R[k]=sum1/float(N-abs(k)) plot(tau,R,'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * R=zeros(2N-1) tau=roll(arange(-(N-1),N)dt,N) X=fft(append(x,zeros(N))) E=dt*2abs(X)**2 RE=real(ifft(E))/float(dt) for k in range(-(N-1),N): R[k]=RE[k]/float(dt)/float(N-abs(k)) plot(tau,R,'o') show()`	direkte Berechnung: `R=zeros(1,2N-1); tau=zeros(1,2N-1); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; for i=1:N-abs(k) sum1=sum1+x(i)x(i+abs(k)); end R(mod(k,2N-1)+1)=sum1/(N-abs(k)); end plot(tau,R,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `R=zeros(1,2N-1); tau=circshift((-(N-1):N-1)dt,[0;N]); X=fft([x,zeros(1,N)]); E=dt^2abs(X).^2; RE=real(ifft(E))/dt; for k=-(N-1):N-1 R(mod(k,2N-1)+1)=RE(mod(k,2*N)+1)/dt/(N-abs(k)); end plot(tau,R,'o')`
(Auto-)Kovarianzfunktion (asymptotisch erwartungstreu)	$C_{k} = C (τ_{k}) = \frac{1}{N - \| k \|} \cdot \sum_{i = 0}^{N - 1 - \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})$ $τ_{k} = k Δ t$ $k = - (N - 1) \dots N - 1$ (außerhalb dieses Bereichs nicht bekannt)	direkte Berechnung: `C=zeros(2N-1) tau=zeros(2N-1) mxe=mean(x) for k in range(-(N-1),N): tau[k]=kdt sum1=0 for i in range(0,N-abs(k)): sum1+=(x[i]-mxe)(x[i+abs(k)]-mxe) C[k]=sum1/float(N-abs(k)) plot(tau,C,'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * C=zeros(2N-1) tau=roll(arange(-(N-1),N)dt,N) mxe=mean(x) X=fft(append(x-mxe,zeros(N))) E=dt*2abs(X)**2 CE=real(ifft(E))/float(dt) for k in range(-(N-1),N): C[k]=CE[k]/float(dt)/float(N-abs(k)) plot(tau,C,'o') show()`	direkte Berechnung: `C=zeros(1,2N-1); tau=zeros(1,2N-1); mxe=mean(x); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; for i=1:N-abs(k) sum1=sum1+(x(i)-mxe)(x(i+abs(k))-mxe); end C(mod(k,2N-1)+1)=sum1/(N-abs(k)); end plot(tau,C,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `C=zeros(1,2N-1); tau=circshift((-(N-1):N-1)dt,[0;N]); mxe=mean(x); X=fft([x-mxe,zeros(1,N)]); E=dt^2abs(X).^2; CE=real(ifft(E))/dt; for k=-(N-1):N-1 C(mod(k,2N-1)+1)=CE(mod(k,2*N)+1)/dt/(N-abs(k)); end plot(tau,C,'o')`
(Auto-)Kovarianzfunktion (mit Bessel-Korrektur für korrelierte Daten, erwartungstreu)	s. H. Nobach: Practical Realization of Bessel's Correction for a Bias-Free Estimation of the Auto-Covariance and the Cross-Covariance Functions	direkte Berechnung: from numpy.linalg import * Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N Cp=zeros(2Nc-1) tau=zeros(2Nc-1) mxe=mean(x) for k in range(-(Nc-1),Nc): tau[k]=kdt sum1=0 for i in range(0,N-abs(k)): sum1+=(x[i]-mxe)(x[i+abs(k)]-mxe) Cp[k]=sum1/float(N-abs(k)) A=zeros([2Nc-1,2Nc-1]) for i in range(-(Nc-1),Nc): for j in range(-(Nc-1),Nc): A[i,j]=int(i==j)+(N-abs(j))/float(N*2)-2(N-maximum(abs(i),maximum(abs(j),minimum(N,abs(i-j)))))/float(N(N-abs(i))) Ainv=inv(A) C=dot(Ainv,Cp) plot(tau,C,'o') show() Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: from numpy.fft import from numpy.linalg import * Nc=200 #bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N Cp=zeros(2Nc-1) tau=roll(arange(-(Nc-1),Nc)dt,Nc) mxe=mean(x) X=fft(append(x-mxe,zeros(N))) E=dt*2abs(X)*2 CE=real(ifft(E))/float(dt) for k in range(-(Nc-1),Nc): Cp[k]=CE[k]/float(dt)/float(N-abs(k)) A=zeros([2Nc-1,2Nc-1]) for i in range(-(Nc-1),Nc): for j in range(-(Nc-1),Nc): A[i,j]=int(i==j)+(N-abs(j))/float(N2)-2(N-maximum(abs(i),maximum(abs(j),minimum(N,abs(i-j)))))/float(N*(N-abs(i))) Ainv=inv(A) C=dot(Ainv,Cp) plot(tau,C,'o') show()	direkte Berechnung: Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N Cp=zeros(1,2Nc-1); tau=zeros(1,2Nc-1); mxe=mean(x); for k=-(Nc-1):Nc-1 tau(mod(k,2Nc-1)+1)=kdt; sum1=0; for i=1:N-abs(k) sum1=sum1+(x(i)-mxe)(x(i+abs(k))-mxe); end Cp(mod(k,2Nc-1)+1)=sum1/(N-abs(k)); end A=zeros(2Nc-1,2Nc-1); for i=-(Nc-1):Nc-1 for j=-(Nc-1):Nc-1 A(mod(i,2Nc-1)+1,mod(j,2Nc-1)+1)=(i==j)+(N-abs(j))/N^2-2(N-max(max(abs(i),abs(j)),min(N,abs(i-j))))/(N(N-abs(i))); end end Ainv=inv(A); C=AinvCp'; plot(tau,C,'o') Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N Cp=zeros(1,2Nc-1); tau=circshift((-(Nc-1):Nc-1)dt,[0;Nc]); mxe=mean(x); X=fft([x-mxe,zeros(1,N)]); E=dt^2abs(X).^2; CE=real(ifft(E))/dt; for k=-(Nc-1):Nc-1 Cp(mod(k,2Nc-1)+1)=CE(mod(k,2N)+1)/dt/(N-abs(k)); end A=zeros(2Nc-1,2Nc-1); for i=-(Nc-1):Nc-1 for j=-(Nc-1):Nc-1 A(mod(i,2Nc-1)+1,mod(j,2Nc-1)+1)=(i==j)+(N-abs(j))/N^2-2(N-max(max(abs(i),abs(j)),min(N,abs(i-j))))/(N(N-abs(i))); end end Ainv=inv(A); C=Ainv*Cp'; plot(tau,C,'o')
(Auto-)Kovarianzfunktion (mit lokaler Normierung, asymptotisch erwartungstreu)	$C_{k} = C (τ_{k}) = \frac{s^{2} \cdot \sum_{i = 0}^{N - 1 - \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})}{\sqrt{[\sum_{i = 0}^{N - 1 - \| k \|} {(x_{i} - \overline{x})}^{2}] [\sum_{i = \| k \|}^{N - 1} {(x_{i} - \overline{x})}^{2}]}}$ $τ_{k} = k Δ t$ $k = - (N - 1) \dots N - 1$ (außerhalb dieses Bereichs nicht bekannt)	direkte Berechnung: `C=zeros(2N-1) tau=zeros(2N-1) mxe=mean(x) vxe=var(x) for k in range(-(N-1),N): tau[k]=kdt sum1=0 sum2=0 sum3=0 for i in range(0,N-abs(k)): sum1+=(x[i]-mxe)(x[i+abs(k)]-mxe) sum2+=(x[i]-mxe)2 sum3+=(x[i+abs(k)]-mxe)2 C[k]=vxesum1/sqrt(sum2sum3) plot(tau,C,'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * tau=roll(arange(-(N-1),N)dt,N) mxe=mean(x) vxe=var(x) X=fft(append(x-mxe,zeros(N))) Q=fft(append((x-mxe)2,zeros(N))) O=fft(append(ones(N),zeros(N))) EX=dt2abs(X)2 EA=dt2conj(Q)O EB=dt*2conj(O)Q CEX=real(ifft(EX))/float(dt) CEA=real(ifft(EA))/float(dt) CEB=real(ifft(EB))/float(dt) C=vxeCEX/sqrt(CEACEB) C=append(C[0:N],C[N+1:2N]) plot(tau,C,'o') show()`	direkte Berechnung: `C=zeros(1,2N-1); tau=zeros(1,2N-1); mxe=mean(x); vxe=var(x,1); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; sum2=0; sum3=0; for i=1:N-abs(k) sum1=sum1+(x(i)-mxe)(x(i+abs(k))-mxe); sum2=sum2+(x(i)-mxe)^2; sum3=sum3+(x(i+abs(k))-mxe)^2; end C(mod(k,2N-1)+1)=vxesum1/sqrt(sum2sum3); end plot(tau,C,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `tau=circshift((-(N-1):N-1)dt,[0;N]); mxe=mean(x); vxe=var(x,1); X=fft([x-mxe,zeros(1,N)]); Q=fft([(x-mxe).^2,zeros(1,N)]); O=fft([ones(1,N),zeros(1,N)]); EX=dt^2abs(X).^2; EA=dt^2conj(Q).O; EB=dt^2conj(O).Q; CEX=real(ifft(EX))/dt; CEA=real(ifft(EA))/dt; CEB=real(ifft(EB))/dt; C=vxeCEX./sqrt(CEA.CEB); C=[C(1:N),C(N+2:end)]; plot(tau,C,'o')`
(Auto-)Korrelationskoeffizientenfunktion (asymptotisch erwartungstreu)	$ρ_{k} = ρ (τ_{k}) = \frac{1}{s^{2} (N - \| k \|)} \cdot \sum_{i = 0}^{N - 1 - \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})$ $τ_{k} = k Δ t$ $k = - (N - 1) \dots N - 1$ (außerhalb dieses Bereichs nicht bekannt)	direkte Berechnung: `rho=zeros(2N-1) tau=zeros(2N-1) mxe=mean(x) vxe=var(x) for k in range(-(N-1),N): tau[k]=kdt sum1=0 for i in range(0,N-abs(k)): sum1+=(x[i]-mxe)(x[i+abs(k)]-mxe) rho[k]=sum1/float(N-abs(k))/vxe plot(tau,rho,'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import * rho=zeros(2N-1) tau=roll(arange(-(N-1),N)dt,N) mxe=mean(x) vxe=var(x) X=fft(append(x-mxe,zeros(N))) E=dt*2abs(X)**2 CE=real(ifft(E))/float(dt) for k in range(-(N-1),N): rho[k]=CE[k]/float(dt)/float(N-abs(k))/vxe plot(tau,rho,'o') show()`	direkte Berechnung: `rho=zeros(1,2N-1); tau=zeros(1,2N-1); mxe=mean(x); vxe=var(x,1); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; for i=1:N-abs(k) sum1=sum1+(x(i)-mxe)(x(i+abs(k))-mxe); end rho(mod(k,2N-1)+1)=sum1/(N-abs(k))/vxe; end plot(tau,rho,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `rho=zeros(1,2N-1); tau=circshift((-(N-1):N-1)dt,[0;N]); mxe=mean(x); vxe=var(x,1); X=fft([x-mxe,zeros(1,N)]); E=dt^2abs(X).^2; CE=real(ifft(E))/dt; for k=-(N-1):N-1 rho(mod(k,2N-1)+1)=CE(mod(k,2*N)+1)/dt/(N-abs(k))/vxe; end plot(tau,rho,'o')`
(Auto-)Korrelationskoeffizientenfunktion (mit lokaler Normierung, asymptotisch erwartungstreu)	$ρ_{k} = ρ (τ_{k}) = \frac{\sum_{i = 0}^{N - 1 - \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})}{\sqrt{[\sum_{i = 0}^{N - 1 - \| k \|} {(x_{i} - \overline{x})}^{2}] [\sum_{i = \| k \|}^{N - 1} {(x_{i} - \overline{x})}^{2}]}}$ $τ_{k} = k Δ t$ $k = - (N - 1) \dots N - 1$ (außerhalb dieses Bereichs nicht bekannt)	direkte Berechnung: `rho=zeros(2N-1) tau=zeros(2N-1) mxe=mean(x) for k in range(-(N-1),N): tau[k]=kdt sum1=0 sum2=0 sum3=0 for i in range(0,N-abs(k)): sum1+=(x[i]-mxe)(x[i+abs(k)]-mxe) sum2+=(x[i]-mxe)2 sum3+=(x[i+abs(k)]-mxe)2 rho[k]=sum1/sqrt(sum2sum3) plot(tau,rho,'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=roll(arange(-(N-1),N)dt,N) mxe=mean(x) X=fft(append(x-mxe,zeros(N))) Q=fft(append((x-mxe)2,zeros(N))) O=fft(append(ones(N),zeros(N))) EX=dt2abs(X)2 EA=dt2conj(Q)O EB=dt*2conj(O)Q CEX=real(ifft(EX))/float(dt) CEA=real(ifft(EA))/float(dt) CEB=real(ifft(EB))/float(dt) rho=CEX/sqrt(CEACEB) rho=append(rho[0:N],rho[N+1:2*N]) plot(tau,rho,'o') show()`	direkte Berechnung: `rho=zeros(1,2N-1); tau=zeros(1,2N-1); mxe=mean(x); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; sum2=0; sum3=0; for i=1:N-abs(k) sum1=sum1+(x(i)-mxe)(x(i+abs(k))-mxe); sum2=sum2+(x(i)-mxe)^2; sum3=sum3+(x(i+abs(k))-mxe)^2; end rho(mod(k,2N-1)+1)=sum1/sqrt(sum2sum3); end plot(tau,rho,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `tau=circshift((-(N-1):N-1)dt,[0;N]); mxe=mean(x); X=fft([x-mxe,zeros(1,N)]); Q=fft([(x-mxe).^2,zeros(1,N)]); O=fft([ones(1,N),zeros(1,N)]); EX=dt^2abs(X).^2; EA=dt^2conj(Q).O; EB=dt^2conj(O).Q; CEX=real(ifft(EX))/dt; CEA=real(ifft(EA))/dt; CEB=real(ifft(EB))/dt; rho=CEX./sqrt(CEA.CEB); rho=[rho(1:N),rho(N+2:end)]; plot(tau,rho,'o')`
(Auto-)Leistungsdichtespektrum (erwartungstreu)	(wegen der nötigen Normierung nur über die Korrelationsfunktion, evtl. Einschränkung des verwendeten Bereichs der Korrelationsfunktion zur Reduktion der Schätzvarianz des Leistungsdichtespektrums, $K \leq 2 N - 1$ ) $S_{j} = S (f_{j}) = Δ t \cdot F F T {R_{k}} = Δ t \cdot \sum_{k = - ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{k} 𝐞^{- 2 π 𝐢 j k / K}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung des primären Spektrums: `from numpy.fft import * X=fft(append(x,zeros(N))) E=dt*2abs(X)*2 RE=real(ifft(E))/float(dt) K=100 K=minimum(2N-1,K) R=zeros(K) for k in range(-(K//2),(K-1)//2+1): R[k]=RE[k]/float(dt)/float(N-abs(k)) f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) plot(f,S,'o') show()` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import * K=100 K=minimum(2N-1,K) R=zeros(K) for k in range(-(K//2),(K-1)//2+1): sum1=0 for i in range(0,N-abs(k)): sum1+=x[i]x[i+abs(k)] R[k]=sum1/float(N-abs(k)) f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) plot(f,S,'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung des primären Spektrums: `X=fft([x,zeros(1,N)]); E=dt^2abs(X).^2; RE=real(ifft(E))/dt; K=100; K=min(2N-1,K); R=zeros(1,K); for k=-floor(K/2):floor((K-1)/2) R(mod(k,K)+1)=RE(mod(k,2N)+1)/dt/(N-abs(k)); end f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); plot(f,S,'o')` Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `K=100; K=min(2N-1,K); R=zeros(1,K); for k=-floor(K/2):floor((K-1)/2) sum1=0; for i=1:N-abs(k) sum1=sum1+x(i)x(i+abs(k)); end R(mod(k,K)+1)=sum1/(N-abs(k)); end f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dt*real(fft(R)); plot(f,S,'o')`