Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für stochastisch abgetastete, reelle Einzelenergiesignale, 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 Messwerten $x_{i}$ (Rechteckimpuls)		`T=100.0 dr=1.0 t=[] x=[] te=0.0 while te<T: tp=exponential(1.0/dr) te+=tp if te<T: t.append(te) if (te>=10) and (te<90): x.append(1.0) else: x.append(0.0) t=array(t) x=array(x) plot(t,x,'o') show()`	`T=100; dr=1; t=[]; x=[]; te=0; while te<T tp=-log(1-rand())/dr; te=te+tp; if te<T t(end+1)=te; if (te>=10) && (te<90) x(end+1)=1; else x(end+1)=0; end end end plot(t,x,'o')`
(Momente verschwinden für Energiesignale)
Autokorrelationsfunktion und Energiedichtespektrum über Slotkorrelation (ohne Selbstprodukte durch $b_{k} (0) = 0$ )	$R_{E, k} = R_{E} (τ_{k}) = (T - \| τ_{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] - 1$ $E_{j} = E (f_{j}) = Δ t \cdot F F T {R_{E, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{E, 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(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) R1=zeros(K) 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]+=x[i]x[j] R0[k]+=1.0 j-=1 i+=1 for k in range(1,(K+1)//2): R1[k]=R1[-k] R0[k]=R0[-k] RE=zeros(K) for k in range(-(K//2),(K+1)//2): if R0[k]!=0: RE[k]=R1[k]/float(R0[k])(T-abs(tau[k])) plot(tau,RE,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) E=dt*real(fft(RE)) plot(f,E,'o') show()	N=length(t); dt=1.0; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); R1=zeros(1,K); 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)+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)=R1(K-k+2); R0(k)=R0(K-k+2); end RE=zeros(1,K); for k=1:K if R0(k)~=0 RE(k)=R1(k)./R0(k).(T-abs(tau(k))); end end plot(tau,RE,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); E=dt*real(fft(RE)); plot(f,E,'o')
(Slotkorrelation mit lokaler Normierung für Energiesignale nicht möglich)
Autokorrelationsfunktion und Energiedichtespektrum ü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}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [2 T / Δ t] - 1$ $E_{j} = E (f_{j}) = Δ t \cdot F F T {R_{E, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{E, 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=arange(0,J//2+1)/float(Jdt) X=zeros(size(fp))+0j for i in range(0,N): X+=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(x2))/float(N*2-N) for k in range(-(J//2-1),0): E[k]=E[-k] RE=real(ifft(E))/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) RE=append(RE[0:K//2+1],RE[-(K//2):len(RE)]) plot(tau,RE,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) E=dt*real(fft(RE)) plot(f,E,'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=(0:fix(J/2))/(Jdt); X=zeros(size(fp))+0j; for i=1:N X=X+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(x.^2))/(N^2-N); end for k=fix(J/2)+2:J E(k)=E(J-k+2); end RE=real(ifft(E))/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); RE=[RE(1:fix(K/2)+1),RE(fix(K/2)+3:end)]; plot(tau,RE,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); E=dt*real(fft(RE)); plot(f,E,'o')
Autokorrelationsfunktion und Energiedichtespektrum ü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_{E, k} = R_{E} (τ_{k}) := \frac{T - \| τ_{k} \|}{R_{E, k}^{'}} R_{E, k}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [2 T / Δ t] - 1$ $E_{j} = E (f_{j}) = Δ t \cdot F F T {R_{E, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{E, 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=arange(0,J//2+1)/float(Jdt) 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=zeros(J) Ep=zeros(J) for k in range(0,J//2+1): E[k]=T2(abs(X[k])2-sum(x2))/float(N2-N) Ep[k]=T2(abs(Xp[k])2-N)/float(N2-N) 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(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) RE=append(RE[0:K//2+1],RE[-(K//2):len(RE)]) REp=append(REp[0:K//2+1],REp[-(K//2):len(REp)]) for k in range(-(K//2),(K+1)//2): if REp[k]!=0: RE[k]/=float(REp[k]) RE[k]=T-abs(tau[k]) plot(tau,RE,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) E=dtreal(fft(RE)) plot(f,E,'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); fp=(0:fix(J/2))/(Jdt); 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=zeros([1,J]); Ep=zeros([1,J]); for k=1:fix(J/2)+1 E(k)=T^2(abs(X(k)).^2-sum(x.^2))/(N^2-N); Ep(k)=T^2(abs(Xp(k)).^2-N)/(N^2-N); 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(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); RE=[RE(1:fix(K/2)+1),RE(fix(K/2)+3:end)]; REp=[REp(1:fix(K/2)+1),REp(fix(K/2)+3:end)]; for k=1:K if REp(k)~=0 RE(k)=RE(k)/REp(k)(T-abs(tau(k))); end end plot(tau,RE,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); E=dtreal(fft(RE)); plot(f,E,'o')
(Direkte Spektralschätzung mit lokaler Normierung für Energiesignale nicht möglich)
(Lomb-Scargle-Methode in Matlab und Python nur für mittelwertfreie Daten; s. The generalised Lomb-Scargle periodogram.)
Autokorrelationsfunktion und Energiedichtespektrum ü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) 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(x2))/float(N*2-N) RE=real(ifft(E))/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) RE=append(RE[0:K//2+1],RE[-(K//2):len(RE)]) plot(tau,RE,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) E=dt*real(fft(RE)) plot(f,E,'o') show()`	`N=length(t); dt=1.0; J=ceil(2T/dt); x1=zeros([1,J]); 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(x.^2))/(N^2-N); RE=real(ifft(E))/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); RE=[RE(1:fix(K/2)+1),RE(fix(K/2)+3:end)]; plot(tau,RE,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); E=dt*real(fft(RE)); plot(f,E,'o')`
Autokorrelationsfunktion und Energiedichtespektrum ü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) 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(x2))/float(N2-N) Ep=T2(abs(Xp)2-N)/float(N2-N) RE=real(ifft(E))/float(dt) REp=real(ifft(Ep))/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) RE=append(RE[0:K//2+1],RE[-(K//2):len(RE)]) REp=append(REp[0:K//2+1],REp[-(K//2):len(REp)]) for k in range(-(K//2),(K+1)//2): if REp[k]!=0: RE[k]/=float(REp[k]) RE[k]=T-abs(tau[k]) plot(tau,RE,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) E=dtreal(fft(RE)) plot(f,E,'o') show()	N=length(t); dt=1.0; J=ceil(2T/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)+1; x1(j)=x1(j)+x(i); end X=fft(x1); Xp=fft(x0); E=T^2(abs(X).^2-sum(x.^2))/(N^2-N); Ep=T^2(abs(Xp).^2-N)/(N^2-N); RE=real(ifft(E))/dt; REp=real(ifft(Ep))/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); RE=[RE(1:fix(K/2)+1),RE(fix(K/2)+3:end)]; REp=[REp(1:fix(K/2)+1),REp(fix(K/2)+3:end)]; for k=1:K if REp(k)~=0 RE(k)=RE(k)/REp(k)(T-abs(tau(k))); end end plot(tau,RE,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); E=dt*real(fft(RE)); plot(f,E,'o')
(Zeitquantisierung mit lokaler Normierung für Energiesignale nicht möglich)
Autokorrelationsfunktion und Energiedichtespektrum ü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) 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=real(ifft(Ep))/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) RE=zeros(K) c=exp(-dd/float(dt))/(1-exp(-dd/float(dt)))2 for k in range(-(K//2),(K+1)//2): if k==0: RE[k]=REp[k] else: RE[k]=(2c+1)REp[k]-cREp[k-1]-cREp[k+1] plot(tau,RE,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) E=dt*real(fft(RE)) plot(f,E,'o') show()	N=length(t); dt=1.0; J=ceil(2T/dt); xp=zeros([1,J]); 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=real(ifft(Ep))/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); RE=zeros([1,K]); c=exp(-dd/dt)/(1-exp(-dd/dt))^2; RE(1)=REp(1); for k=2:fix(K/2)+1 RE(k)=((2c+1)REp(k)-cREp(k-1)-cREp(k+1)); end for k=fix(K/2)+2:K RE(k)=RE(K-k+2); end plot(tau,RE,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); E=dt*real(fft(RE)); plot(f,E,'o')