Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für stochastisch abgetastete, komplexe Einzelleistungssignale, ohne Gewichtung

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von Wertepaaren $(t_{i}, x_{i})$ mit den Messzeitpunkten $t_{i}$ und den komplexen Messwerten $x_{i}$ (zwei AR1-Prozess als Beispiel, mit Erwartungswert verschieden von null)		`T=10000.0 Ti=10.0 dr=1.0 mxs=3.0+1.0j vxs=1.0 t=[] x=[] te=0.0 xe=normal(real(mxs),sqrt(vxs/2.0))+1jnormal(imag(mxs),sqrt(vxs/2.0)) while te<T: tp=exponential(1.0/dr) te+=tp if te<T: t.append(te) phi=exp(-tp/Ti) theta=sqrt((1-phi2)vxs/2.0) xe=(xe-mxs)phi+normal(real(mxs),theta)+1jnormal(imag(mxs),theta) x.append(xe) t=array(t) x=array(x) plot(t,real(x),'o',t,imag(x),'o') show()`	`T=10000; Ti=10; dr=1; mxs=3+1j; vxs=1; t=[]; x=[]; te=0; xe=real(mxs)+sqrt(vxs/2)randn()+1j(imag(mxs)+sqrt(vxs/2)randn()); while te<T tp=-log(1-rand())/dr; te=te+tp; if te<T t(end+1)=te; phi=exp(-tp/Ti); theta=sqrt((1-phi^2)vxs/2); xe=(xe-mxs)phi+real(mxs)+thetarandn()+1j(imag(mxs)+thetarandn()); x(end+1)=xe; end end plot(t,real(x),'o',t,imag(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})}^{} (x_{i} - \overline{x}) = \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}^{} x_{i}) - {\overline{x}}^{*} \overline{x} = (\frac{1}{N} \cdot \sum_{i = 0}^{N - 1} {\| x_{i} \|}^{2}) - {\| \overline{x} \|}^{2}$	`var(x)`	`var(x,1)`
Autokorrelationsfunktion und Leistungsdichtespektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ )	$R_{k} = R (τ_{k}) = \frac{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) x_{i}^{*} x_{j}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i})}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K)+0j R0=zeros(K) i=1 while i<N: j=i-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/float(dt))) R1[k]+=conj(x[i])x[j] R0[k]+=1.0 j-=1 i+=1 for k in range(1,(K+1)//2): R1[k]=conj(R1[-k]) R0[k]=R0[-k] R=R1/R0 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K)+0j; R0=zeros(1,K); i=2; while i<=N j=i-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+conj(x(i))x(j); R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+1; j=j-1; end i=i+1; end for k=2:fix((K+1)/2) R1(k)=conj(R1(K-k+2)); R0(k)=R0(K-k+2); end R=R1./R0; plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ , mit lokaler Normierung)	$R_{k} = R (τ_{k}) = \frac{s^{2} [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) {(x_{i} - \overline{x})}^{*} (x_{j} - \overline{x})]}{\sqrt{[\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) {\| x_{i} - \overline{x} \|}^{2}] [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) {\| x_{j} - \overline{x} \|}^{2}]}} + {\| \overline{x} \|}^{2}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) mxe=mean(x) vxe=var(x) R1=zeros(K)+0j R2=zeros(K) R3=zeros(K) i=1 while i<N: j=i-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/float(dt))) R1[k]+=conj(x[i]-mxe)(x[j]-mxe) R2[k]+=abs(x[i]-mxe)2 R3[k]+=abs(x[j]-mxe)2 j-=1 i+=1 for k in range(1,(K+1)//2): R1[k]=conj(R1[-k]) R2[k]=R3[-k] R3[k]=R2[-k] R=vxeR1/sqrt(R2R3)+abs(mxe)*2 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); mxe=mean(x); vxe=var(x,1); R1=zeros(1,K)+0j; R2=zeros(1,K); R3=zeros(1,K); i=2; while i<=N j=i-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+conj(x(i)-mxe)(x(j)-mxe); R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+abs(x(i)-mxe)^2; R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+abs(x(j)-mxe)^2; j=j-1; end i=i+1; end for k=2:fix((K+1)/2) R1(k)=conj(R1(K-k+2)); R2(k)=R3(K-k+2); R3(k)=R2(K-k+2); end R=vxeR1./sqrt(R2.R3)+abs(mxe)^2; plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, ohne Normierung)	$E_{j} = E (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} {\| x_{i} \|}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, k} = R_{E} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{k} = R (τ_{k}) = \frac{R_{E, k}}{T - \| τ_{k} \|}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j for i in range(0,N): X+=x[i]exp(-2jpifpt[i]) E=T*2(abs(X)2-sum(abs(x)2))/float(N*2-N) RE=ifft(E)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T-abs(tau[k])) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; for i=1:N X=X+x(i)exp(-2jpifpt(i)); end E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); RE=ifft(E)/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/(T-abs(tau(mod(k+K,K)+1))); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit Normierung)	$E_{j} = E (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} {\| x_{i} \|}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j}^{'} \|}^{2} - N}$ $X_{j} = \sum_{i = 0}^{N - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, k} = R_{E} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, k}^{'} = R_{E}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{k} = R (τ_{k}) = \frac{R_{E, k}}{R_{E, k}^{'}}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j for i in range(0,N): X+=x[i]exp(-2jpifpt[i]) Xp+=exp(-2jpifpt[i]) E=T2(abs(X)2-sum(abs(x)2))/float(N2-N) Ep=T2(abs(Xp)2-N)/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; for i=1:N X=X+x(i)exp(-2jpifpt(i)); Xp=Xp+exp(-2jpifpt(i)); end E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); Ep=T^2(abs(Xp).^2-N)/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/REp(mod(k+length(REp),length(REp))+1); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit lokaler Normierung)	$E_{j} = E (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = \frac{T^{2}}{N (N - 1)} {X_{j}^{'' } X_{j}^{'} - \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2}}$ $E_{j}^{''} = E^{''} (f_{j}) = \frac{T^{2}}{N (N - 1)} {X_{j}^{' } X_{j}^{''} - \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} (x_{i} - \overline{x}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{''} = \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, k} = R_{E} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, k}^{'} = R_{E}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, k}^{''} = R_{E}^{''} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}^{''}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j}^{''} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{k} = R (τ_{k}) = \frac{s^{2} R_{E, k}}{\sqrt{R_{E, k}^{'} R_{E, k}^{''}}} + {\| \overline{x} \|}^{2}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 mxe=mean(x) vxe=var(x) J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j Xpp=zeros(size(fp))+0j for i in range(0,N): X+=(x[i]-mxe)exp(-2jpifpt[i]) Xp+=exp(-2jpifpt[i]) Xpp+=abs(x[i]-mxe)2exp(-2jpifpt[i]) E=T2(abs(X)2-sum(abs(x-mxe)2))/float(N2-N) Ep=T2(conj(Xpp)Xp-sum(abs(x-mxe)2))/float(N2-N) Epp=T*2(conj(Xp)Xpp-sum(abs(x-mxe)2))/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+abs(mxe)2 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; mxe=mean(x); vxe=var(x,1); J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; Xpp=zeros(size(fp))+0j; for i=1:N X=X+(x(i)-mxe)exp(-2jpifpt(i)); Xp=Xp+exp(-2jpifpt(i)); Xpp=Xpp+abs(x(i)-mxe)^2exp(-2jpifpt(i)); end E=T^2(abs(X).^2-sum(abs(x-mxe).^2))/(N^2-N); Ep=T^2(conj(Xpp).Xp-sum(abs(x-mxe).^2))/(N^2-N); Epp=T^2(conj(Xp).Xpp-sum(abs(x-mxe).^2))/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=vxeRE(mod(k+length(RE),length(RE))+1)/sqrt(REp(mod(k+length(REp),length(REp))+1)REpp(mod(k+length(REpp),length(REpp))+1))+abs(mxe)^2; end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
(keine Lomb-Scargle-Methode für komplexe Werte in Matlab und Python)
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x1[j]+=x[i]; X=fft(x1) E=T*2(abs(X)2-sum(abs(x)2))/float(N*2-N) RE=ifft(E)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T-abs(tau[k])) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); x1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x1(j)=x1(j)+x(i); end X=fft(x1); E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); RE=ifft(E)/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/(T-abs(tau(mod(k+K,K)+1))); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x0=zeros(J) x1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=1.0; x1[j]+=x[i]; X=fft(x1) Xp=fft(x0) E=T*2(abs(X)2-sum(abs(x)2))/float(N2-N) Ep=T2(abs(Xp)2-N)/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); x0=zeros([1,J]); x1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+1; x1(j)=x1(j)+x(i); end X=fft(x1); Xp=fft(x0); E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); Ep=T^2(abs(Xp).^2-N)/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/REp(mod(k+length(REp),length(REp))+1); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit lokaler Normierung)		from numpy.fft import * N=len(t) dt=1.0 mxe=mean(x) vxe=var(x) J=int(ceil(2T/float(dt))) x0=zeros(J) x1=zeros(J)+0j x2=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=1.0; x1[j]+=x[i]-mxe; x2[j]+=abs(x[i]-mxe)2; X=fft(x1) Xp=fft(x0) Xpp=fft(x2) E=T2(abs(X)2-sum(abs(x-mxe)2))/float(N2-N) Ep=T2(conj(Xpp)Xp-sum(abs(x-mxe)2))/float(N2-N) Epp=T2(conj(Xp)Xpp-sum(abs(x-mxe)2))/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+abs(mxe)2 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; mxe=mean(x); vxe=var(x,1); J=ceil(2T/dt); x0=zeros([1,J]); x1=zeros([1,J])+0j; x2=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+1; x1(j)=x1(j)+x(i)-mxe; x2(j)=x2(j)+abs(x(i)-mxe)^2; end X=fft(x1); Xp=fft(x0); Xpp=fft(x2); E=T^2(abs(X).^2-sum(abs(x-mxe).^2))/(N^2-N); Ep=T^2(conj(Xpp).Xp-sum(abs(x-mxe).^2))/(N^2-N); Epp=T^2(conj(Xp).Xpp-sum(abs(x-mxe).^2))/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=vxeRE(mod(k+length(RE),length(RE))+1)/sqrt(REp(mod(k+length(REp),length(REp))+1)REpp(mod(k+length(REpp),length(REpp))+1))+abs(mxe)^2; end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) xp=zeros(J)+0j for i in range(0,N): for j in range(int(floor(t[i]/float(dt)))+1,int(floor((t[(i+1)%N]+T((i+1)//N))/float(dt)))+1): xp[j]=x[i] Xp=fft(xp) Ep=dt*2abs(Xp)*2 REp=ifft(Ep)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j c=exp(-dd/float(dt))/(1-exp(-dd/float(dt)))2 for k in range(-(K//2),(K+1)//2): if k==0: R[k]=REp[k]/float(T) else: R[k]=((2c+1)REp[k]-cREp[k-1]-cREp[k+1])/float(T-abs(tau[k])) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); xp=zeros([1,J])+0j; for i=1:N for j=fix(t(i)/dt)+2:fix((t(mod(i,N)+1)+Tfix(i/N))/dt)+1 xp(j)=x(i); end end Xp=fft(xp); Ep=dt^2abs(Xp).^2; REp=ifft(Ep)/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; c=exp(-dd/dt)/(1-exp(-dd/dt))^2; for k=-fix(K/2):fix((K-1)/2) if k==0 R(mod(k+K,K)+1)=REp(mod(k+length(REp),length(REp))+1)/T; else R(mod(k+K,K)+1)=((2c+1)REp(mod(k+length(REp),length(REp))+1)-cREp(mod(k-1+length(REp),length(REp))+1)-cREp(mod(k+1+length(REp),length(REp))+1))/(T-abs(tau(mod(k+K,K)+1))); end end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von Wertepaaren $(t_{i}, x_{i})$ mit den Messzeitpunkten $t_{i}$ und den komplexen Messwerten $x_{i}$ (zwei AR1-Prozess als Beispiel, mit Erwartungswert verschieden von null)		`T=10000.0 Ti=10.0 dr=1.0 mxs=3.0+1.0j vxs=1.0 t=[] x=[] te=0.0 xe=normal(real(mxs),sqrt(vxs/2.0))+1jnormal(imag(mxs),sqrt(vxs/2.0)) while te<T: tp=exponential(1.0/dr) te+=tp if te<T: t.append(te) phi=exp(-tp/Ti) theta=sqrt((1-phi2)vxs/2.0) xe=(xe-mxs)phi+normal(real(mxs),theta)+1jnormal(imag(mxs),theta) x.append(xe) t=array(t) x=array(x) plot(t,real(x),'o',t,imag(x),'o') show()`	`T=10000; Ti=10; dr=1; mxs=3+1j; vxs=1; t=[]; x=[]; te=0; xe=real(mxs)+sqrt(vxs/2)randn()+1j(imag(mxs)+sqrt(vxs/2)randn()); while te<T tp=-log(1-rand())/dr; te=te+tp; if te<T t(end+1)=te; phi=exp(-tp/Ti); theta=sqrt((1-phi^2)vxs/2); xe=(xe-mxs)phi+real(mxs)+thetarandn()+1j(imag(mxs)+thetarandn()); x(end+1)=xe; end end plot(t,real(x),'o',t,imag(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})}^{} (x_{i} - \overline{x}) = \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}^{} x_{i}) - {\overline{x}}^{*} \overline{x} = (\frac{1}{N} \cdot \sum_{i = 0}^{N - 1} {\| x_{i} \|}^{2}) - {\| \overline{x} \|}^{2}$	`var(x)`	`var(x,1)`
Autokorrelationsfunktion und Leistungsdichtespektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ )	$R_{k} = R (τ_{k}) = \frac{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) x_{i}^{*} x_{j}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i})}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K)+0j R0=zeros(K) i=1 while i<N: j=i-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/float(dt))) R1[k]+=conj(x[i])x[j] R0[k]+=1.0 j-=1 i+=1 for k in range(1,(K+1)//2): R1[k]=conj(R1[-k]) R0[k]=R0[-k] R=R1/R0 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K)+0j; R0=zeros(1,K); i=2; while i<=N j=i-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+conj(x(i))x(j); R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+1; j=j-1; end i=i+1; end for k=2:fix((K+1)/2) R1(k)=conj(R1(K-k+2)); R0(k)=R0(K-k+2); end R=R1./R0; plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ , mit lokaler Normierung)	$R_{k} = R (τ_{k}) = \frac{s^{2} [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) {(x_{i} - \overline{x})}^{*} (x_{j} - \overline{x})]}{\sqrt{[\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) {\| x_{i} - \overline{x} \|}^{2}] [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) {\| x_{j} - \overline{x} \|}^{2}]}} + {\| \overline{x} \|}^{2}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t \| < Δ t / 2 \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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 ⌋$	from numpy.fft import * N=len(t) dt=1.0 K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) mxe=mean(x) vxe=var(x) R1=zeros(K)+0j R2=zeros(K) R3=zeros(K) i=1 while i<N: j=i-1 while (j>=0) and (t[j]>t[i]-(K//2+0.5)dt): k=int(round((t[j]-t[i])/float(dt))) R1[k]+=conj(x[i]-mxe)(x[j]-mxe) R2[k]+=abs(x[i]-mxe)2 R3[k]+=abs(x[j]-mxe)2 j-=1 i+=1 for k in range(1,(K+1)//2): R1[k]=conj(R1[-k]) R2[k]=R3[-k] R3[k]=R2[-k] R=vxeR1/sqrt(R2R3)+abs(mxe)*2 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); mxe=mean(x); vxe=var(x,1); R1=zeros(1,K)+0j; R2=zeros(1,K); R3=zeros(1,K); i=2; while i<=N j=i-1; while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)dt) k=round((t(j)-t(i))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+conj(x(i)-mxe)(x(j)-mxe); R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+abs(x(i)-mxe)^2; R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+abs(x(j)-mxe)^2; j=j-1; end i=i+1; end for k=2:fix((K+1)/2) R1(k)=conj(R1(K-k+2)); R2(k)=R3(K-k+2); R3(k)=R2(K-k+2); end R=vxeR1./sqrt(R2.R3)+abs(mxe)^2; plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, ohne Normierung)	$E_{j} = E (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} {\| x_{i} \|}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, k} = R_{E} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{k} = R (τ_{k}) = \frac{R_{E, k}}{T - \| τ_{k} \|}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j for i in range(0,N): X+=x[i]exp(-2jpifpt[i]) E=T*2(abs(X)2-sum(abs(x)2))/float(N*2-N) RE=ifft(E)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T-abs(tau[k])) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; for i=1:N X=X+x(i)exp(-2jpifpt(i)); end E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); RE=ifft(E)/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/(T-abs(tau(mod(k+K,K)+1))); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit Normierung)	$E_{j} = E (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} {\| x_{i} \|}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j}^{'} \|}^{2} - N}$ $X_{j} = \sum_{i = 0}^{N - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, k} = R_{E} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, k}^{'} = R_{E}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{k} = R (τ_{k}) = \frac{R_{E, k}}{R_{E, k}^{'}}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j for i in range(0,N): X+=x[i]exp(-2jpifpt[i]) Xp+=exp(-2jpifpt[i]) E=T2(abs(X)2-sum(abs(x)2))/float(N2-N) Ep=T2(abs(Xp)2-N)/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; for i=1:N X=X+x(i)exp(-2jpifpt(i)); Xp=Xp+exp(-2jpifpt(i)); end E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); Ep=T^2(abs(Xp).^2-N)/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/REp(mod(k+length(REp),length(REp))+1); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit lokaler Normierung)	$E_{j} = E (f_{j}) = \frac{T^{2}}{N (N - 1)} {{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = \frac{T^{2}}{N (N - 1)} {X_{j}^{'' } X_{j}^{'} - \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2}}$ $E_{j}^{''} = E^{''} (f_{j}) = \frac{T^{2}}{N (N - 1)} {X_{j}^{' } X_{j}^{''} - \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} (x_{i} - \overline{x}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{''} = \sum_{i = 0}^{N - 1} {\| x_{i} - \overline{x} \|}^{2} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ 2 T / Δ t ⌉$ $R_{E, k} = R_{E} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, k}^{'} = R_{E}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, k}^{''} = R_{E}^{''} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{j}^{''}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{j}^{''} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{k} = R (τ_{k}) = \frac{s^{2} R_{E, k}}{\sqrt{R_{E, k}^{'} R_{E, k}^{''}}} + {\| \overline{x} \|}^{2}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K < 2 T / Δ t$ $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}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * N=len(t) dt=1.0 mxe=mean(x) vxe=var(x) J=int(ceil(2T/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j Xpp=zeros(size(fp))+0j for i in range(0,N): X+=(x[i]-mxe)exp(-2jpifpt[i]) Xp+=exp(-2jpifpt[i]) Xpp+=abs(x[i]-mxe)2exp(-2jpifpt[i]) E=T2(abs(X)2-sum(abs(x-mxe)2))/float(N2-N) Ep=T2(conj(Xpp)Xp-sum(abs(x-mxe)2))/float(N2-N) Epp=T*2(conj(Xp)Xpp-sum(abs(x-mxe)2))/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+abs(mxe)2 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; mxe=mean(x); vxe=var(x,1); J=ceil(2T/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; Xpp=zeros(size(fp))+0j; for i=1:N X=X+(x(i)-mxe)exp(-2jpifpt(i)); Xp=Xp+exp(-2jpifpt(i)); Xpp=Xpp+abs(x(i)-mxe)^2exp(-2jpifpt(i)); end E=T^2(abs(X).^2-sum(abs(x-mxe).^2))/(N^2-N); Ep=T^2(conj(Xpp).Xp-sum(abs(x-mxe).^2))/(N^2-N); Epp=T^2(conj(Xp).Xpp-sum(abs(x-mxe).^2))/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=vxeRE(mod(k+length(RE),length(RE))+1)/sqrt(REp(mod(k+length(REp),length(REp))+1)REpp(mod(k+length(REpp),length(REpp))+1))+abs(mxe)^2; end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
(keine Lomb-Scargle-Methode für komplexe Werte in Matlab und Python)
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x1[j]+=x[i]; X=fft(x1) E=T*2(abs(X)2-sum(abs(x)2))/float(N*2-N) RE=ifft(E)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T-abs(tau[k])) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); x1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x1(j)=x1(j)+x(i); end X=fft(x1); E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); RE=ifft(E)/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/(T-abs(tau(mod(k+K,K)+1))); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x0=zeros(J) x1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=1.0; x1[j]+=x[i]; X=fft(x1) Xp=fft(x0) E=T*2(abs(X)2-sum(abs(x)2))/float(N2-N) Ep=T2(abs(Xp)2-N)/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); x0=zeros([1,J]); x1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+1; x1(j)=x1(j)+x(i); end X=fft(x1); Xp=fft(x0); E=T^2(abs(X).^2-sum(abs(x).^2))/(N^2-N); Ep=T^2(abs(Xp).^2-N)/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=RE(mod(k+length(RE),length(RE))+1)/REp(mod(k+length(REp),length(REp))+1); end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit lokaler Normierung)		from numpy.fft import * N=len(t) dt=1.0 mxe=mean(x) vxe=var(x) J=int(ceil(2T/float(dt))) x0=zeros(J) x1=zeros(J)+0j x2=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=1.0; x1[j]+=x[i]-mxe; x2[j]+=abs(x[i]-mxe)2; X=fft(x1) Xp=fft(x0) Xpp=fft(x2) E=T2(abs(X)2-sum(abs(x-mxe)2))/float(N2-N) Ep=T2(conj(Xpp)Xp-sum(abs(x-mxe)2))/float(N2-N) Epp=T2(conj(Xp)Xpp-sum(abs(x-mxe)2))/float(N2-N) RE=ifft(E)/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+abs(mxe)2 plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; mxe=mean(x); vxe=var(x,1); J=ceil(2T/dt); x0=zeros([1,J]); x1=zeros([1,J])+0j; x2=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+1; x1(j)=x1(j)+x(i)-mxe; x2(j)=x2(j)+abs(x(i)-mxe)^2; end X=fft(x1); Xp=fft(x0); Xpp=fft(x2); E=T^2(abs(X).^2-sum(abs(x-mxe).^2))/(N^2-N); Ep=T^2(conj(Xpp).Xp-sum(abs(x-mxe).^2))/(N^2-N); Epp=T^2(conj(Xp).Xpp-sum(abs(x-mxe).^2))/(N^2-N); RE=ifft(E)/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; for k=-fix(K/2):fix((K-1)/2) R(mod(k+K,K)+1)=vxeRE(mod(k+length(RE),length(RE))+1)/sqrt(REp(mod(k+length(REp),length(REp))+1)REpp(mod(k+length(REpp),length(REpp))+1))+abs(mxe)^2; end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) xp=zeros(J)+0j for i in range(0,N): for j in range(int(floor(t[i]/float(dt)))+1,int(floor((t[(i+1)%N]+T((i+1)//N))/float(dt)))+1): xp[j]=x[i] Xp=fft(xp) Ep=dt*2abs(Xp)*2 REp=ifft(Ep)/float(dt) K=200 K=minimum(int(ceil(2T/float(dt)))-1,K) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K)+0j c=exp(-dd/float(dt))/(1-exp(-dd/float(dt)))2 for k in range(-(K//2),(K+1)//2): if k==0: R[k]=REp[k]/float(T) else: R[k]=((2c+1)REp[k]-cREp[k-1]-cREp[k+1])/float(T-abs(tau[k])) plot(tau,real(R),'o',arange(-5K,5K)dt/10.0,vxsexp(-abs(arange(-5K,5K)dt/10.0/float(Ti)))+abs(mxs)*2,tau,imag(R),'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) S=dtreal(fft(R)) p1=1-dt/float(Ti) loglog(f,S,'o',arange(1,5K)/float(10Kdt),vxs(1-p12)dt/(1+p1*2-2p1cos(2piarange(1,5K)/float(10Kdt)*dt))) show()	N=length(t); dt=1.0; J=ceil(2T/dt); xp=zeros([1,J])+0j; for i=1:N for j=fix(t(i)/dt)+2:fix((t(mod(i,N)+1)+Tfix(i/N))/dt)+1 xp(j)=x(i); end end Xp=fft(xp); Ep=dt^2abs(Xp).^2; REp=ifft(Ep)/dt; K=200; K=min(ceil(2T/dt)-1,K); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K])+0j; c=exp(-dd/dt)/(1-exp(-dd/dt))^2; for k=-fix(K/2):fix((K-1)/2) if k==0 R(mod(k+K,K)+1)=REp(mod(k+length(REp),length(REp))+1)/T; else R(mod(k+K,K)+1)=((2c+1)REp(mod(k+length(REp),length(REp))+1)-cREp(mod(k-1+length(REp),length(REp))+1)-cREp(mod(k+1+length(REp),length(REp))+1))/(T-abs(tau(mod(k+K,K)+1))); end end plot(tau,real(R),'o',(-5K:5K)dt/10,vxsexp(-abs((-5K:5K)dt/10/Ti))+abs(mxs)^2,tau,imag(R),'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); S=dtreal(fft(R)); p1=1-dt/Ti; loglog(f,S,'o',(1:5K)/(10Kdt),vxs(1-p1^2)dt./(1+p1^2-2p1cos(2pi(1:5K)/(10Kdt)*dt)))