Signal- und Messdatenverarbeitung

Codes für reelle 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 Werten xi mit i=0…N−1 (Rechteckimpuls als Beispiel)
N=100
dt=1.0
t=arange(0,N)*dt
x=zeros(N)
for i in range(10,90):
  x[i]=1

plot(t,x,'o')
show()
N=100;
dt=1.0;
t=(0:N-1)*dt;
x=zeros(1,N);
for i=10:89
x(i)=1;
end
plot(t,x,'o')
(Momente verschwinden für Energiesignale)
(Auto-)Korrelationsfunktion RE,k=RE(τk)=Δt⋅∑i=0N−1−|k|xixi+|k|
τk=kΔt
k=−(N−1)…N−1
(außerhalb dieses Bereichs null) 		direkte Berechnung:
RE=zeros(2*N-1)
tau=zeros(2*N-1)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0
  for i in range(0,N-abs(k)):
    sum1+=x[i]*x[i+abs(k)]
  
  RE[k]=dt*sum1

plot(tau,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=real(ifft(E))/float(dt)
RE=append(RE[0:N],RE[N+1:2*N])
plot(tau,RE,'o')
show()
 direkte Berechnung:
RE=zeros(1,2*N-1);
tau=zeros(1,2*N-1);
for k=-(N-1):N-1
tau(mod(k,2*N-1)+1)=k*dt;
sum1=0;
for i=1:N-abs(k)
sum1=sum1+x(i)*x(i+abs(k));
end
RE(mod(k,2*N-1)+1)=dt*sum1;
end
plot(tau,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=real(ifft(E))/dt;
RE=[RE(1:N),RE(N+2:end)];
plot(tau,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)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)
for k in range(-N,N):
  sum1=0
  for i in range(0,N-abs(k)):
    sum1+=x[i]*x[i+abs(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);
for k=-N:N-1
sum1=0;
for i=1:N-abs(k)
sum1=sum1+x(i)*x(i+abs(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')