Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für stochastisch abgetastete, 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 Wertepaaren $(t_{i}, x_{i})$ mit den Messzeitpunkten $t_{i}$ und den Messwerten $x_{i}$ (Superposition zweier harmonischer Signale als Beispiel, mit einem Erwartungswert verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)		`T=1000.0 dr=5.0 mxs=3.0 axs=1.5 t=[] x=[] te=0.0 while te<T: tp=exponential(1.0/(2dr)) te+=tp xe=axssin(0.1pite)+axssin(0.05pi*te)+mxs if (te<T) and (random()<1/(1+exp(-(xe-mxs)))): t.append(te) x.append(xe) t=array(t) x=array(x) plot(t,x,'o') show()`	`T=1000; dr=5; mxs=3; axs=1.5; t=[]; x=[]; te=0; while te<T tp=-log(1-rand())/(2dr); te=te+tp; xe=axssin(0.1pite)+axssin(0.05pi*te)+mxs; if (te<T) && (rand<1/(1+exp(-(xe-mxs)))) t(end+1)=te; x(end+1)=xe; end end plot(t,x,'o')`
Generieren von Gewichten $g_{i}$ passend zu den $(t_{i}, x_{i})$ (Inverse der verschobenen Sigmoidfunktion)		`g=1+exp(-(x-mxs)) plot(t,g,'o') show()`	`g=1+exp(-(x-mxs)); 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`
Autokorrelationsfunktion und Leistungsspektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ , auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$R_{k} = R (τ_{k}) = \frac{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} x_{i} x_{j}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j}}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t + c T \| < Δ t / 2 für irgendein ganzes c \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K) R0=zeros(K) for i in range(1,N): for j in range(0,i): k=int(round((t[j]-t[i]-Tround((t[j]-t[i])/float(T)))/float(dt))) R1[k]+=g[i]g[j]x[i]x[j] R0[k]+=g[i]g[j] for k in range(-(K//2),1): R1[k]+=R1[-k] R0[k]+=R0[-k] for k in range(1,(K+1)//2): R1[k]=R1[-k] R0[k]=R0[-k] R=R1/R0 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K); R0=zeros(1,K); for i=2:N for j=1:i-1 k=round((t(j)-t(i)-Tround((t(j)-t(i))/T))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)g(j)x(i)x(j); R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+g(i)g(j); end end R1(1)=2R1(1); R0(1)=2R0(1); for k=fix((K+1)/2)+1:K R1(k)=R1(k)+R1(K-k+2); R0(k)=R0(k)+R0(K-k+2); end for k=2:fix((K+1)/2) R1(k)=R1(K-k+2); R0(k)=R0(K-k+2); end R=R1./R0; plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ , mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$R_{k} = R (τ_{k}) = \frac{s^{2} [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} (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}) g_{i} g_{j} {(x_{i} - \overline{x})}^{2}] [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {(x_{j} - \overline{x})}^{2}]}} + {\overline{x}}^{2}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t + c T \| < Δ t / 2 für irgendein ganzes c \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) mxe=sum(gx)/float(sum(g)) vxe=sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 R1=zeros(K) R2=zeros(K) R3=zeros(K) for i in range(1,N): for j in range(0,i): k=int(round((t[j]-t[i]-Tround((t[j]-t[i])/float(T)))/float(dt))) R1[k]+=g[i]g[j](x[i]-mxe)(x[j]-mxe) R2[k]+=g[i]g[j](x[i]-mxe)2 R3[k]+=g[i]g[j](x[j]-mxe)2 for k in range(-(K//2),1): R1[k]+=R1[-k] R2[k]+=R3[-k] if k%K!=(-k)%K: R3[k]+=R2[-k] else: R3[k]=R2[-k] for k in range(1,(K+1)//2): R1[k]=R1[-k] R2[k]=R3[-k] R3[k]=R2[-k] R=vxeR1/sqrt(R2R3)+mxe2 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; R1=zeros(1,K); R2=zeros(1,K); R3=zeros(1,K); for i=2:N for j=1:i-1 k=round((t(j)-t(i)-Tround((t(j)-t(i))/T))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)g(j)(x(i)-mxe)(x(j)-mxe); R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+g(i)g(j)(x(i)-mxe)^2; R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+g(i)g(j)(x(j)-mxe)^2; end end R1(1)=2R1(1); R2(1)=R2(1)+R3(1); R3(1)=R2(1); for k=fix((K+1)/2)+1:K R1(k)=R1(k)+R1(K-k+2); R2(k)=R2(k)+R3(K-k+2); if mod(k,K)~=mod(K-k+2,K) R3(k)=R3(k)+R2(K-k+2); else R3(k)=R2(K-k+2); end end for k=2:fix((K+1)/2) R1(k)=R1(K-k+2); R2(k)=R3(K-k+2); R3(k)=R2(K-k+2); end R=vxeR1./sqrt(R2.R3)+mxe^2; plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, ohne Normierung)	$E_{j} = E (f_{j}) = T^{2} \frac{{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{2}}{{(\sum_{i = 0}^{N - 1} g_{i})}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J = [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{1}{T} R_{E, k}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = J = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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(round(T/float(dt))) fp=arange(0,J//2+1)/float(Jdt) X=zeros(size(fp))+0j for i in range(0,N): X+=g[i]x[i]exp(-2jpifpt[i]) E=zeros(J) for k in range(0,J//2+1): E[k]=T*2(abs(X[k])*2-sum((gx)2))/float(sum(g)2-sum(g*2)) for k in range(-(J//2-1),0): E[k]=E[-k] RE=real(ifft(E))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; J=round(T/dt); fp=(0:fix(J/2))/(Jdt); X=zeros(size(fp))+0j; for i=1:N X=X+g(i)x(i)exp(-2jpifpt(i)); end E=zeros([1,J]); for k=1:fix(J/2)+1 E(k)=T^2(abs(X(k)).^2-sum((g.x).^2))/(sum(g)^2-sum(g.^2)); end for k=fix(J/2)+2:J E(k)=E(J-k+2); end RE=real(ifft(E))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/T; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$E_{j} = E (f_{j}) = T^{2} \frac{{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = T^{2} \frac{{\| X_{j}^{'} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J = [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 = J = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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(round(T/float(dt))) fp=arange(0,J//2+1)/float(Jdt) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j for i in range(0,N): X+=g[i]x[i]exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) E=zeros(J) Ep=zeros(J) for k in range(0,J//2+1): E[k]=T*2(abs(X[k])*2-sum((gx)2))/(abs(Xp[0])2-sum(g2)) Ep[k]=T2(abs(Xp[k])2-sum(g2))/(abs(Xp[0])2-sum(g2)) for k in range(-(J//2-1),0): E[k]=E[-k] Ep[k]=Ep[-k] RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; J=round(T/dt); fp=(0:fix(J/2))/(Jdt); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; for i=1:N X=X+g(i)x(i)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(i)); end E=zeros([1,J]); Ep=zeros([1,J]); for k=1:fix(J/2)+1 E(k)=T^2(abs(X(k)).^2-sum((g.x).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep(k)=T^2(abs(Xp(k)).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2)); end for k=fix(J/2)+2:J E(k)=E(J-k+2); Ep(k)=Ep(J-k+2); end RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/REp(k); end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$E_{j} = E (f_{j}) = T^{2} \frac{{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = T^{2} \frac{X_{j}^{'' } X_{j}^{'} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{j}^{''} = E^{''} (f_{j}) = T^{2} \frac{X_{j}^{' } X_{j}^{''} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} (x_{i} - \overline{x}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{''} = \sum_{i = 0}^{N - 1} g_{i} {(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 = [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 = J = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 J=int(round(T/float(dt))) fp=arange(0,J//2+1)/float(Jdt) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j Xpp=zeros(size(fp))+0j for i in range(0,N): X+=g[i](x[i]-mxe)exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) Xpp+=g[i](x[i]-mxe)*2exp(-2jpifpt[i]) E=zeros(J) Ep=zeros(J) Epp=zeros(J) for k in range(0,J//2+1): E[k]=T2(abs(X[k])*2-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Ep[k]=T2(conj(Xpp[k])Xp[k]-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Epp[k]=T2(conj(Xp[k])Xpp[k]-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g*2)) for k in range(-(J//2-1),0): E[k]=E[-k] Ep[k]=Ep[-k] Epp[k]=Epp[-k] RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+mxe*2 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; J=round(T/dt); fp=(0:fix(J/2))/(Jdt); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; Xpp=zeros(size(fp))+0j; for i=1:N X=X+g(i)(x(i)-mxe)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(i)); Xpp=Xpp+g(i)(x(i)-mxe)^2exp(-2jpifpt(i)); end E=zeros([1,J]); Ep=zeros([1,J]); Epp=zeros([1,J]); for k=1:fix(J/2)+1 E(k)=T^2(abs(X(k)).^2-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep(k)=T^2(conj(Xpp(k))Xp(k)-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Epp(k)=T^2(conj(Xp(k))Xpp(k)-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); end for k=fix(J/2)+2:J E(k)=E(J-k+2); Ep(k)=Ep(J-k+2); Epp(k)=Epp(J-k+2); end RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=vxeRE(k)/sqrt(REp(k)REpp(k))+mxe^2; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
(keine Gewichtung für Lomb-Scargle-Methode in Matlab und Python; s. The generalised Lomb-Scargle periodogram.)
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(round(T/float(dt))) x1=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x1[j]+=g[i]x[i]; X=fft(x1) E=T*2(abs(X)*2-sum((gx)2))/float(sum(g)2-sum(g*2)) RE=real(ifft(E))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	`N=length(t); dt=1.0; J=round(T/dt); x1=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x1(j)=x1(j)+g(i)x(i); end X=fft(x1); E=T^2(abs(X).^2-sum((g.x).^2))/(sum(g)^2-sum(g.^2)); RE=real(ifft(E))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/T; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')`
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)		from numpy.fft import * N=len(t) dt=1.0 J=int(round(T/float(dt))) x0=zeros(J) x1=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=g[i]; x1[j]+=g[i]x[i]; X=fft(x1) Xp=fft(x0) E=T*2(abs(X)*2-sum((gx)2))/(abs(Xp[0])2-sum(g2)) Ep=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; J=round(T/dt); x0=zeros([1,J]); x1=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+g(i); x1(j)=x1(j)+g(i)x(i); end X=fft(x1); Xp=fft(x0); E=T^2(abs(X).^2-sum((g.x).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep=T^2(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2)); RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/REp(k); end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)		from numpy.fft import * N=len(t) dt=1.0 mxe=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 J=int(round(T/float(dt))) x0=zeros(J) x1=zeros(J) x2=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=g[i]; x1[j]+=g[i](x[i]-mxe); x2[j]+=g[i](x[i]-mxe)2; X=fft(x1) Xp=fft(x0) Xpp=fft(x2) E=T2(abs(X)2-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Ep=T2(conj(Xpp)Xp-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Epp=T2(conj(Xp)Xpp-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g*2)) RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+mxe*2 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; J=round(T/dt); x0=zeros([1,J]); x1=zeros([1,J]); x2=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+g(i); x1(j)=x1(j)+g(i)(x(i)-mxe); x2(j)=x2(j)+g(i)(x(i)-mxe)^2; end X=fft(x1); Xp=fft(x0); Xpp=fft(x2); E=T^2(abs(X).^2-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep=T^2(conj(Xpp).Xp-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Epp=T^2(conj(Xp).Xpp-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=vxeRE(k)/sqrt(REp(k)REpp(k))+mxe^2; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
(keine Gewichtung für Interpolation)

		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 Messwerten $x_{i}$ (Superposition zweier harmonischer Signale als Beispiel, mit einem Erwartungswert verschieden von null) mit einer vom Wert abhängigen Annahmewahrscheinlichkeit (verschobene Sigmoidfunktion als Beispiel)		`T=1000.0 dr=5.0 mxs=3.0 axs=1.5 t=[] x=[] te=0.0 while te<T: tp=exponential(1.0/(2dr)) te+=tp xe=axssin(0.1pite)+axssin(0.05pi*te)+mxs if (te<T) and (random()<1/(1+exp(-(xe-mxs)))): t.append(te) x.append(xe) t=array(t) x=array(x) plot(t,x,'o') show()`	`T=1000; dr=5; mxs=3; axs=1.5; t=[]; x=[]; te=0; while te<T tp=-log(1-rand())/(2dr); te=te+tp; xe=axssin(0.1pite)+axssin(0.05pi*te)+mxs; if (te<T) && (rand<1/(1+exp(-(xe-mxs)))) t(end+1)=te; x(end+1)=xe; end end plot(t,x,'o')`
Generieren von Gewichten $g_{i}$ passend zu den $(t_{i}, x_{i})$ (Inverse der verschobenen Sigmoidfunktion)		`g=1+exp(-(x-mxs)) plot(t,g,'o') show()`	`g=1+exp(-(x-mxs)); 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`
Autokorrelationsfunktion und Leistungsspektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ , auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$R_{k} = R (τ_{k}) = \frac{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} x_{i} x_{j}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j}}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t + c T \| < Δ t / 2 für irgendein ganzes c \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K) R0=zeros(K) for i in range(1,N): for j in range(0,i): k=int(round((t[j]-t[i]-Tround((t[j]-t[i])/float(T)))/float(dt))) R1[k]+=g[i]g[j]x[i]x[j] R0[k]+=g[i]g[j] for k in range(-(K//2),1): R1[k]+=R1[-k] R0[k]+=R0[-k] for k in range(1,(K+1)//2): R1[k]=R1[-k] R0[k]=R0[-k] R=R1/R0 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K); R0=zeros(1,K); for i=2:N for j=1:i-1 k=round((t(j)-t(i)-Tround((t(j)-t(i))/T))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)g(j)x(i)x(j); R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+g(i)g(j); end end R1(1)=2R1(1); R0(1)=2R0(1); for k=fix((K+1)/2)+1:K R1(k)=R1(k)+R1(K-k+2); R0(k)=R0(k)+R0(K-k+2); end for k=2:fix((K+1)/2) R1(k)=R1(K-k+2); R0(k)=R0(K-k+2); end R=R1./R0; plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ , mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$R_{k} = R (τ_{k}) = \frac{s^{2} [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} (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}) g_{i} g_{j} {(x_{i} - \overline{x})}^{2}] [\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} {(x_{j} - \overline{x})}^{2}]}} + {\overline{x}}^{2}$ $b_{k} (t) = {\begin{matrix} 1 & falls 0 < \| t + c T \| < Δ t / 2 für irgendein ganzes c \\ 0 & sonst \end{matrix}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) mxe=sum(gx)/float(sum(g)) vxe=sum(gx2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 R1=zeros(K) R2=zeros(K) R3=zeros(K) for i in range(1,N): for j in range(0,i): k=int(round((t[j]-t[i]-Tround((t[j]-t[i])/float(T)))/float(dt))) R1[k]+=g[i]g[j](x[i]-mxe)(x[j]-mxe) R2[k]+=g[i]g[j](x[i]-mxe)2 R3[k]+=g[i]g[j](x[j]-mxe)2 for k in range(-(K//2),1): R1[k]+=R1[-k] R2[k]+=R3[-k] if k%K!=(-k)%K: R3[k]+=R2[-k] else: R3[k]=R2[-k] for k in range(1,(K+1)//2): R1[k]=R1[-k] R2[k]=R3[-k] R3[k]=R2[-k] R=vxeR1/sqrt(R2R3)+mxe2 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; R1=zeros(1,K); R2=zeros(1,K); R3=zeros(1,K); for i=2:N for j=1:i-1 k=round((t(j)-t(i)-Tround((t(j)-t(i))/T))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+g(i)g(j)(x(i)-mxe)(x(j)-mxe); R2(mod(K+k,K)+1)=R2(mod(K+k,K)+1)+g(i)g(j)(x(i)-mxe)^2; R3(mod(K+k,K)+1)=R3(mod(K+k,K)+1)+g(i)g(j)(x(j)-mxe)^2; end end R1(1)=2R1(1); R2(1)=R2(1)+R3(1); R3(1)=R2(1); for k=fix((K+1)/2)+1:K R1(k)=R1(k)+R1(K-k+2); R2(k)=R2(k)+R3(K-k+2); if mod(k,K)~=mod(K-k+2,K) R3(k)=R3(k)+R2(K-k+2); else R3(k)=R2(K-k+2); end end for k=2:fix((K+1)/2) R1(k)=R1(K-k+2); R2(k)=R3(K-k+2); R3(k)=R2(K-k+2); end R=vxeR1./sqrt(R2.R3)+mxe^2; plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, ohne Normierung)	$E_{j} = E (f_{j}) = T^{2} \frac{{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{2}}{{(\sum_{i = 0}^{N - 1} g_{i})}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J = [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{1}{T} R_{E, k}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = J = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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(round(T/float(dt))) fp=arange(0,J//2+1)/float(Jdt) X=zeros(size(fp))+0j for i in range(0,N): X+=g[i]x[i]exp(-2jpifpt[i]) E=zeros(J) for k in range(0,J//2+1): E[k]=T*2(abs(X[k])*2-sum((gx)2))/float(sum(g)2-sum(g*2)) for k in range(-(J//2-1),0): E[k]=E[-k] RE=real(ifft(E))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; J=round(T/dt); fp=(0:fix(J/2))/(Jdt); X=zeros(size(fp))+0j; for i=1:N X=X+g(i)x(i)exp(-2jpifpt(i)); end E=zeros([1,J]); for k=1:fix(J/2)+1 E(k)=T^2(abs(X(k)).^2-sum((g.x).^2))/(sum(g)^2-sum(g.^2)); end for k=fix(J/2)+2:J E(k)=E(J-k+2); end RE=real(ifft(E))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/T; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$E_{j} = E (f_{j}) = T^{2} \frac{{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = T^{2} \frac{{\| X_{j}^{'} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J = [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 = J = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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(round(T/float(dt))) fp=arange(0,J//2+1)/float(Jdt) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j for i in range(0,N): X+=g[i]x[i]exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) E=zeros(J) Ep=zeros(J) for k in range(0,J//2+1): E[k]=T*2(abs(X[k])*2-sum((gx)2))/(abs(Xp[0])2-sum(g2)) Ep[k]=T2(abs(Xp[k])2-sum(g2))/(abs(Xp[0])2-sum(g2)) for k in range(-(J//2-1),0): E[k]=E[-k] Ep[k]=Ep[-k] RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; J=round(T/dt); fp=(0:fix(J/2))/(Jdt); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; for i=1:N X=X+g(i)x(i)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(i)); end E=zeros([1,J]); Ep=zeros([1,J]); for k=1:fix(J/2)+1 E(k)=T^2(abs(X(k)).^2-sum((g.x).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep(k)=T^2(abs(Xp(k)).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2)); end for k=fix(J/2)+2:J E(k)=E(J-k+2); Ep(k)=Ep(J-k+2); end RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/REp(k); end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)	$E_{j} = E (f_{j}) = T^{2} \frac{{\| X_{j} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{j}^{'} = E^{'} (f_{j}) = T^{2} \frac{X_{j}^{'' } X_{j}^{'} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{j}^{''} = E^{''} (f_{j}) = T^{2} \frac{X_{j}^{' } X_{j}^{''} - \sum_{i = 0}^{N - 1} g_{i}^{2} {(x_{i} - \overline{x})}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} (x_{i} - \overline{x}) 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{''} = \sum_{i = 0}^{N - 1} g_{i} {(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 = [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 = J = [T / Δ t]$ $P_{j} = P (f_{j}) = \frac{1}{K} \cdot F F T {R_{k}} = \frac{1}{K} \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=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 J=int(round(T/float(dt))) fp=arange(0,J//2+1)/float(Jdt) X=zeros(size(fp))+0j Xp=zeros(size(fp))+0j Xpp=zeros(size(fp))+0j for i in range(0,N): X+=g[i](x[i]-mxe)exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) Xpp+=g[i](x[i]-mxe)*2exp(-2jpifpt[i]) E=zeros(J) Ep=zeros(J) Epp=zeros(J) for k in range(0,J//2+1): E[k]=T2(abs(X[k])*2-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Ep[k]=T2(conj(Xpp[k])Xp[k]-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Epp[k]=T2(conj(Xp[k])Xpp[k]-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g*2)) for k in range(-(J//2-1),0): E[k]=E[-k] Ep[k]=Ep[-k] Epp[k]=Epp[-k] RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+mxe*2 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; J=round(T/dt); fp=(0:fix(J/2))/(Jdt); X=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; Xpp=zeros(size(fp))+0j; for i=1:N X=X+g(i)(x(i)-mxe)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(i)); Xpp=Xpp+g(i)(x(i)-mxe)^2exp(-2jpifpt(i)); end E=zeros([1,J]); Ep=zeros([1,J]); Epp=zeros([1,J]); for k=1:fix(J/2)+1 E(k)=T^2(abs(X(k)).^2-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep(k)=T^2(conj(Xpp(k))Xp(k)-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Epp(k)=T^2(conj(Xp(k))Xpp(k)-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); end for k=fix(J/2)+2:J E(k)=E(J-k+2); Ep(k)=Ep(J-k+2); Epp(k)=Epp(J-k+2); end RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=vxeRE(k)/sqrt(REp(k)REpp(k))+mxe^2; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
(keine Gewichtung für Lomb-Scargle-Methode in Matlab und Python; s. The generalised Lomb-Scargle periodogram.)
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(round(T/float(dt))) x1=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x1[j]+=g[i]x[i]; X=fft(x1) E=T*2(abs(X)*2-sum((gx)2))/float(sum(g)2-sum(g*2)) RE=real(ifft(E))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(T) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	`N=length(t); dt=1.0; J=round(T/dt); x1=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x1(j)=x1(j)+g(i)x(i); end X=fft(x1); E=T^2(abs(X).^2-sum((g.x).^2))/(sum(g)^2-sum(g.^2)); RE=real(ifft(E))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/T; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')`
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)		from numpy.fft import * N=len(t) dt=1.0 J=int(round(T/float(dt))) x0=zeros(J) x1=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=g[i]; x1[j]+=g[i]x[i]; X=fft(x1) Xp=fft(x0) E=T*2(abs(X)*2-sum((gx)2))/(abs(Xp[0])2-sum(g2)) Ep=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=RE[k]/float(REp[k]) plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; J=round(T/dt); x0=zeros([1,J]); x1=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+g(i); x1(j)=x1(j)+g(i)x(i); end X=fft(x1); Xp=fft(x0); E=T^2(abs(X).^2-sum((g.x).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep=T^2(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2)); RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=RE(k)/REp(k); end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
Autokorrelationsfunktion und Leistungsspektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit lokaler Normierung, auch für den Fall, dass die anzuwendende Periode kleiner als das Beobachtungsintervall ist)		from numpy.fft import * N=len(t) dt=1.0 mxe=sum(gx)/float(sum(g)) vxe=sum(gx*2)/float(sum(g))-(sum(gx)/float(sum(g)))*2 J=int(round(T/float(dt))) x0=zeros(J) x1=zeros(J) x2=zeros(J) for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=g[i]; x1[j]+=g[i](x[i]-mxe); x2[j]+=g[i](x[i]-mxe)2; X=fft(x1) Xp=fft(x0) Xpp=fft(x2) E=T2(abs(X)2-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Ep=T2(conj(Xpp)Xp-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g2)) Epp=T2(conj(Xp)Xpp-sum((g(x-mxe))2))/(abs(Xp[0])2-sum(g*2)) RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) REpp=real(ifft(Epp))/float(dt) K=int(round(T/float(dt))) tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R=zeros(K) for k in range(-(K//2),(K+1)//2): R[k]=vxeRE[k]/sqrt(REp[k]REpp[k])+mxe*2 plot(tau,R,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) P=real(fft(R))/float(K) plot(f,P,'o') show()	N=length(t); dt=1.0; mxe=sum(g.x)/sum(g); vxe=sum(g.x.^2)/sum(g)-(sum(g.x)/sum(g))^2; J=round(T/dt); x0=zeros([1,J]); x1=zeros([1,J]); x2=zeros([1,J]); for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+g(i); x1(j)=x1(j)+g(i)(x(i)-mxe); x2(j)=x2(j)+g(i)(x(i)-mxe)^2; end X=fft(x1); Xp=fft(x0); Xpp=fft(x2); E=T^2(abs(X).^2-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Ep=T^2(conj(Xpp).Xp-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); Epp=T^2(conj(Xp).Xpp-sum((g.(x-mxe)).^2))/(abs(Xp(1))^2-sum(g.^2)); RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; REpp=real(ifft(Epp))/dt; K=round(T/dt); tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R=zeros([1,K]); for k=1:fix(K/2)+1 R(k)=vxeRE(k)/sqrt(REp(k)REpp(k))+mxe^2; end for k=fix(K/2)+2:K R(k)=R(K-k+2); end plot(tau,R,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); P=real(fft(R))/K; plot(f,P,'o')
(keine Gewichtung für Interpolation)