Signal- und Messdatenverarbeitung

Testcodes für Kapitel 10: Spektren

systematische Zusammenstellung von kombinierbaren Programmabschnitten zur Signal- und Messdatenverarbeitung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#pip install --upgrade numpy
#pip install --upgrade matplotlib
from numpy import *
from numpy.random import *
from numpy.fft import *
from matplotlib.pyplot import *
%Matlab: erfordert die Statistics and Machine Learning Toolbox
%Octave: pkg load statistics
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
end
periodisches Signal (Überlagerung zweier harmonischer Signale)
N=100
dt=1.0
t=arange(0,N)*dt
x=sin(2*pi*arange(0,N)/20.0)+sin(4*pi*arange(0,N)/20.0+0.6)
plot(t,x,'.')
show()
N=100;
dt=1.0;
t=(0:N-1)*dt;
x=sin(2*pi*(0:N-1)/20)+sin(4*pi*(0:N-1)/20+0.6);
plot(t,x,'.')
Leistungsspektrum
X=fft(x)
df=1/float(N*dt)
f=((arange(0,N)+N//2)%N-(N//2))*df
P=abs(X)**2/float(N**2)
plot(f,P,'.')
show()
X=fft(x);
df=1/(N*dt);
f=(mod((0:N-1)+fix(N/2),N)-fix(N/2))*df;
P=abs(X).^2/N^2;
plot(f,P,'.')
Prüfen anhand der Signalleistung
print(sum(x**2)/float(N))
print(var(x)+mean(x)**2)
print(sum(P))
sum(x.^2)/N
var(x,1)+mean(x)^2
sum(P)
Autokorrelationsfunktion aus Spektrum
R=N*real(ifft(P))
tau=arange(0,N)*dt
plot(tau,R,'.')
show()
R=N*real(ifft(P));
tau=(0:N-1)*dt;
plot(tau,R,'.')
Prüfen anhand der Signalleistung
print(sum(x**2)/float(N))
print(var(x)+mean(x)**2)
print(R[0])
sum(x.^2)/N
var(x,1)+mean(x)^2
R(1)
Energiesignal (Rechteckimpuls)
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,'.')
show()
N=100;
dt=1.0;
t=(0:N-1)*dt;
x=zeros(1,N);
for i=11:90
x(i)=1;
end
plot(t,x,'.')
Energiedichtespektrum via DFT mit Interpretation als periodisch fortgesetztes Signal und anschließende Bestimmung von Zwischenwerten mittels Fourier-Transformation
X=fft(x)
P=abs(X)**2/float(N**2)
Eext=zeros(2*N)
for j in range(0,N):
  Eext[2*j]=(N*dt)**2*P[j]

df=1/float(2*N*dt)
f=((arange(0,2*N)+N)%(2*N)-N)*df
for j in range(0,N):
  Xd=0
  for jj in range(0,N):
    Xd+=x[jj]*exp(-2*pi*1j*f[2*j+1]*dt*jj)
  
  Eext[2*j+1]=dt**2*abs(Xd)**2

plot(f,Eext,'.')
show()
X=fft(x);
P=abs(X).^2/N^2;
Eext=zeros(1,2*N);
for j=1:N
Eext(2*j-1)=(N*dt)^2*P(j);
end
df=1/(2*N*dt);
f=(mod((0:2*N-1)+N,2*N)-N)*df;
for j=1:N
Xd=0;
for jj=1:N
Xd=Xd+x(jj)*exp(-2*pi*1j*f(2*j)*dt*jj);
end
Eext(2*j)=dt**2*abs(Xd)^2;
end
plot(f,Eext,'.')
Energiedichtespektrum via DFT mit Interpretation als periodisch fortgesetztes Signal und anschließende Bestimmung von Zwischenwerten mittels Interpolation der DFT
X=fft(x)
P=abs(X)**2/float(N**2)
df=1/float(N*dt)
f=((arange(0,N)+N//2)%N-(N//2))*df
#Mutiplikation mit rotierendem Zeiger / Verschiebung des Signals um -T/2
Xp=X*exp(pi*1j*f*N*dt)
Eint=zeros(2*N)
for j in range(0,N):
  Eint[2*j]=(N*dt)**2*P[j]

for j in range(-(N//2),(N+1)//2):
  Xf=0
  for jj in range(-(N//2),(N+1)//2):
    Xf+=Xp[jj]*sinc(j+0.5-jj)
  
  Pint=abs(Xf)**2/float(N**2)
  Eint[2*j+1]=(N*dt)**2*Pint

df=1/float(2*N*dt)
f=((arange(0,2*N)+N)%(2*N)-N)*df
plot(f,Eint,'.')
show()
X=fft(x);
P=abs(X).^2/N^2;
df=1/(N*dt);
f=(mod((0:N-1)+fix(N/2),N)-fix(N/2))*df;
%Mutiplikation mit rotierendem Zeiger / Verschiebung des Signals um -T/2
Xp=X.*exp(pi*1j*f*N*dt);
Eint=zeros(1,2*N);
for j=1:N
Eint(2*j-1)=(N*dt)^2*P(j);
end
for j=-floor(N/2):floor((N-1)/2)
Xf=0;
for jj=-floor(N/2):floor((N-1)/2)
Xf=Xf+Xp(mod(jj,N)+1)*sinc(j+0.5-jj);
end
Pint=abs(Xf)^2/N^2;
Eint(mod(2*j+1,2*N)+1)=(N*dt)^2*Pint;
end
df=1/(2*N*dt);
f=(mod((0:2*N-1)+N,2*N)-N)*df;
plot(f,Eint,'.')
Energiedichtespektrum via Zero Padding und Vergleich der Resultate
X=fft(append(x,zeros(N)))
df=1/float(2*N*dt)
f=((arange(0,2*N)+N)%(2*N)-N)*df
E=abs(X)**2*dt**2
plot(f,Eext,'o',f,Eint,'x',f,E,'.')
show()
X=fft([x,zeros(1,N)]);
df=1/(2*N*dt);
f=(mod((0:2*N-1)+N,2*N)-N)*df;
E=abs(X).^2*dt^2;
plot(f,Eext,'o',f,Eint,'x',f,E,'.')
Prüfen anhand der Signalenergie
print(sum(x**2)*dt)
print(sum(E)*df)
sum(x.^2)*dt
sum(E)*df
(Energie-)Autokorrelationsfunktion aus Spektrum
RE=real(ifft(E))/float(dt)
tau=((arange(0,2*N)+N)%(2*N)-N)*dt
plot(tau,RE,'.')
show()
RE=real(ifft(E))/dt;
tau=(mod((0:2*N-1)+N,2*N)-N)*dt;
plot(tau,RE,'.')
Prüfen anhand der Signalenergie
print(sum(x**2)*dt)
print(RE[0])
sum(x.^2)*dt
RE(1)
stationäres Leistungssignal (AR1-Prozess)
N=1000
phi=0.95
ex=1
vx=1
theta=sqrt((1-phi**2)*vx)
x=zeros(N)
x[0]=normal(ex,sqrt(vx))
for i in range(1,N):
  x[i]=phi*(x[i-1]-ex)+normal(ex,theta)

dt=1.0
t=arange(0,N)*dt
plot(t,x,'.')
show()
N=1000;
phi=0.95;
ex=1;
vx=1;
theta=sqrt((1-phi^2)*vx);
x=zeros(1,N);
x(1)=normrnd(ex,sqrt(vx));
for i=2:N
x(i)=phi*(x(i-1)-ex)+normrnd(ex,theta);
end
dt=1.0;
t=(0:N-1)*dt;
plot(t,x,'.')
Leistungsdichtespektrum
X=fft(append(x,zeros(N)))
df=1/float(2*N*dt)
f=((arange(0,2*N)+N)%(2*N)-N)*df
E=abs(X)**2*dt**2
RE=real(ifft(E))/float(dt)
tau=((arange(0,2*N)+N)%(2*N)-N)*dt
R=zeros(len(RE))
for k in range(0,len(RE)):
  if N*dt-abs(tau[k])>0:
    R[k]=RE[k]/(N*dt-abs(tau[k]))
  

S=real(fft(R))*dt
plot(f,S,'.')
show()
X=fft([x,zeros(1,N)]);
df=1/(2*N*dt);
f=(mod((0:2*N-1)+N,2*N)-N)*df;
E=abs(X).^2*dt^2;
RE=real(ifft(E))/dt;
tau=(mod((0:2*N-1)+N,2*N)-N)*dt;
R=zeros(1,length(RE));
for k=1:length(RE)
if N*dt-abs(tau(k))>0
R(k)=RE(k)/(N*dt-abs(tau(k)));
end
end
S=real(fft(R))*dt;
plot(f,S,'.')
Prüfen anhand der Signalleistung
print(sum(x**2)/float(N))
print(var(x)+mean(x)**2)
print(sum(S)*df)
sum(x.^2)/N
var(x,1)+mean(x)^2
sum(S)*df
Autokorrelationsfunktion aus Spektrum
R=real(ifft(S))/float(dt)
tau=((arange(0,2*N)+N)%(2*N)-N)*dt
plot(tau,R,'.')
show()
R=real(ifft(S))/dt;
tau=(mod((0:2*N-1)+N,2*N)-N)*dt;
plot(tau,R,'.')
Prüfen anhand der Signalleistung
print(sum(x**2)/float(N))
print(var(x)+mean(x)**2)
print(R[0])
sum(x.^2)/N
var(x,1)+mean(x)^2
R(1)
Kürzen der Autokorrelationsfunktion auf N′′<2N−1 Werte
Npp=100
Rpp=zeros(Npp)
taupp=roll(arange(-(Npp//2),(Npp+1)//2),(Npp+1)//2)*dt
for k in range(-(Npp//2),(Npp+1)//2):
  Rpp[k]=R[k]

plot(taupp,Rpp,'.')
show()
fpp=roll(arange(-(Npp//2),(Npp+1)//2)/float(Npp*dt),(Npp+1)//2)
Spp=real(fft(Rpp))*dt
plot(fpp,Spp,'.')
show()
Npp=100;
Rpp=zeros(1,Npp);
taupp=circshift((-fix(Npp/2):fix((Npp-1)/2)),[0;fix((Npp+1)/2)])*dt;
for k=-floor(Npp/2):floor((Npp-1)/2)
Rpp(mod(k,Npp)+1)=R(mod(k,length(R))+1);
end;
plot(taupp,Rpp,'.')
waitforbuttonpress
fpp=circshift((-fix(Npp/2):fix((Npp-1)/2))/(Npp*dt),[0;fix((Npp+1)/2)]);
Spp=real(fft(Rpp))*dt;
plot(fpp,Spp,'.')