Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für komplexe 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$ komplexen Werten $x_{i}$ mit $i = 0 \dots N - 1$ (zwei verschobene Rechteckimpulse als Beispiel)		`N=100 dt=1.0 t=arange(0,N)*dt x=zeros(N)+0j for i in range(10,70): x[i]+=1 for i in range(30,90): x[i]+=1j plot(t,real(x),'o',t,imag(x),'o') show()`	`N=100; dt=1.0; t=(0:N-1)*dt; x=zeros(1,N)+0j; for i=10:69 x(i)=x(i)+1; end for i=30:89 x(i)=x(i)+1j; end plot(t,real(x),'o',t,imag(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 = max (0, - k)}^{N - 1 - max (0, k)} g_{i} g_{i + k} x_{i}^{*} x_{i + k}}{\sum_{i = max (0, - k)}^{N - 1 - max (0, 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)+0j tau=zeros(2N-1) for k in range(-(N-1),N): tau[k]=kdt sum1=0+0j sum0=0 for i in range(maximum(0,-k),N-maximum(0,k)): sum1+=g[i]g[i+k]conj(x[i])x[i+k] sum0+=g[i]g[i+k] if float(sum0)!=0: RE[k]=(Ndt-abs(tau[k]))sum1/float(sum0) plot(tau,real(RE),'o',tau,imag(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=ifft(E)/float(dt) RE0=real(ifft(Ep))/float(dt) RE=zeros(2N-1)+0j for k in range(-(N-1),N): if RE0[k]!=0: RE[k]=(N-abs(k))dtRE1[k]/RE0[k] plot(tau,real(RE),'o',tau,imag(RE),'o') show()`	direkte Berechnung: `RE=zeros(1,2N-1)+0j; tau=zeros(1,2N-1); for k=-(N-1):N-1 tau(mod(k,2N-1)+1)=kdt; sum1=0+0j; sum0=0; for i=1+max(0,-k):N-max(0,k) sum1=sum1+g(i)g(i+k)conj(x(i))x(i+k); sum0=sum0+g(i)g(i+k); end if sum0~=0 RE(mod(k,2N-1)+1)=(N-abs(k))dtsum1/sum0; end end plot(tau,real(RE),'o',tau,imag(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=ifft(E)/dt; RE0=real(ifft(Ep))/dt; RE=zeros(1,2N-1)+0j; 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,real(RE),'o',tau,imag(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=ifft(E1)/float(dt) RE0=real(ifft(E0))/float(dt) RE=zeros(2N)+0j 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)+0j for k in range(-N,N): sum1=0+0j sum0=0 for i in range(maximum(0,-k),N-maximum(0,k)): sum1+=g[i]g[i+k]conj(x[i])x[i+k] sum0+=g[i]g[i+k] if sum0!=0: RE[k]=(N-abs(k))dtsum1/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=ifft(E1)/dt; RE0=real(ifft(E0))/dt; RE=zeros(2N,1)+0j; 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)+0j; for k=-N:N-1 sum1=0+0j; sum0=0; for i=1+max(0,-k):N-max(0,k) sum1=sum1+g(i)g(i+k)conj(x(i))x(i+k); sum0=sum0+g(i)g(i+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')`