Signal- und Messdatenverarbeitung

Codes für komplexe Energiesignalpaare, ohne Gewichtung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren von Nx komplexen Werten xi mit i=0…Nx−1 und Ny komplexen Werten yi mit i=0…Ny−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')
(Momente verschwinden für Energiesignale)
Kreuzkorrelationsfunktion RE,yx,k=RE,yx(τk)=RE,xy*(−τk)=Δt⋅∑i=max(0,−k)min(Nx−1,Ny−1−k)xi*yi+k
τk=kΔt
k=−(Nx−1)…Ny−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]=k*dt
  sum1=0+0j
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i])*y[i+k]
  
  REyx[k]=dt*sum1

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)
X=fft(append(x,zeros(Ny)))
Y=fft(append(y,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
REyx=ifft(Eyx)/float(dt)
REyx=append(REyx[0:Ny],REyx[Ny+1:Nx+Ny])
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)=k*dt;
sum1=0+0j;
for i=1+max(0,-k):min(Nx-1,Ny-1-k)
sum1=sum1+conj(x(i))*y(i+k);
end
REyx(mod(k,Nx+Ny-1)+1)=dt*sum1;
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]);
X=fft([x,zeros(1,Ny)]);
Y=fft([y,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
REyx=ifft(Eyx)/dt;
REyx=[REyx(1:Ny-1),REyx(Ny+1:end)];
plot(tau,real(REyx),'o',tau,imag(REyx),'o')
Kreuzenergiedichtespektrum Xj=FFT{xi|0}=∑i=0Nx−1xi𝐞−2π𝐢ij/(Nx+Ny)
Yj=FFT{yi|0}=∑i=0Ny−1yi𝐞−2π𝐢ij/(Nx+Ny)
(imaginäre Einheit 𝐢, implizites Zeropadding mit Ny bzw. Nx Werten ergibt Nx+Ny Frequenzen)
Eyx,j=Eyx(fj)=(Δt)2Xj*Yj
fj=j(Nx+Ny)Δt
j=−⌊(Nx+Ny)/2⌋…⌊(Nx+Ny−1)/2⌋ 		direkte Berechnung:
from numpy.fft import *
f=roll(arange(-((Nx+Ny)//2),(Nx+Ny+1)//2)/float((Nx+Ny)*dt),(Nx+Ny+1)//2)
X=fft(append(x,zeros(Ny)))
Y=fft(append(y,zeros(Nx)))
Eyx=dt**2*conj(X)*Y
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
  for i in range(maximum(0,-k),minimum(Nx,Ny-k)):
    sum1+=conj(x[i])*y[i+k]
  
  REyx[k]=dt*sum1

f=roll(arange(-((Nx+Ny)//2),(Nx+Ny+1)//2)/float((Nx+Ny)*dt),(Nx+Ny+1)//2)
Eyx=dt*fft(REyx)
plot(f,real(Eyx),'o',f,imag(Eyx),'o')
show()
 direkte Berechnung:
f=circshift((-fix((Nx+Ny)/2):fix((Nx+Ny-1)/2))/((Nx+Ny)*dt),[0;fix((Nx+Ny+1)/2)]);
X=fft([x,zeros(1,Ny)]);
Y=fft([y,zeros(1,Nx)]);
Eyx=dt^2*conj(X).*Y;
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;
for i=1+max(0,-k):min(Nx-1,Ny-1-k)
sum1=sum1+conj(x(i))*y(i+k);
end
REyx(mod(k,Nx+Ny)+1)=dt*sum1;
end
f=circshift((-fix((Nx+Ny)/2):fix((Nx+Ny-1)/2))/((Nx+Ny)*dt),[0;fix((Nx+Ny+1)/2)]);
Eyx=dt*fft(REyx);
plot(f,real(Eyx),'o',f,imag(Eyx),'o')