Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für reelle, stochastische Einzelleistungssignale, 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$ (AR1-Prozess als Beispiel, mit einem Erwartungswert verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)		`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)) while random()>1/(1+exp(-(x[0]-mxs))): x[0]=normal(mxs,sqrt(vxs)) for i in range(1,N): x[i]=phi(x[i-1]-mxs)+normal(mxs,theta) while random()>1/(1+exp(-(x[i]-mxs))): 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(); while rand>1/(1+exp(-(x(1)-mxs))) x(1)=mxs+sqrt(vxs)randn(); end for i=2:N x(i)=phi(x(i-1)-mxs)+mxs+thetarandn(); while rand>1/(1+exp(-(x(i)-mxs))) x(i)=phi(x(i-1)-mxs)+mxs+thetarandn(); end end 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$ (Inverse der verschobenen Sigmoidfunktion)		`g=zeros(N) for i in range(0,N): g[i]=(1+exp(-(x[i]-mxs)))**(1+phi) plot(t,g,'o') show()`	`g=zeros(1,N); for i=1:N g(i)=(1+exp(-(x(i)-mxs)))^(1+phi); end plot(t,g,'o')`
Mittelwert	$\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 Bessel-Korrektur, asymptotisch erwartungstreu)	$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`
Varianz (mit Bessel-Korrektur für korrelierte Daten, erwartungstreu)	$s^{2} = \frac{\sum_{i = 0}^{N - 1} g_{i} {(x_{i} - \overline{x})}^{2}}{\sum_{i = 0}^{N - 1} g_{i}} + s_{\overline{x}}^{2} = \frac{\sum_{i = 0}^{N - 1} g_{i} x_{i}^{2}}{\sum_{i = 0}^{N - 1} g_{i}} - {\overline{x}}^{2} + s_{\overline{x}}^{2}$ mit der erwartungstreuen Schätzung $s_{\overline{x}}^{2}$ der Schätzvarianz des Mittelwertschätzers für korrelierte Daten aus der erwartungstreuen Schätzung 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 mxe=sum(gx)/float(sum(g)) xe=append(g(x-mxe),zeros(N)) ge=append(g,zeros(N)) X=real(ifft(abs(fft(xe))2)) Y=real(ifft(abs(fft(ge))2)) Cp=zeros(2Nc-1) for k in range(-(Nc-1),Nc): Cp[k]=X[k]/Y[k] D=sum(g) A=zeros([2Nc-1,2Nc-1]) for k in range(-(Nc-1),Nc): Gk=real(ifft(conj(fft(geroll(ge,-k)))fft(ge))) Hk=real(ifft(conj(fft(ge))fft(geroll(ge,k)))) for j in range(-(Nc-1),Nc): A[k,j]=int(k==j)-(Gk[j]+Hk[j])/float(Y[k]D)+Y[j]/float(D*2) Ainv=inv(A) Ce=dot(Ainv,Cp) sC=0 for k in range(-(Nc-1),Nc): sC+=(N-abs(k))Ce[k] sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))2+sC/float(N2)	Nc=200; %bitte so waehlen, dass die Korrelationsfunktion ausreichend abgeklungen ist, aber deutlich kleiner als N mxe=sum(g.x)/sum(g); xe=[g.(x-mxe),zeros(1,N)]; ge=[g,zeros(1,N)]; X=real(ifft(abs(fft(xe)).^2)); Y=real(ifft(abs(fft(ge)).^2)); Cp=zeros(1,2Nc-1); for k=-(Nc-1):Nc-1 Cp(mod(k,2Nc-1)+1)=X(mod(k,2N)+1)/Y(mod(k,2N)+1); end D=sum(g); A=zeros(2Nc-1,2Nc-1); for k=-(Nc-1):Nc-1 Gk=real(ifft(conj(fft(ge.circshift(ge,[0;-k]))).fft(ge))); Hk=real(ifft(conj(fft(ge)).fft(ge.circshift(ge,[0;k])))); for j=-(Nc-1):Nc-1 A(mod(k,2Nc-1)+1,mod(j,2Nc-1)+1)=(k==j)-(Gk(mod(j,2N)+1)+Hk(mod(j,2N)+1))/(Y(mod(k,2N)+1)D)+Y(mod(j,2N)+1)/(D^2); end end Ainv=inv(A); Ce=AinvCp'; sC=0; for k=-(Nc-1):Nc-1 sC=sC+(N-abs(k))Ce(mod(k,2Nc-1)+1); end sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2+sC/(N^2)
(Auto-)Korrelationsfunktion (erwartungstreu)	$R_{k} = R (τ_{k}) = \frac{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} x_{i} x_{i + \| k \|}}{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{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 sum0=0 for i in range(0,N-abs(k)): sum1+=g[i]g[i+abs(k)]x[i]x[i+abs(k)] sum0+=g[i]g[i+abs(k)] R[k]=sum1/float(sum0) plot(tau,R,'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=roll(arange(-(N-1),N)dt,N) X1=fft(append(gx,zeros(N))) E1=dt*2abs(X1)2 X0=fft(append(g,zeros(N))) E0=dt2abs(X0)2 R=real(ifft(E1))/real(ifft(E0)) R=append(R[0:N],R[N+1:2N]) 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; sum0=0; for i=1:N-abs(k) sum1=sum1+g(i)g(i+abs(k))x(i)x(i+abs(k)); sum0=sum0+g(i)g(i+abs(k)); end R(mod(k,2N-1)+1)=sum1/sum0; end plot(tau,R,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `tau=circshift((-(N-1):N-1)dt,[0;N]); X1=fft([g.x,zeros(1,N)]); E1=dt^2abs(X1).^2; X0=fft([g,zeros(1,N)]); E0=dt^2*abs(X0).^2; R=real(ifft(E1))./real(ifft(E0)); R=[R(1:N),R(N+2:end)]; plot(tau,R,'o')`
(Auto-)Kovarianzfunktion (asymptotisch erwartungstreu)	$C_{k} = C (τ_{k}) = \frac{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})}{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|}}$ $τ_{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=sum(gx)/float(sum(g)) for k in range(-(N-1),N): tau[k]=kdt sum1=0 sum0=0 for i in range(0,N-abs(k)): sum1+=g[i]g[i+abs(k)](x[i]-mxe)(x[i+abs(k)]-mxe) sum0+=g[i]g[i+abs(k)] C[k]=sum1/float(sum0) 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=sum(gx)/float(sum(g)) X1=fft(append(g(x-mxe),zeros(N))) E1=dt2abs(X1)2 X0=fft(append(g,zeros(N))) E0=dt2abs(X0)2 C=real(ifft(E1))/real(ifft(E0)) 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=sum(g.x)/sum(g); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; sum0=0; for i=1:N-abs(k) sum1=sum1+g(i)g(i+abs(k))(x(i)-mxe)(x(i+abs(k))-mxe); sum0=sum0+g(i)g(i+abs(k)); end C(mod(k,2N-1)+1)=sum1/sum0; end plot(tau,C,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `tau=circshift((-(N-1):N-1)dt,[0;N]); mxe=sum(g.x)/sum(g); X1=fft([g.(x-mxe),zeros(1,N)]); E1=dt^2abs(X1).^2; X0=fft([g,zeros(1,N)]); E0=dt^2*abs(X0).^2; C=real(ifft(E1))./real(ifft(E0)); C=[C(1:N),C(N+2: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) S=zeros(2Nc-1) tau=zeros(2Nc-1) mxe=sum(gx)/float(sum(g)) for k in range(-(Nc-1),Nc): tau[k]=kdt sum1=0 sum0=0 for i in range(0,N-abs(k)): sum1+=g[i]g[i+abs(k)](x[i]-mxe)(x[i+abs(k)]-mxe) sum0+=g[i]g[i+abs(k)] Cp[k]=sum1/float(sum0) S[k]=sum0 D=sum(g) A=zeros([2Nc-1,2Nc-1]) for k in range(-(Nc-1),Nc): for j in range(-(Nc-1),Nc): s=0 for i in range(maximum(0,maximum(-j,-k)),minimum(N,minimum(N-j,N-k))): s+=g[i]g[i+j]g[i+k] for i in range(maximum(0,maximum(-j,k-j)),minimum(N,minimum(N-j,N+k-j))): s+=g[i]g[i+j]g[i+j-k] A[k,j]=int(k==j)-s/float(S[k]D)+S[j]/float(D*2) 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 tau=roll(arange(-(Nc-1),Nc)dt,Nc) mxe=sum(gx)/float(sum(g)) xe=append(g(x-mxe),zeros(N)) ge=append(g,zeros(N)) X=real(ifft(abs(fft(xe))2)) Y=real(ifft(abs(fft(ge))2)) Cp=zeros(2Nc-1) for k in range(-(Nc-1),Nc): Cp[k]=X[k]/Y[k] D=sum(g) A=zeros([2Nc-1,2Nc-1]) for k in range(-(Nc-1),Nc): Gk=real(ifft(conj(fft(geroll(ge,-k)))fft(ge))) Hk=real(ifft(conj(fft(ge))fft(geroll(ge,k)))) for j in range(-(Nc-1),Nc): A[k,j]=int(k==j)-(Gk[j]+Hk[j])/float(Y[k]D)+Y[j]/float(D*2) 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); S=zeros(1,2Nc-1); tau=zeros(1,2Nc-1); mxe=sum(g.x)/sum(g); for k=-(Nc-1):Nc-1 tau(mod(k,2Nc-1)+1)=kdt; sum1=0; sum0=0; for i=1:N-abs(k) sum1=sum1+g(i)g(i+abs(k))(x(i)-mxe)(x(i+abs(k))-mxe); sum0=sum0+g(i)g(i+abs(k)); end Cp(mod(k,2Nc-1)+1)=sum1/sum0; S(mod(k,2Nc-1)+1)=sum0; end D=sum(g); A=zeros(2Nc-1,2Nc-1); for k=-(Nc-1):Nc-1 for j=-(Nc-1):Nc-1 s=0; for i=max(0,max(-j,-k))+1:min(N,min(N-j,N-k)) s=s+g(i)g(i+j)g(i+k); end for i=max(0,max(-j,k-j))+1:min(N,min(N-j,N+k-j)) s=s+g(i)g(i+j)g(i+j-k); end A(mod(k,2Nc-1)+1,mod(j,2Nc-1)+1)=(k==j)-s/(S(mod(k,2Nc-1)+1)D)+S(mod(j,2Nc-1)+1)/(D^2); 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 tau=circshift((-(Nc-1):Nc-1)dt,[0;Nc]); mxe=sum(g.x)/sum(g); xe=[g.(x-mxe),zeros(1,N)]; ge=[g,zeros(1,N)]; X=real(ifft(abs(fft(xe)).^2)); Y=real(ifft(abs(fft(ge)).^2)); Cp=zeros(1,2Nc-1); for k=-(Nc-1):Nc-1 Cp(mod(k,2Nc-1)+1)=X(mod(k,2N)+1)/Y(mod(k,2N)+1); end D=sum(g); A=zeros(2Nc-1,2Nc-1); for k=-(Nc-1):Nc-1 Gk=real(ifft(conj(fft(ge.circshift(ge,[0;-k]))).fft(ge))); Hk=real(ifft(conj(fft(ge)).fft(ge.circshift(ge,[0;k])))); for j=-(Nc-1):Nc-1 A(mod(k,2Nc-1)+1,mod(j,2Nc-1)+1)=(k==j)-(Gk(mod(j,2N)+1)+Hk(mod(j,2N)+1))/(Y(mod(k,2N)+1)D)+Y(mod(j,2N)+1)/(D^2); 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 \|} g_{i} g_{i + \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})}{\sqrt{[\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} {(x_{i} - \overline{x})}^{2}] [\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} {(x_{i + \| k \|} - \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=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 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+=g[i]g[i+abs(k)](x[i]-mxe)(x[i+abs(k)]-mxe) sum2+=g[i]g[i+abs(k)](x[i]-mxe)2 sum3+=g[i]g[i+abs(k)](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=sum(gx)/float(sum(g)) vxe=sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 X=fft(append(g(x-mxe),zeros(N))) Q=fft(append(g(x-mxe)2,zeros(N))) O=fft(append(g,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=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; 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+g(i)g(i+abs(k))(x(i)-mxe)(x(i+abs(k))-mxe); sum2=sum2+g(i)g(i+abs(k))(x(i)-mxe)^2; sum3=sum3+g(i)g(i+abs(k))(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=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; X=fft([g.(x-mxe),zeros(1,N)]); Q=fft([g.(x-mxe).^2,zeros(1,N)]); O=fft([g,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{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})}{s^{2} \cdot \sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|}}$ $τ_{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=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 for k in range(-(N-1),N): tau[k]=kdt sum1=0 sum0=0 for i in range(0,N-abs(k)): sum1+=g[i]g[i+abs(k)](x[i]-mxe)(x[i+abs(k)]-mxe) sum0+=g[i]g[i+abs(k)] rho[k]=sum1/sum0/vxe 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=sum(gx)/float(sum(g)) vxe=sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 X1=fft(append(g(x-mxe),zeros(N))) E1=dt*2abs(X1)2 X0=fft(append(g,zeros(N))) E0=dt2abs(X0)2 rho=real(ifft(E1))/real(ifft(E0))/vxe rho=append(rho[0:N],rho[N+1:2N]) plot(tau,rho,'o') show()`	direkte Berechnung: `rho=zeros(1,2N-1); tau=zeros(1,2N-1); mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; sum0=0; for i=1:N-abs(k) sum1=sum1+g(i)g(i+abs(k))(x(i)-mxe)(x(i+abs(k))-mxe); sum0=sum0+g(i)g(i+abs(k)); end rho(mod(k,2N-1)+1)=sum1/sum0/vxe; end plot(tau,rho,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `tau=circshift((-(N-1):N-1)dt,[0;N]); mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; X1=fft([g.(x-mxe),zeros(1,N)]); E1=dt^2abs(X1).^2; X0=fft([g,zeros(1,N)]); E0=dt^2*abs(X0).^2; rho=real(ifft(E1))./real(ifft(E0))/vxe; rho=[rho(1:N),rho(N+2:end)]; plot(tau,rho,'o')`
(Auto-)Korrelationskoeffizientenfunktion (mit lokaler Normierung, asymptotisch erwartungstreu)	$ρ_{k} = ρ (τ_{k}) = \frac{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} (x_{i} - \overline{x}) (x_{i + \| k \|} - \overline{x})}{\sqrt{[\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} {(x_{i} - \overline{x})}^{2}] [\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} {(x_{i + \| k \|} - \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=sum(gx)/float(sum(g)) 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+=g[i]g[i+abs(k)](x[i]-mxe)(x[i+abs(k)]-mxe) sum2+=g[i]g[i+abs(k)](x[i]-mxe)2 sum3+=g[i]g[i+abs(k)](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=sum(gx)/float(sum(g)) X=fft(append(g(x-mxe),zeros(N))) Q=fft(append(g(x-mxe)2,zeros(N))) O=fft(append(g,zeros(N))) EX=dt2abs(X)2 EA=dt2conj(Q)O EB=dt2conj(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=sum(g.x)/sum(g); 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+g(i)g(i+abs(k))(x(i)-mxe)(x(i+abs(k))-mxe); sum2=sum2+g(i)g(i+abs(k))(x(i)-mxe)^2; sum3=sum3+g(i)g(i+abs(k))(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=sum(g.x)/sum(g); X=fft([g.(x-mxe),zeros(1,N)]); Q=fft([g.(x-mxe).^2,zeros(1,N)]); O=fft([g,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=REX./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 der primären Spektren: `from numpy.fft import * X1=fft(append(gx,zeros(N))) E1=dt2abs(X1)2 X0=fft(append(g,zeros(N))) E0=dt2abs(X0)2 RA=real(ifft(E1))/real(ifft(E0)) K=100 K=minimum(2N-1,K) R=zeros(K) for k in range(-(K//2),(K-1)//2+1): R[k]=RA[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 sum0=0 for i in range(0,N-abs(k)): sum1+=g[i]g[i+abs(k)]x[i]x[i+abs(k)] sum0+=g[i]g[i+abs(k)] R[k]=sum1/float(sum0) f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dt*real(fft(R)) plot(f,S,'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `X1=fft([g.x,zeros(1,N)]); E1=dt^2abs(X1).^2; X0=fft([g,zeros(1,N)]); E0=dt^2abs(X0).^2; RA=real(ifft(E1))./real(ifft(E0)); 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)=RA(mod(k,2N)+1); 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; sum0=0; for i=1:N-abs(k) sum1=sum1+g(i)g(i+abs(k))x(i)x(i+abs(k)); sum0=sum0+g(i)g(i+abs(k)); end R(mod(k,K)+1)=sum1/sum0; end f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); plot(f,S,'o')`