Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für komplexe Energiesignalpaare, mit Gewichtung

		Python	Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)		`from numpy import * from numpy.random import * from matplotlib.pyplot import *`
Generieren von $N_{x}$ komplexen Werten $x_{i}$ mit $i = 0 \dots N_{x} - 1$ und $N_{y}$ komplexen Werten $y_{i}$ mit $i = 0 \dots N_{y} - 1$ (vier verschobene Rechteckimpulse als Beispiel)		`Nx=100 Ny=150 dt=1.0 tx=arange(0,Nx)dt ty=arange(0,Ny)dt x=zeros(Nx)+0j y=zeros(Ny)+0j for i in range(10,70): x[i]+=1 for i in range(30,90): x[i]+=1j for i in range(20,110): y[i]+=1 for i in range(50,130): y[i]+=1j plot(tx,real(x),'o',tx,imag(x),'o',ty,real(y),'o',ty,imag(y),'o') show()`	`Nx=100; Ny=150; dt=1.0; tx=(0:Nx-1)dt; ty=(0:Ny-1)dt; x=zeros(1,Nx)+0j; y=zeros(1,Ny)+0j; for i=10:69 x(i)=x(i)+1; end for i=30:89 x(i)=x(i)+1j; end for i=20:109 y(i)=y(i)+1; end for i=50:129 y(i)=y(i)+1j; end plot(tx,real(x),'o',tx,imag(x),'o',ty,real(y),'o',ty,imag(y),'o')`
Generieren von $N_{x}$ Gewichten $g_{i}$ mit $i = 0 \dots N_{x} - 1$ , passend zu den $x_{i}$ mit $i = 0 \dots N_{x} - 1$ und $N_{y}$ Gewichten $h_{i}$ mit $i = 0 \dots N_{y} - 1$ , passend zu den $y_{i}$ mit $i = 0 \dots N_{y} - 1$ (Kosinusfenster als Beispiel)		`g=sin(pi(arange(0,Nx)+0.5)/float(Nx)) h=sin(pi(arange(0,Ny)+0.5)/float(Ny)) plot(tx,g,'o',ty,h,'o') show()`	`g=sin(pi((0:Nx-1)+0.5)/Nx); h=sin(pi((0:Ny-1)+0.5)/Ny); plot(tx,g,'o',ty,h,'o')`
(Momente verschwinden für Energiesignale)
Kreuzkorrelationsfunktion	$R_{E, y x, k} = R_{E, y x} (τ_{k}) = R_{E, x y}^{} (- τ_{k}) = min (N_{x}, N_{y}, N_{x} + k, N_{y} - k) Δ t \frac{\sum_{i = max (0, - k)}^{min (N_{x} - 1, N_{y} - 1 - k)} g_{i} h_{i + k} x_{i}^{} y_{i + k}}{\sum_{i = max (0, - k)}^{min (N_{x} - 1, N_{y} - 1 - k)} g_{i} h_{i + k}}$ $τ_{k} = k Δ t$ $k = - (N_{x} - 1) \dots N_{y} - 1$ (außerhalb dieses Bereichs null)	direkte Berechnung: `REyx=zeros(Nx+Ny-1)+0j tau=zeros(Nx+Ny-1) for k in range(-(Nx-1),Ny): tau[k]=kdt sum1=0+0j sum0=0 for i in range(maximum(0,-k),minimum(Nx,Ny-k)): sum1+=g[i]h[i+k]conj(x[i])y[i+k] sum0+=g[i]h[i+k] if sum0!=0: REyx[k]=minimum(Nx,minimum(Ny,minimum(Nx+k,Ny-k)))dtsum1/float(sum0) plot(tau,real(REyx),'o',tau,imag(REyx),'o') show()` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `from numpy.fft import tau=roll(arange(-(Nx-1),Ny)dt,Ny) X1=fft(append(gx,zeros(Ny))) Y1=fft(append(hy,zeros(Nx))) X0=fft(append(g,zeros(Ny))) Y0=fft(append(h,zeros(Nx))) Eyx1=dt2conj(X1)Y1 Eyx0=dt2conj(X0)Y0 REyx1=ifft(Eyx1)/float(dt) REyx0=real(ifft(Eyx0))/float(dt) for k in range(-(Nx-1),(Ny-1)): if REyx0[k]!=0: REyx[k]=minimum(Nx,minimum(Ny,minimum(Nx+k,Ny-k)))dt*REyx1[k]/REyx0[k] plot(tau,real(REyx),'o',tau,imag(REyx),'o') show()`	direkte Berechnung: `REyx=zeros(1,Nx+Ny-1)+0j; tau=zeros(1,Nx+Ny-1); for k=-(Nx-1):Ny-1 tau(mod(k,Nx+Ny-1)+1)=kdt; sum1=0+0j; sum0=0; for i=1+max(0,-k):min(Nx-1,Ny-1-k) sum1=sum1+g(i)h(i+k)conj(x(i))y(i+k); sum0=sum0+g(i)h(i+k); end if sum0~=0 REyx(mod(k,Nx+Ny-1)+1)=min(min(Nx,Ny),min(Nx+k,Ny-k))dtsum1/sum0; end end plot(tau,real(REyx),'o',tau,imag(REyx),'o')` Berechnung aus Energiedichtespektrum mittels Wiener-Chintschin-Theorem: `tau=circshift((-(Nx-1):Ny-1)dt,[0;Ny]); X1=fft([g.x,zeros(1,Ny)]); Y1=fft([h.y,zeros(1,Nx)]); X0=fft([g,zeros(1,Ny)]); Y0=fft([h,zeros(1,Nx)]); Eyx1=dt^2conj(X1).Y1; Eyx0=dt^2conj(X0).Y0; REyx1=ifft(Eyx1)/dt; REyx0=real(ifft(Eyx0))/dt; for k=-(Nx-1):Ny-1 if REyx0(mod(k,Nx+Ny-1)+1)~=0 REyx(mod(k,Nx+Ny-1)+1)=min(min(Nx,Ny),min(Nx+k,Ny-k))dtREyx1(mod(k,Nx+Ny)+1)/REyx0(mod(k,Nx+Ny)+1); end end plot(tau,real(REyx),'o',tau,imag(REyx),'o')`
Kreuzenergiedichtespektrum	(wegen der nötigen Entfaltung mit dem Energiedichtespektrum der Gewichte nur über die Korrelationsfunktion) $E_{y x, j} = E_{y x} (f_{j}) = Δ t \cdot F F T {R_{E, y x, k}} = Δ t \cdot \sum_{k = - (N_{x} - 1)}^{N_{y} - 1} R_{E, y x, k} 𝐞^{- 2 π 𝐢 j k / (N_{x} + N_{y})}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{(N_{x} + N_{y}) Δ t}$ $j = - ⌊ (N_{x} + N_{y}) / 2 ⌋ \dots ⌊ (N_{x} + N_{y} - 1) / 2 ⌋$	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(Ny))) Y1=fft(append(hy,zeros(Nx))) Eyx1=dt*2conj(X1)Y1 X0=fft(append(g,zeros(Ny))) Y0=fft(append(h,zeros(Nx))) Eyx0=dt2conj(X0)Y0 REyx1=ifft(Eyx1)/float(dt) REyx0=real(ifft(Eyx0))/float(dt) REyx=zeros(Nx+Ny)+0j for k in range(-Nx,Ny): if REyx0[k]!=0: REyx[k]=minimum(Nx,minimum(Ny,minimum(Nx+k,Ny-k)))dtREyx1[k]/REyx0[k] f=roll(arange(-((Nx+Ny)//2),(Nx+Ny+1)//2)/float((Nx+Ny)dt),(Nx+Ny+1)//2) Eyx=dtfft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show() Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `from numpy.fft import REyx=zeros(Nx+Ny)+0j for k in range(-Nx,Ny): sum1=0+0j sum0=0 for i in range(maximum(0,-k),minimum(Nx,Ny-k)): sum1+=g[i]h[i+k]conj(x[i])y[i+k] sum0+=g[i]h[i+k] if sum0!=0: REyx[k]=minimum(Nx,minimum(Ny,minimum(Nx+k,Ny-k)))dtsum1/float(sum0) f=roll(arange(-((Nx+Ny)//2),(Nx+Ny+1)//2)/float((Nx+Ny)dt),(Nx+Ny+1)//2) Eyx=dtfft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show()`	Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren: X1=fft([g.x,zeros(1,Ny)]); Y1=fft([h.y,zeros(1,Nx)]); Eyx1=dt^2conj(X1).Y1; X0=fft([g,zeros(1,Ny)]); Y0=fft([h,zeros(1,Nx)]); Eyx0=dt^2conj(X0).Y0; REyx1=ifft(Eyx1)/dt; REyx0=real(ifft(Eyx0))/dt; REyx=zeros(Nx+Ny,1)+0j; for k=-Nx:Ny-1 if REyx0(mod(k,Nx+Ny)+1)~=0 REyx(mod(k,Nx+Ny)+1)=min(min(Nx,Ny),min(Nx+k,Ny-k))dtREyx1(mod(k,Nx+Ny)+1)/REyx0(mod(k,Nx+Ny)+1); end end f=circshift(((-fix((Nx+Ny)/2):fix((Nx+Ny-1)/2))/(Nx+Ny)dt),[0;Nx]); Eyx=dtfft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o') Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem: `REyx=zeros(Nx+Ny,1)+0j; for k=-Nx:Ny-1 sum1=0+0j; sum0=0; for i=1+max(0,-k):min(Nx-1,Ny-1-k) sum1=sum1+g(i)h(i+k)conj(x(i))y(i+k); sum0=sum0+g(i)h(i+k); end if sum0~=0 REyx(mod(k,Nx+Ny)+1)=min(min(Nx,Ny),min(Nx+k,Ny-k))dtsum1/sum0; end end f=circshift(((-fix((Nx+Ny)/2):fix((Nx+Ny-1)/2))/(Nx+Ny)dt),[0;Nx]); Eyx=dtfft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o')`