Signal- und Messdatenverarbeitung

Codes für komplexe Einzelenergiesignale, ohne 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 xi mit i=0…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')
(Momente verschwinden für Energiesignale)
(Auto-)Korrelationsfunktion RE,k=RE(τk)=Δt⋅∑i=max(0,−k)N−1−max(0,k)xi*xi+k
τk=kΔt
k=−(N−1)…N−1
(außerhalb dieses Bereichs null) 		direkte Berechnung:
RE=zeros(2*N-1)+0j
tau=zeros(2*N-1)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0+0j
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=conj(x[i])*x[i+k]
  
  RE[k]=dt*sum1

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(x,zeros(N)))
E=dt**2*abs(X)**2
RE=ifft(E)/float(dt)
RE=append(RE[0:N],RE[N+1:2*N])
plot(tau,real(RE),'o',tau,imag(RE),'o')
show()
 direkte Berechnung:
RE=zeros(1,2*N-1)+0j;
tau=zeros(1,2*N-1);
for k=-(N-1):N-1
tau(mod(k,2*N-1)+1)=k*dt;
sum1=0+0j;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+conj(x(i))*x(i+k);
end
RE(mod(k,2*N-1)+1)=dt*sum1;
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([x,zeros(1,N)]);
E=dt^2*abs(X).^2;
RE=ifft(E)/dt;
RE=[RE(1:N),RE(N+2:end)];
plot(tau,real(RE),'o',tau,imag(RE),'o')
(Auto-)Energiedichtespektrum Xj=FFT{xi|0}=∑i=0N−1xi𝐞−π𝐢ij/N
(imaginäre Einheit 𝐢, implizites Zeropadding mit N Werten ergibt 2N Frequenzen)
Ej=E(fj)=(Δt)2Xj*Xj=(Δt)2|Xj|2
fj=j2NΔt
j=−N…N−1 		direkte Berechnung:
from numpy.fft import *
f=roll(arange(-N,N)/float(2*N*dt),N)
X=fft(append(x,zeros(N)))
E=dt**2*abs(X)**2
plot(f,E,'o')
show()
Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
RE=zeros(2*N)+0j
for k in range(-N,N):
  sum1=0+0j
  for i in range(maximum(0,-k),N-maximum(0,k)):
    sum1+=conj(x[i])*x[i+k]
  
  RE[k]=dt*sum1

f=roll(arange(-N,N)/float(2*N*dt),N)
E=dt*real(fft(RE))
plot(f,E,'o')
show()
 direkte Berechnung:
f=circshift((-N:N-1)/(2*N*dt),[0;N]);
X=fft([x,zeros(1,N)]);
E=dt^2*abs(X).^2;
plot(f,E,'o')
Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
RE=zeros(2*N,1)+0j;
for k=-N:N-1
sum1=0+0j;
for i=1+max(0,-k):N-max(0,k)
sum1=sum1+conj(x(i))*x(i+k);
end
RE(mod(k,2*N)+1)=dt*sum1;
end
f=circshift((-N:N-1)/(2*N*dt),[0;N]);
E=dt*real(fft(RE));
plot(f,E,'o')