Holger Nobach - nambis.de

Signal- und Messdatenverarbeitung

Codes für stochastisch und gemeinsam abgetastete, 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 Wertetripeln $(t_{i}, x_{i}, y_{i})$ mit den Messzeitpunkten $t_{i}$ und den komplexen Messwerten $x_{i}$ und $y_{i}$ (Rechteckimpulse)		T=100.0 dr=1.0 t=[] x=[] y=[] te=0.0 while te<T: tp=exponential(1.0/dr) te+=tp if te<T: t.append(te) xe=0.0+0j if (te>=10) and (te<70): xe+=1.0 if (te>=30) and (te<90): xe+=1.0j ye=0.0+0j if (te>=20) and (te<90): ye+=1.0 if (te>=10) and (te<90): ye+=1.0j x.append(xe) y.append(ye) t=array(t) x=array(x) y=array(y) plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o') show()	`T=100; dr=1; t=[]; x=[]; y=[]; te=0; while te<T tp=-log(1-rand())/dr; te=te+tp; if te<T t(end+1)=te; xe=0+0j; if (te>=10) && (te<70) xe=xe+1; end if (te>=30) && (te<90) xe=xe+1j; end ye=0+0j; if (te>=20) && (te<90) ye=ye+1; end if (te>=10) && (te<90) ye=ye+1j; end x(end+1)=xe; y(end+1)=ye; end end plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')`
Generieren von reellen Gewichten $g_{i}$ passend zu den $(t_{i}, x_{i}, y_{i})$ (Kosinusfenster als Beispiel)		`g=sin(pi*(t/float(T))) plot(t,g,'o') show()`	`g=sin(pi*(t/T)); plot(t,g,'o')`
(Momente verschwinden für Energiesignale)
Kreuzkorrelationsfunktion und Kreuzenergiedichtespektrum über Slotkorrelation (ohne koinzidente Produkte)	$R_{E, y x, k} = R_{E, y x} (τ_{k}) = R_{E, x y}^{} (- τ_{k}) = (T - \| τ_{k} \|) \frac{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j} x_{i}^{} y_{j}}{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} b_{k} (t_{j} - t_{i}) g_{i} g_{j}}$ $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_{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 * 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)+0j R0=zeros(K) i=0 jb=i while i<N: j=jb while (jb<N) and (t[jb]<t[i]+((K-1)//2+0.5)dt): jb+=1 j=jb-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]+=g[i]g[j]conj(x[i])y[j] R0[k]+=g[i]g[j] j-=1 i+=1 REyx=zeros(K)+0j for k in range(-(K//2),(K+1)//2): if R0[k]!=0: REyx[k]=R1[k]/float(R0[k])(T-abs(tau[k])) plot(tau,real(REyx),'o',tau,imag(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()	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)+0j; R0=zeros(1,K); i=2; jb=i; while i<=N j=jb; while (jb<=N) && (t(jb)<t(i)+(fix((K-1)/2)+0.5)dt) jb=jb+1; end j=jb-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)+g(i)g(j)conj(x(i))y(j); R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+g(i)g(j); j=j-1; end i=i+1; end REyx=zeros(1,K)+0j; for k=1:K if R0(k)~=0 REyx(k)=R1(k)./R0(k).(T-abs(tau(k))); end end plot(tau,real(REyx),'o',tau,imag(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, mit Abzug der koinzidenten Produkte, ohne Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{} y_{i}}{{(\sum_{i = 0}^{N - 1} g_{i})}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} y_{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, 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 = [2 T / Δ 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 * N=len(t) dt=1.0 J=int(ceil(2T/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,N): X+=g[i]x[i]exp(-2jpifpt[i]) Y+=g[i]y[i]exp(-2jpifpt[i]) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/float(sum(g)2-sum(g2)) REyx=ifft(Eyx)/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) REyx=append(REyx[0:K//2+1],REyx[-(K//2):len(REyx)]) plot(tau,real(REyx),'o',tau,imag(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()	N=length(t); dt=1.0; J=ceil(2T/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:N X=X+g(i)x(i)exp(-2jpifpt(i)); Y=Y+g(i)y(i)exp(-2jpifpt(i)); end Eyx=T^2(conj(X).Y-sum(g.^2.conj(x).y))/(sum(g)^2-sum(g.^2)); REyx=ifft(Eyx)/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); REyx=[REyx(1:fix(K/2)+1),REyx(fix(K/2)+3:end)]; plot(tau,real(REyx),'o',tau,imag(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 Abzug der koinzidenten Produkte, mit Normierung)	$E_{y x, j} = E_{y x} (f_{j}) = T^{2} \frac{X_{j}^{} Y_{j} - \sum_{i = 0}^{N - 1} g_{i}^{2} x_{i}^{} y_{i}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $E_{y x, j}^{'} = E_{y x}^{'} (f_{j}) = T^{2} \frac{{\| X_{j}^{'} \|}^{2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}{X_{0}^{' 2} - \sum_{i = 0}^{N - 1} g_{i}^{2}}$ $X_{j} = \sum_{i = 0}^{N - 1} g_{i} x_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $Y_{j} = \sum_{i = 0}^{N - 1} g_{i} y_{i} 𝐞^{- 2 π 𝐢 f_{j} t_{i}}$ $X_{j}^{'} = \sum_{i = 0}^{N - 1} g_{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, 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{T - \| τ_{k} \|}{R_{E, y x, k}^{'}} R_{E, y x, k}$ $τ_{k} = k Δ t$ $k = - ⌊ K / 2 ⌋ \dots ⌊ (K - 1) / 2 ⌋$ $K = [2 T / Δ 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 * N=len(t) dt=1.0 J=int(ceil(2T/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 for i in range(0,N): X+=g[i]x[i]exp(-2jpifpt[i]) Y+=g[i]y[i]exp(-2jpifpt[i]) Xp+=g[i]exp(-2jpifpt[i]) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) 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]=T-abs(tau[k]) plot(tau,real(REyx),'o',tau,imag(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()	N=length(t); dt=1.0; J=ceil(2T/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; for i=1:N X=X+g(i)x(i)exp(-2jpifpt(i)); Y=Y+g(i)y(i)exp(-2jpifpt(i)); Xp=Xp+g(i)exp(-2jpifpt(i)); end Eyx=T^2(conj(X).Y-sum(g.^2.conj(x).y))/(abs(Xp(1))^2-sum(g.^2)); Eyxp=T^2(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2)); REyx=ifft(Eyx)/dt; REyxp=real(ifft(Eyxp))/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); 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)(T-abs(tau(k))); end end plot(tau,real(REyx),'o',tau,imag(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')
(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 (mit Abzug der koinzidenten Produkte, ohne Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x1=zeros(J)+0j y1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x1[j]+=g[i]x[i]; y1[j]+=g[i]y[i]; X=fft(x1) Y=fft(y1) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/float(sum(g)2-sum(g2)) REyx=ifft(Eyx)/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) REyx=append(REyx[0:K//2+1],REyx[-(K//2):len(REyx)]) plot(tau,real(REyx),'o',tau,imag(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()	N=length(t); dt=1.0; J=ceil(2T/dt); x1=zeros([1,J])+0j; y1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x1(j)=x1(j)+g(i)x(i); y1(j)=y1(j)+g(i)y(i); end X=fft(x1); Y=fft(y1); Eyx=T^2(conj(X).Y-sum(g.^2.conj(x).y))/(sum(g)^2-sum(g.^2)); REyx=ifft(Eyx)/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); REyx=[REyx(1:fix(K/2)+1),REyx(fix(K/2)+3:end)]; plot(tau,real(REyx),'o',tau,imag(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 Zeitquantisierung (mit Abzug der koinzidenten Produkte, mit Normierung)		from numpy.fft import * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) x0=zeros(J) x1=zeros(J)+0j y1=zeros(J)+0j for i in range(0,N): j=int(floor((t[i]-Tfloor(t[i]/float(T)))/float(dt))) x0[j]+=g[i]; x1[j]+=g[i]x[i]; y1[j]+=g[i]y[i]; X=fft(x1) Y=fft(y1) Xp=fft(x0) Eyx=T*2(conj(X)Y-sum(g2conj(x)y))/(abs(Xp[0])2-sum(g2)) Eyxp=T2(abs(Xp)2-sum(g2))/(abs(Xp[0])2-sum(g2)) REyx=ifft(Eyx)/float(dt) REyxp=real(ifft(Eyxp))/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) 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]=T-abs(tau[k]) plot(tau,real(REyx),'o',tau,imag(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()	N=length(t); dt=1.0; J=ceil(2T/dt); x0=zeros([1,J]); x1=zeros([1,J])+0j; y1=zeros([1,J])+0j; for i=1:N j=fix((t(i)-Tfix(t(i)/T))/dt)+1; x0(j)=x0(j)+g(i); x1(j)=x1(j)+g(i)x(i); y1(j)=y1(j)+g(i)y(i); end X=fft(x1); Y=fft(y1); Xp=fft(x0); Eyx=T^2(conj(X).Y-sum(g.^2.conj(x).y))/(abs(Xp(1))^2-sum(g.^2)); Eyxp=T^2(abs(Xp).^2-sum(g.^2))/(abs(Xp(1))^2-sum(g.^2)); REyx=ifft(Eyx)/dt; REyxp=real(ifft(Eyxp))/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); 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)(T-abs(tau(k))); end end plot(tau,real(REyx),'o',tau,imag(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 * N=len(t) dt=1.0 J=int(ceil(2T/float(dt))) xp=zeros(J)+0j yp=zeros(J)+0j 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] yp[j]=y[i] Xp=fft(xp) Yp=fft(yp) Eyxp=dt*2conj(Xp)Yp REyxp=ifft(Eyxp)/float(dt) K=int(round(2T/float(dt)))-1 tau=roll(arange(-(K//2),(K+1)//2)dt,(K+1)//2) REyx=zeros(K)+0j c=exp(-dd/float(dt))/(1-exp(-dd/float(dt)))2 for k in range(-(K//2),(K+1)//2): if k==0: REyx[k]=REyxp[k] else: REyx[k]=(2c+1)REyxp[k]-cREyxp[k-1]-cREyxp[k+1] plot(tau,real(REyx),'o',tau,imag(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()	N=length(t); dt=1.0; J=ceil(2T/dt); xp=zeros([1,J])+0j; yp=zeros([1,J])+0j; 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); yp(j)=y(i); end end Xp=fft(xp); Yp=fft(yp); Eyxp=dt^2conj(Xp).Yp; REyxp=ifft(Eyxp)/dt; K=round(2T/dt)-1; tau=circshift((-fix(K/2):fix((K-1)/2))dt,[0;fix((K+1)/2)]); REyx=zeros([1,K])+0j; c=exp(-dd/dt)/(1-exp(-dd/dt))^2; for k=-fix(K/2):fix((K-1)/2) if k==0 REyx(mod(k+K,K)+1)=REyxp(mod(k+length(REyxp),length(REyxp))+1); else REyx(mod(k+K,K)+1)=(2c+1)REyxp(mod(k+length(REyxp),length(REyxp))+1)-cREyxp(mod(k-1+length(REyxp),length(REyxp))+1)-cREyxp(mod(k+1+length(REyxp),length(REyxp))+1); end end plot(tau,real(REyx),'o',tau,imag(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')