Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für reelle Einzelenergiesignale, mit Gewichtung

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von $N$ Werten $x_{i}$ mit $i = 0 \dots N - 1$ (Rechteckimpuls als Beispiel)		`N=100 dt=1.0 t=arange(0,N)*dt x=zeros(N) for i in range(10,90): x[i]=1 plot(t,x,'o') show()`	`N=100; dt=1.0; t=(0:N-1)*dt; x=zeros(1,N); for i=10:89 x(i)=1; end plot(t,x,'o')`
Generieren von $N$ Gewichten $g_{i}$ mit $i = 0 \dots N - 1$ , passend zu den $x_{i}$ mit $i = 0 \dots N - 1$ (Kosinusfenster als Beispiel)		`g=sin(pi*(arange(0,N)+0.5)/float(N)) plot(t,g,'o') show()`	`g=sin(pi*((0:N-1)+0.5)/N); plot(t,g,'o')`
(Momente verschwinden für Energiesignale)
(Auto-)Korrelationsfunktion	$R_{E, k} = R_{E} (τ_{k}) = (N Δ t - \| τ_{k} \|) \frac{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|} x_{i} x_{i + \| k \|}}{\sum_{i = 0}^{N - 1 - \| k \|} g_{i} g_{i + \| k \|}}$ $τ_{k} = k Δ t$ $k = - (N - 1) \dots N - 1$ (außerhalb dieses Bereichs null)	direkte Berechnung: `RE=zeros(2N-1) tau=zeros(2N-1) for k in range(-(N-1),N): tau[k]=kdt sum1=0 sum0=0 for i in range(0,N-abs(k)): sum1+=g[i]g[i+abs(k)]x[i]x[i+abs(k)] sum0+=g[i]g[i+abs(k)] if sum0!=0: RE[k]=(Ndt-abs(tau[k]))sum1/float(sum0) plot(tau,RE,'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=roll(arange(-(N-1),N)dt,N) X=fft(append(gx,zeros(N))) Xp=fft(append(g,zeros(N))) E=dt*2abs(X)2 Ep=dt2abs(Xp)2 RE1=real(ifft(E))/float(dt) RE0=real(ifft(Ep))/float(dt) RE=zeros(2N-1) for k in range(-(N-1),N): if RE0[k]!=0: RE[k]=(N-abs(k))dtRE1[k]/RE0[k] plot(tau,RE,'o') show()`	direkte Berechnung: `RE=zeros(1,2N-1); tau=zeros(1,2N-1); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0; sum0=0; for i=1:N-abs(k) sum1=sum1+g(i)g(i+abs(k))x(i)x(i+abs(k)); sum0=sum0+g(i)g(i+abs(k)); end if sum0~=0 RE(mod(k,2N-1)+1)=(N-abs(k))dtsum1/sum0; end end plot(tau,RE,'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `tau=circshift((-(N-1):N-1)dt,[0;N]); X=fft([g.x,zeros(1,N)]); Xp=fft([g,zeros(1,N)]); E=dt^2abs(X).^2; Ep=dt^2abs(Xp).^2; RE1=real(ifft(E))/dt; RE0=real(ifft(Ep))/dt; RE=zeros(1,2N-1); for k=-(N-1):N-1 if RE0(mod(k,2N)+1)~=0 RE(mod(k,2N-1)+1)=(N-abs(k))dtRE1(mod(k,2N)+1)/RE0(mod(k,2N)+1); end end plot(tau,RE,'o')`
(Auto-)Energiedichtespektrum	(wegen der nötigen Entfaltung mit dem Energiedichtespektrum der Gewichte nur über die Korrelationsfunktion) $E_{j} = E (f_{j}) = Δ t \cdot F F T {R_{E, k}} = Δ t \cdot \sum_{k = - (N - 1)}^{N - 1} R_{E, k} 𝐞^{- π 𝐢 j k / N}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{2 N Δ t}$ $j = - N \dots N - 1$	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `from numpy.fft import * X1=fft(append(gx,zeros(N))) E1=dt2abs(X1)2 X0=fft(append(g,zeros(N))) E0=dt2abs(X0)2 RE1=real(ifft(E1))/float(dt) RE0=real(ifft(E0))/float(dt) RE=zeros(2N) for k in range(-N,N): if RE0[k]!=0: RE[k]=(N-abs(k))dtRE1[k]/RE0[k] f=roll(arange(-N,N)/float(2Ndt),N) E=dtreal(fft(RE)) plot(f,E,'o') show()` Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import RE=zeros(2N) for k in range(-N,N): sum1=0 sum0=0 for i in range(0,N-abs(k)): sum1+=g[i]g[i+abs(k)]x[i]x[i+abs(k)] sum0+=g[i]g[i+abs(k)] if sum0!=0: RE[k]=(N-abs(k))dtsum1/float(sum0) f=roll(arange(-N,N)/float(2Ndt),N) E=dtreal(fft(RE)) plot(f,E,'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: `X1=fft([g.x,zeros(1,N)]); E1=dt^2abs(X1).^2; X0=fft([g,zeros(1,N)]); E0=dt^2abs(X0).^2; RE1=real(ifft(E1))/dt; RE0=real(ifft(E0))/dt; RE=zeros(2N,1); for k=-N:N-1 if RE0(mod(k,2N)+1)~=0 RE(mod(k,2N)+1)=(N-abs(k))dtRE1(mod(k,2N)+1)/RE0(mod(k,2N)+1); end end f=circshift((-N:N-1)/(2Ndt),[0;N]); E=dtreal(fft(RE)); plot(f,E,'o')` Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `RE=zeros(2N,1); for k=-N:N-1 sum1=0; sum0=0; for i=1:N-abs(k) sum1=sum1+g(i)g(i+abs(k))x(i)x(i+abs(k)); sum0=sum0+g(i)g(i+abs(k)); end if sum0~=0 RE(mod(k,2N)+1)=(N-abs(k))dtsum1/sum0; end end f=circshift((-N:N-1)/(2Ndt),[0;N]); E=dtreal(fft(RE)); plot(f,E,'o')`