Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für stochastisch und unabhängig abgetastete, reelle Energiesignalpaare, 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_{x, i}, x_{i})$ mit den Messzeitpunkten $t_{x, i}$ und den Messwerten $x_{i}$ mit $i = 0 \dots N_{x} - 1$ sowie den Wertepaaren $(t_{y, i}, y_{i})$ mit den Messzeitpunkten $t_{y, i}$ und den Messwerten $y_{i}$ mit $i = 0 \dots N_{y} - 1$ (zwei Rechteckimpulse)		Tx=100.0 Ty=150.0 drx=1.0 dry=1.5 tx=[] ty=[] x=[] y=[] te=0.0 while te<maximum(Tx,Ty): tp=exponential(1.0/(drx+dry)) te+=tp if te<maximum(Tx,Ty): if random()<drx/float(drx+dry): if te<Tx: tx.append(te) if (te>=10) and (te<70): x.append(1.0) else: x.append(0.0) else: if te<Ty: ty.append(te) if (te>=20) and (te<110): y.append(1.0) else: y.append(0.0) tx=array(tx) ty=array(ty) x=array(x) y=array(y) plot(tx,x,'o',ty,y,'o') show()	`Tx=100; Ty=150; drx=1; dry=1.5; tx=[]; ty=[]; x=[]; y=[]; te=0; while te<max(Tx,Ty) tp=-log(1-rand())/(drx+dry); te=te+tp; if te<max(Tx,Ty) if rand<drx/(drx+dry) if te<Tx tx(end+1)=te; if (te>=10) && (te<70) x(end+1)=1; else x(end+1)=0; end end else if te<Ty ty(end+1)=te; if (te>=20) && (te<110) y(end+1)=1; else y(end+1)=0; end end end end end plot(tx,x,'o',ty,y,'o')`
(Momente verschwinden für Energiesignale)
Kreuzkorrelationsfunktion und Kreuzenergiedichtespektrum über Slotkorrelation	$R_{E, y x, k} = R_{E, y x} (τ_{k}) = R_{E, x y} (- τ_{k}) = min (T_{x}, T_{y}, T_{x} + τ_{k}, T_{y} - τ_{k}) \frac{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 1} b_{k} (t_{j} - t_{i}) x_{i} y_{j}}{\sum_{i = 0}^{N_{x} - 1} \sum_{j = 0}^{N_{y} - 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 = [(T_{x} + T_{y}) / Δ t] - 1$ $E_{y x, j} = E_{y x} (f_{j}) = Δ t \cdot F F T {R_{E, y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{E, y x, 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 * Nx=len(tx) Ny=len(ty) dt=1.0 K=int(round((Tx+Ty)/float(dt)))-1 tau=roll(arange(-int(round(Tx/float(dt)))+1,K-int(round(Tx/float(dt)))+1)dt,K-int(round(Tx/float(dt)))+1) R1=zeros(K) R0=zeros(K) i=0 jb=i while i<Nx: j=jb while (jb<Ny) and (ty[jb]<tx[i]+(K-int(round(Tx/float(dt)))+0.5)dt): jb+=1 j=jb-1 while (j>=0) and (ty[j]>tx[i]-(K-int(round(Tx/float(dt)))-0.5)dt): k=int(round((ty[j]-tx[i])/float(dt))) R1[k]+=x[i]y[j] R0[k]+=1.0 j-=1 i+=1 REyx=zeros(K) for k in range(-(K//2),(K+1)//2): if R0[k]!=0: REyx[k]=R1[k]/float(R0[k])minimum(Tx,minimum(Ty,minimum(Tx+tau[k],Ty-tau[k]))) plot(tau,REyx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Eyx=dt*fft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; K=round((Tx+Ty)/dt)-1; tau=circshift((-round(Tx/dt)+1:K-round(Tx/dt))dt,[0;K-round(Tx/dt)+1]); R1=zeros(1,K); R0=zeros(1,K); i=1; jb=i; while i<=Nx j=jb; while (jb<=Ny) && (ty(jb)<tx(i)+(K-round(Tx/dt)+0.5)dt) jb=jb+1; end j=jb-1; while (j>0) && (ty(j)>tx(i)-(round(Tx/dt)-0.5)dt) k=round((ty(j)-tx(i))/dt); R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+x(i)y(j); R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+1; j=j-1; end i=i+1; end REyx=zeros(1,K); for k=1:K if R0(k)~=0 REyx(k)=R1(k)./R0(k).min(min(Tx,Ty),min(Tx+tau(k),Ty-tau(k))); end end plot(tau,REyx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Eyx=dt*fft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o')
(Slotkorrelation mit lokaler Normierung für Energiesignale nicht möglich)
Kreuzkorrelationsfunktion und Kreuzenergiedichtespektrum über direkte Spektralschätzung (Fourier-Transformation, ohne Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = \frac{T_{x} T_{y} X_{j}^{*} Y_{j}}{N_{x} N_{y}}$ $X_{j} = \sum_{i = 0}^{N_{x} - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j} = \sum_{i = 0}^{N_{y} - 1} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ (T_{x} + T_{y}) / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = R_{E, x y} (- τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [(T_{x} + T_{y}) / Δ t] - 1$ $E_{y x, j} = E_{y x} (f_{j}) = Δ t \cdot F F T {R_{E, y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{E, y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(ceil((Tx+Ty)/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j Y=zeros(size(fp))+0j for i in range(0,Nx): X+=x[i]exp(-2jpifptx[i]) for i in range(0,Ny): Y+=y[i]exp(-2jpifpty[i]) Eyx=TxTyconj(X)Y/float(NxNy) REyx=real(ifft(Eyx))/float(dt) K=int(round((Tx+Ty)/float(dt)))-1 tau=roll(arange(-int(round(Tx/float(dt)))+1,K-int(round(Tx/float(dt)))+1)dt,K-int(round(Tx/float(dt)))+1) REyx=append(REyx[0:K//2+1],REyx[-(K//2):len(REyx)]) plot(tau,REyx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Eyx=dtfft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=ceil((Tx+Ty)/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Y=zeros(size(fp))+0j; for i=1:Nx X=X+x(i)exp(-2jpifptx(i)); end for i=1:Ny Y=Y+y(i)exp(-2jpifpty(i)); end Eyx=TxTyconj(X).Y/(NxNy); REyx=real(ifft(Eyx))/dt; K=round((Tx+Ty)/dt)-1; tau=circshift((-round(Tx/dt)+1:K-round(Tx/dt))dt,[0;K-round(Tx/dt)+1]); REyx=[REyx(1:fix(K/2)+1),REyx(fix(K/2)+3:end)]; plot(tau,REyx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Eyx=dtfft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o')
Kreuzkorrelationsfunktion und Kreuzenergiedichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = \frac{T_{x} T_{y} X_{j}^{} Y_{j}}{N_{x} N_{y}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = \frac{T_{x} T_{y} X_{j}^{' } Y_{j}^{'}}{N_{x} N_{y}}$ $X_{j} = \sum_{i = 0}^{N_{x} - 1} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j} = \sum_{i = 0}^{N_{y} - 1} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ $X_{j}^{'} = \sum_{i = 0}^{N_{x} - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{x, i}}$ $Y_{j}^{'} = \sum_{i = 0}^{N_{y} - 1} 𝐞^{- 2 π 𝐢 f_{j} t_{y, i}}$ (imaginäre Einheit $𝐢$ ) $f_{j} = \frac{j}{J Δ t}$ $j = - ⌊ J / 2 ⌋ \dots ⌊ (J - 1) / 2 ⌋$ $J \geq ⌈ (T_{x} + T_{y}) / Δ t ⌉$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) = R_{E, x y} (- τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k}^{'} = R_{E, y x}^{'} (τ_{k}) = \frac{1}{Δ t} \cdot I F F T {E_{y x, j}^{'}} = \frac{1}{J Δ t} \cdot \sum_{j = - ⌊ J / 2 ⌋}^{⌊ (J - 1) / 2 ⌋} E_{y x, j}^{'} 𝐞^{2 π 𝐢 f_{j} τ_{k} / J}$ $R_{E, y x, k} = R_{E, y x} (τ_{k}) := \frac{min (T_{x}, T_{y}, T_{x} + τ_{k}, T_{y} - τ_{k})}{R_{E, y x, k}^{'}} R_{E, y x, k}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [(T_{x} + T_{y}) / Δ t] - 1$ $E_{y x, j} = E_{y x} (f_{j}) = Δ t \cdot F F T {R_{E, y x, k}} = Δ t \cdot \sum_{k = ⌊ K / 2 ⌋}^{⌊ (K - 1) / 2 ⌋} R_{E, y x, k} 𝐞^{- 2 π 𝐢 j k / K}$ $f_{j} = \frac{j}{K Δ t}$ $j = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$	from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(ceil((Tx+Ty)/float(dt))) fp=roll(arange(-(J//2),(J+1)//2)/float(Jdt),(J+1)//2) X=zeros(size(fp))+0j Y=zeros(size(fp))+0j Xp=zeros(size(fp))+0j Yp=zeros(size(fp))+0j for i in range(0,Nx): X+=x[i]exp(-2jpifptx[i]) Xp+=exp(-2jpifptx[i]) for i in range(0,Ny): Y+=y[i]exp(-2jpifpty[i]) Yp+=exp(-2jpifpty[i]) Eyx=TxTyconj(X)Y/float(NxNy) Eyxp=TxTyconj(Xp)Yp/float(NxNy) REyx=real(ifft(Eyx))/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=int(round((Tx+Ty)/float(dt)))-1 tau=roll(arange(-int(round(Tx/float(dt)))+1,K-int(round(Tx/float(dt)))+1)dt,K-int(round(Tx/float(dt)))+1) REyx=append(REyx[0:K//2+1],REyx[-(K//2):len(REyx)]) REyxp=append(REyxp[0:K//2+1],REyxp[-(K//2):len(REyxp)]) for k in range(-(K//2),(K+1)//2): if REyxp[k]!=0: REyx[k]/=float(REyxp[k]) REyx[k]=minimum(Tx,minimum(Ty,minimum(Tx+tau[k],Ty-tau[k]))) plot(tau,REyx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Eyx=dt*fft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=ceil((Tx+Ty)/dt); fp=circshift((-fix(J/2):fix((J-1)/2))/(Jdt),[0;fix((J+1)/2)]); X=zeros(size(fp))+0j; Y=zeros(size(fp))+0j; Xp=zeros(size(fp))+0j; Yp=zeros(size(fp))+0j; for i=1:Nx X=X+x(i)exp(-2jpifptx(i)); Xp=Xp+exp(-2jpifptx(i)); end for i=1:Ny Y=Y+y(i)exp(-2jpifpty(i)); Yp=Yp+exp(-2jpifpty(i)); end Eyx=TxTyconj(X).Y/(NxNy); Eyxp=TxTyconj(Xp).Yp/(NxNy); REyx=real(ifft(Eyx))/dt; REyxp=real(ifft(Eyxp))/dt; K=round((Tx+Ty)/dt)-1; tau=circshift((-round(Tx/dt)+1:K-round(Tx/dt))dt,[0;K-round(Tx/dt)+1]); REyx=[REyx(1:fix(K/2)+1),REyx(fix(K/2)+3:end)]; REyxp=[REyxp(1:fix(K/2)+1),REyxp(fix(K/2)+3:end)]; for k=1:K if REyxp(k)~=0 REyx(k)=REyx(k)/REyxp(k)min(min(Tx,Ty),min(Tx+tau(k),Ty-tau(k))); end end plot(tau,REyx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Eyx=dt*fft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o')
(Direkte Spektralschätzung mit lokaler Normierung für Energiesignale nicht möglich)
(keine Lomb-Scargle-Methode für Kreuzspektren in Matlab und Python)
Kreuzkorrelationsfunktion und Kreuzenergiedichtespektrum über Zeitquantisierung (ohne Normierung)		from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(ceil((Tx+Ty)/float(dt))) x1=zeros(J) y1=zeros(J) for i in range(0,Nx): j=int(floor((tx[i]-Tfloor(tx[i]/float(T)))/float(dt))) x1[j]+=x[i]; for i in range(0,Ny): j=int(floor((ty[i]-Tfloor(ty[i]/float(T)))/float(dt))) y1[j]+=y[i]; X=fft(x1) Y=fft(y1) Eyx=TxTyconj(X)Y/float(NxNy) REyx=real(ifft(Eyx))/float(dt) K=int(round((Tx+Ty)/float(dt)))-1 tau=roll(arange(-int(round(Tx/float(dt)))+1,K-int(round(Tx/float(dt)))+1)dt,K-int(round(Tx/float(dt)))+1) REyx=append(REyx[0:K//2+1],REyx[-(K//2):len(REyx)]) plot(tau,REyx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Eyx=dt*fft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=ceil((Tx+Ty)/dt); x1=zeros([1,J]); y1=zeros([1,J]); for i=1:Nx j=fix((tx(i)-Tfix(tx(i)/T))/dt)+1; x1(j)=x1(j)+x(i); end for i=1:Ny j=fix((ty(i)-Tfix(ty(i)/T))/dt)+1; y1(j)=y1(j)+y(i); end X=fft(x1); Y=fft(y1); Eyx=TxTyconj(X).Y/(NxNy); REyx=real(ifft(Eyx))/dt; K=round((Tx+Ty)/dt)-1; tau=circshift((-round(Tx/dt)+1:K-round(Tx/dt))dt,[0;K-round(Tx/dt)+1]); REyx=[REyx(1:fix(K/2)+1),REyx(fix(K/2)+3:end)]; plot(tau,REyx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Eyx=dt*fft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o')
Kreuzkorrelationsfunktion und Kreuzenergiedichtespektrum über Zeitquantisierung (mit Normierung)		from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(ceil((Tx+Ty)/float(dt))) x0=zeros(J) y0=zeros(J) x1=zeros(J) y1=zeros(J) for i in range(0,Nx): j=int(floor((tx[i]-Tfloor(tx[i]/float(T)))/float(dt))) x0[j]+=1.0; x1[j]+=x[i]; for i in range(0,Ny): j=int(floor((ty[i]-Tfloor(ty[i]/float(T)))/float(dt))) y0[j]+=1.0; y1[j]+=y[i]; X=fft(x1) Y=fft(y1) Xp=fft(x0) Yp=fft(y0) Eyx=TxTyconj(X)Y/float(NxNy) Eyxp=TxTyconj(Xp)Yp/float(NxNy) REyx=real(ifft(Eyx))/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=int(round((Tx+Ty)/float(dt)))-1 tau=roll(arange(-int(round(Tx/float(dt)))+1,K-int(round(Tx/float(dt)))+1)dt,K-int(round(Tx/float(dt)))+1) REyx=append(REyx[0:K//2+1],REyx[-(K//2):len(REyx)]) REyxp=append(REyxp[0:K//2+1],REyxp[-(K//2):len(REyxp)]) for k in range(-(K//2),(K+1)//2): if REyxp[k]!=0: REyx[k]/=float(REyxp[k]) REyx[k]=minimum(Tx,minimum(Ty,minimum(Tx+tau[k],Ty-tau[k]))) plot(tau,REyx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Eyx=dtfft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=ceil((Tx+Ty)/dt); x0=zeros([1,J]); y0=zeros([1,J]); x1=zeros([1,J]); y1=zeros([1,J]); for i=1:Nx j=fix((tx(i)-Tfix(tx(i)/T))/dt)+1; x0(j)=x0(j)+1; x1(j)=x1(j)+x(i); end for i=1:Ny j=fix((ty(i)-Tfix(ty(i)/T))/dt)+1; y0(j)=y0(j)+1; y1(j)=y1(j)+y(i); end X=fft(x1); Y=fft(y1); Xp=fft(x0); Yp=fft(y0); Eyx=TxTyconj(X).Y/(NxNy); Eyxp=TxTyconj(Xp).Yp/(NxNy); REyx=real(ifft(Eyx))/dt; REyxp=real(ifft(Eyxp))/dt; K=round((Tx+Ty)/dt)-1; tau=circshift((-round(Tx/dt)+1:K-round(Tx/dt))dt,[0;K-round(Tx/dt)+1]); REyx=[REyx(1:fix(K/2)+1),REyx(fix(K/2)+3:end)]; REyxp=[REyxp(1:fix(K/2)+1),REyxp(fix(K/2)+3:end)]; for k=1:K if REyxp(k)~=0 REyx(k)=REyx(k)/REyxp(k)min(min(Tx,Ty),min(Tx+tau(k),Ty-tau(k))); end end plot(tau,REyx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Eyx=dtfft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o')
(Zeitquantisierung mit lokaler Normierung für Energiesignale nicht möglich)
Kreuzkorrelationsfunktion und Kreuzenergiedichtespektrum über Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)		from numpy.fft import * Nx=len(tx) Ny=len(ty) dt=1.0 J=int(ceil((Tx+Ty)/float(dt))) xp=zeros(J) yp=zeros(J) for i in range(0,Nx): for j in range(int(floor(tx[i]/float(dt)))+1,int(floor((tx[(i+1)%Nx]+Tx((i+1)//Nx))/float(dt)))+1): xp[j]=x[i] for i in range(0,Ny): for j in range(int(floor(ty[i]/float(dt)))+1,int(floor((ty[(i+1)%Ny]+Ty((i+1)//Ny))/float(dt)))+1): yp[j]=y[i] Xp=fft(xp) Yp=fft(yp) Eyxp=dt*2conj(Xp)Yp REyxp=real(ifft(Eyxp))/float(dt) K=int(round((Tx+Ty)/float(dt)))-1 tau=roll(arange(-int(round(Tx/float(dt)))+1,K-int(round(Tx/float(dt)))+1)dt,K-int(round(Tx/float(dt)))+1) REyx=zeros(K) cx=exp(-ddx/float(dt))/((1-exp(-ddx/float(dt)))(1-exp(-ddy/float(dt)))) cy=exp(-ddy/float(dt))/((1-exp(-ddx/float(dt)))(1-exp(-ddy/float(dt)))) for k in range(-(K//2),(K+1)//2): REyx[k]=(cx+cy+1)REyxp[k]-cyREyxp[k-1]-cxREyxp[k+1] plot(tau,REyx,'o') show() f=roll(arange(-(K//2),(K+1)//2)/float(Kdt),(K+1)//2) Eyx=dt*fft(REyx) plot(f,real(Eyx),'o',f,imag(Eyx),'o') show()	Nx=length(tx); Ny=length(ty); dt=1.0; J=ceil((Tx+Ty)/dt); xp=zeros([1,J]); yp=zeros([1,J]); for i=1:Nx for j=fix(tx(i)/dt)+2:fix((tx(mod(i,Nx)+1)+Txfix(i/Nx))/dt)+1 xp(j)=x(i); end end for i=1:Ny for j=fix(ty(i)/dt)+2:fix((ty(mod(i,Ny)+1)+Tyfix(i/Ny))/dt)+1 yp(j)=y(i); end end Xp=fft(xp); Yp=fft(yp); Eyxp=dt^2conj(Xp).Yp; REyxp=real(ifft(Eyxp))/dt; K=round((Tx+Ty)/dt)-1; tau=circshift((-round(Tx/dt)+1:K-round(Tx/dt))dt,[0;K-round(Tx/dt)+1]); REyx=zeros([1,K]); cx=exp(-ddx/dt)/((1-exp(-ddx/dt))(1-exp(-ddy/dt))); cy=exp(-ddy/dt)/((1-exp(-ddx/dt))(1-exp(-ddy/dt))); for k=-fix(K/2):fix((K-1)/2) REyx(mod(k+K,K)+1)=(cx+cy+1)REyxp(mod(k+length(REyxp),length(REyxp))+1)-cyREyxp(mod(k-1+length(REyxp),length(REyxp))+1)-cxREyxp(mod(k+1+length(REyxp),length(REyxp))+1); end plot(tau,REyx,'o') f=circshift((-fix(K/2):fix((K-1)/2))/(Kdt),[0;fix((K+1)/2)]); Eyx=dtfft(REyx); plot(f,real(Eyx),'o',f,imag(Eyx),'o')