Signal- und Messdatenverarbeitung

Testcodes für Kapitel 11: Fensterfunktionen

weitere Codes zu Fensterfunktionen

systematische Zusammenstellung von kombinierbaren Programmabschnitten zur Signal- und Messdatenverarbeitung

Fensterfunktionen und ihre Spektren (Energiedichtespektren von Energiesignalen)

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.fft import *
from matplotlib.pyplot import *
verschiedene Fensterfunktionen …
N=1000
dt=1.0
t=arange(0,N)*dt
M=800
w1=ones(M)
w2=1-abs((2*arange(0,M)+1)/float(M)-1) # Bartlett
w3=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
w1=append(zeros((N-M)//2),append(w1,zeros((N-M+1)//2)))
w2=append(zeros((N-M)//2),append(w2,zeros((N-M+1)//2)))
w3=append(zeros((N-M)//2),append(w3,zeros((N-M+1)//2)))
plot(t,w1,t,w2,t,w3)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
M=800;
w1=ones(1,M);
w2=1-abs((2*(0:M-1)+1)/M-1); % Bartlett
w3=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
w1=[zeros(1,fix((N-M)/2)),w1,zeros(1,fix((N-M+1)/2))];
w2=[zeros(1,fix((N-M)/2)),w2,zeros(1,fix((N-M+1)/2))];
w3=[zeros(1,fix((N-M)/2)),w3,zeros(1,fix((N-M+1)/2))];
plot(t,w1,t,w2,t,w3)
… und ihre Spektren (kleinere Fenster, Energiedichtespektrum eines Energiesignals, unterschiedlich breite Hauptkeule/Verschmierung und unterschiedlich hohe Nebenzipfel)
N=1000
dt=1.0
t=arange(0,N)*dt
M=20
w1=ones(M)
w2=1-abs((2*arange(0,M)+1)/float(M)-1) # Bartlett
w3=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
w1=append(zeros((N-M)//2),append(w1,zeros((N-M+1)//2)))
w2=append(zeros((N-M)//2),append(w2,zeros((N-M+1)//2)))
w3=append(zeros((N-M)//2),append(w3,zeros((N-M+1)//2)))
plot(t,w1,t,w2,t,w3)
show()
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
W1=fft(append(w1,zeros(N-1)))
W2=fft(append(w2,zeros(N-1)))
W3=fft(append(w3,zeros(N-1)))
E1=dt**2*abs(W1)**2
E2=dt**2*abs(W2)**2
E3=dt**2*abs(W3)**2
semilogy(roll(f,N-1),roll(E1,N-1),roll(f,N-1),roll(E2,N-1),roll(f,N-1),roll(E3,N-1))
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
M=20;
w1=ones(1,M);
w2=1-abs((2*(0:M-1)+1)/M-1); % Bartlett
w3=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
w1=[zeros(1,fix((N-M)/2)),w1,zeros(1,fix((N-M+1)/2))];
w2=[zeros(1,fix((N-M)/2)),w2,zeros(1,fix((N-M+1)/2))];
w3=[zeros(1,fix((N-M)/2)),w3,zeros(1,fix((N-M+1)/2))];
plot(t,w1,t,w2,t,w3)
waitforbuttonpress
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
W1=fft([w1,zeros(1,N-1)]);
W2=fft([w2,zeros(1,N-1)]);
W3=fft([w3,zeros(1,N-1)]);
E1=dt^2*abs(W1).^2;
E2=dt^2*abs(W2).^2;
E3=dt^2*abs(W3).^2;
semilogy(circshift(f,N-1),circshift(E1,N-1),circshift(f,N-1),circshift(E2,N-1),circshift(f,N-1),circshift(E3,N-1))
verschieden breite Fensterfunktionen (Rechteckfenster) …
N=1000
dt=1.0
t=arange(0,N)*dt
w1=ones(16)
w2=ones(26)
w3=ones(42)
w1=append(zeros((N-len(w1))//2),append(w1,zeros((N-len(w1)+1)//2)))
w2=append(zeros((N-len(w2))//2),append(w2,zeros((N-len(w2)+1)//2)))
w3=append(zeros((N-len(w3))//2),append(w3,zeros((N-len(w3)+1)//2)))
plot(t,w1,t,w2,t,w3)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
w1=ones(1,16);
w2=ones(1,26);
w3=ones(1,42);
w1=[zeros(1,fix((N-length(w1))/2)),w1,zeros(1,fix((N-length(w1)+1)/2))];
w2=[zeros(1,fix((N-length(w2))/2)),w2,zeros(1,fix((N-length(w2)+1)/2))];
w3=[zeros(1,fix((N-length(w3))/2)),w3,zeros(1,fix((N-length(w3)+1)/2))];
plot(t,w1,t,w2,t,w3)
… und ihre Spektren (gemeinsame Hüllkurve, unterschiedlich starkte Verschmierung)
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
W1=fft(append(w1,zeros(N-1)))
W2=fft(append(w2,zeros(N-1)))
W3=fft(append(w3,zeros(N-1)))
E1=dt**2*abs(W1)**2
E2=dt**2*abs(W2)**2
E3=dt**2*abs(W3)**2
semilogy(roll(f,N-1),roll(E1,N-1),roll(f,N-1),roll(E2,N-1),roll(f,N-1),roll(E3,N-1))
show()
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
W1=fft([w1,zeros(1,N-1)]);
W2=fft([w2,zeros(1,N-1)]);
W3=fft([w3,zeros(1,N-1)]);
E1=dt^2*abs(W1).^2;
E2=dt^2*abs(W2).^2;
E3=dt^2*abs(W3).^2;
semilogy(circshift(f,N-1),circshift(E1,N-1),circshift(f,N-1),circshift(E2,N-1),circshift(f,N-1),circshift(E3,N-1))
verschieden breite Fensterfunktionen (Blackman-Fenster) …
N=1000
dt=1.0
t=arange(0,N)*dt
M=32
w1=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M))
M=52
w2=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M))
M=84
w3=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M))
w1=append(zeros((N-len(w1))//2),append(w1,zeros((N-len(w1)+1)//2)))
w2=append(zeros((N-len(w2))//2),append(w2,zeros((N-len(w2)+1)//2)))
w3=append(zeros((N-len(w3))//2),append(w3,zeros((N-len(w3)+1)//2)))
plot(t,w1,t,w2,t,w3)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
M=32;
w1=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M);
M=52;
w2=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M);
M=84;
w3=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M);
w1=[zeros(1,fix((N-length(w1))/2)),w1,zeros(1,fix((N-length(w1)+1)/2))];
w2=[zeros(1,fix((N-length(w2))/2)),w2,zeros(1,fix((N-length(w2)+1)/2))];
w3=[zeros(1,fix((N-length(w3))/2)),w3,zeros(1,fix((N-length(w3)+1)/2))];
plot(t,w1,t,w2,t,w3)
… und ihre Spektren (ähnliches Verhalten, unterschiedlich starkte Verschmierung)
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
W1=fft(append(w1,zeros(N-1)))
W2=fft(append(w2,zeros(N-1)))
W3=fft(append(w3,zeros(N-1)))
E1=dt**2*abs(W1)**2
E2=dt**2*abs(W2)**2
E3=dt**2*abs(W3)**2
semilogy(roll(f,N-1),roll(E1,N-1),roll(f,N-1),roll(E2,N-1),roll(f,N-1),roll(E3,N-1))
show()
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
W1=fft([w1,zeros(1,N-1)]);
W2=fft([w2,zeros(1,N-1)]);
W3=fft([w3,zeros(1,N-1)]);
E1=dt^2*abs(W1).^2;
E2=dt^2*abs(W2).^2;
E3=dt^2*abs(W3).^2;
semilogy(circshift(f,N-1),circshift(E1,N-1),circshift(f,N-1),circshift(E2,N-1),circshift(f,N-1),circshift(E3,N-1))
verschieden breite Fensterfunktionen (Gauß-Fenster) …
N=1000
dt=1.0
t=arange(0,N)*dt
w1=exp(-(2*arange(0,N)-N-1)**2/float(8*10**2))
w2=exp(-(2*arange(0,N)-N-1)**2/float(8*16**2))
w3=exp(-(2*arange(0,N)-N-1)**2/float(8*26**2))
plot(t,w1,t,w2,t,w3)
show()
semilogy(t,w1,t,w2,t,w3)
ylim(1e-9,1)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
w1=exp(-(2*(0:N-1)-N-1).^2/(8*10^2));
w2=exp(-(2*(0:N-1)-N-1).^2/(8*16^2));
w3=exp(-(2*(0:N-1)-N-1).^2/(8*26^2));
plot(t,w1,t,w2,t,w3)
waitforbuttonpress
semilogy(t,w1,t,w2,t,w3)
ylim([1e-9,1])
… und ihre Spektren (Korrespondenz der Gauß-Form zwischen Zeitsignal und seinem Spektrum)
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
W1=fft(append(w1,zeros(N-1)))
W2=fft(append(w2,zeros(N-1)))
W3=fft(append(w3,zeros(N-1)))
E1=dt**2*abs(W1)**2
E2=dt**2*abs(W2)**2
E3=dt**2*abs(W3)**2
plot(roll(f,N-1),roll(E1,N-1),roll(f,N-1),roll(E2,N-1),roll(f,N-1),roll(E3,N-1))
show()
semilogy(roll(f,N-1),roll(E1,N-1),roll(f,N-1),roll(E2,N-1),roll(f,N-1),roll(E3,N-1))
ylim(1e-18,1e6)
show()
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
W1=fft([w1,zeros(1,N-1)]);
W2=fft([w2,zeros(1,N-1)]);
W3=fft([w3,zeros(1,N-1)]);
E1=dt^2*abs(W1).^2;
E2=dt^2*abs(W2).^2;
E3=dt^2*abs(W3).^2;
plot(circshift(f,N-1),circshift(E1,N-1),circshift(f,N-1),circshift(E2,N-1),circshift(f,N-1),circshift(E3,N-1))
waitforbuttonpress
semilogy(circshift(f,N-1),circshift(E1,N-1),circshift(f,N-1),circshift(E2,N-1),circshift(f,N-1),circshift(E3,N-1))
ylim([1e-18,1e6])

Anwendung von Fensterfunktionen

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.fft import *
from matplotlib.pyplot import *
konstantes Signal mit angewendeter Fensterfunktion …
N=1000
dt=1.0
t=arange(0,N)*dt
x=ones(N)
M=128
w=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
wext=append(zeros((N-M)//2),append(w,zeros((N-M+1)//2)))
xw=x*wext
plot(t,x,t,wext,t,xw,'--')
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=ones(1,N);
M=128;
w=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
wext=[zeros(1,fix((N-M)/2)),w,zeros(1,fix((N-M+1)/2))];
xw=x.*wext;
plot(t,x,t,wext,t,xw,'--')
… und die Spektren
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
xext=ones(2*N-1)
X=fft(xext)
W=fft(append(wext,zeros(N-1)))
XW=fft(append(xw,zeros(N-1)))
Ex=dt**2*abs(X)**2
Ew=dt**2*abs(W)**2
Exw=dt**2*abs(XW)**2
plot(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1))
show()
semilogy(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1),'--')
ylim(1e-18,1e6)
show()
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
xext=ones(1,2*N-1);
X=fft(xext);
W=fft([wext,zeros(1,N-1)]);
XW=fft([xw,zeros(1,N-1)]);
Ex=dt^2*abs(X).^2;
Ew=dt^2*abs(W).^2;
Exw=dt^2*abs(XW).^2;
plot(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
waitforbuttonpress
semilogy(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
ylim([1e-18,1e6])
harmonisches Signal mit angewendeter Fensterfunktion …
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(160/float(2*N-1)*pi*arange(0,N))
M=128
w=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
wext=append(zeros((N-M)//2),append(w,zeros((N-M+1)//2)))
xw=x*wext
plot(t,x,t,wext,t,xw,'--')
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(160/(2*N-1)*pi*(0:N-1));
M=128;
w=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
wext=[zeros(1,fix((N-M)/2)),w,zeros(1,fix((N-M+1)/2))];
xw=x.*wext;
plot(t,x,t,wext,t,xw,'--')
… und die Spektren
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
xext=sin(160/float(2*N-1)*pi*arange(0,2*N-1))
X=fft(xext)
W=fft(append(wext,zeros(N-1)))
XW=fft(append(xw,zeros(N-1)))
Ex=dt**2*abs(X)**2
Ew=dt**2*abs(W)**2
Exw=dt**2*abs(XW)**2
plot(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1))
show()
semilogy(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1),'--')
ylim(1e-18,1e6)
show()
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
xext=sin(160/(2*N-1)*pi*(0:2*N-2));
X=fft(xext);
W=fft([wext,zeros(1,N-1)]);
XW=fft([xw,zeros(1,N-1)]);
Ex=dt^2*abs(X).^2;
Ew=dt^2*abs(W).^2;
Exw=dt^2*abs(XW).^2;
plot(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
waitforbuttonpress
semilogy(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
ylim([1e-18,1e6])
harmonisches Signal mit Gleichanteil und angewendeter Fensterfunktion …
N=1000
dt=1.0
t=arange(0,N)*dt
x=0.5*(1+sin(160/float(2*N-1)*pi*arange(0,N)))
M=128
w=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
wext=append(zeros((N-M)//2),append(w,zeros((N-M+1)//2)))
xw=x*wext
plot(t,x,t,wext,t,xw,'--')
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=0.5*(1+sin(160/(2*N-1)*pi*(0:N-1)));
M=128;
w=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
wext=[zeros(1,fix((N-M)/2)),w,zeros(1,fix((N-M+1)/2))];
xw=x.*wext;
plot(t,x,t,wext,t,xw,'--')
… und die Spektren
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
xext=0.5*(1+sin(160/float(2*N-1)*pi*arange(0,2*N-1)))
X=fft(xext)
W=fft(append(wext,zeros(N-1)))
XW=fft(append(xw,zeros(N-1)))
Ex=dt**2*abs(X)**2
Ew=dt**2*abs(W)**2
Exw=dt**2*abs(XW)**2
plot(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1))
show()
semilogy(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1),'--')
ylim(1e-18,1e6)
show()
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
xext=0.5*(1+sin(160/(2*N-1)*pi*(0:2*N-2)));
X=fft(xext);
W=fft([wext,zeros(1,N-1)]);
XW=fft([xw,zeros(1,N-1)]);
Ex=dt^2*abs(X).^2;
Ew=dt^2*abs(W).^2;
Exw=dt^2*abs(XW).^2;
plot(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
waitforbuttonpress
semilogy(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
ylim([1e-18,1e6])
Superposition zweier harmonischer Signale mit angewendeter Fensterfunktion …
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(90/float(2*N-1)*pi*arange(0,N))+sin(240/float(2*N-1)*pi*arange(0,N))
M=128
w=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
wext=append(zeros((N-M)//2),append(w,zeros((N-M+1)//2)))
xw=x*wext
plot(t,x,t,wext,t,xw,'--')
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(90/(2*N-1)*pi*(0:N-1))+sin(240/(2*N-1)*pi*(0:N-1));
M=128;
w=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
wext=[zeros(1,fix((N-M)/2)),w,zeros(1,fix((N-M+1)/2))];
xw=x.*wext;
plot(t,x,t,wext,t,xw,'--')
… und die Spektren
f=roll(arange(-(N-1),N)/float((2*N-1)*dt),N)
xext=sin(90/float(2*N-1)*pi*arange(0,2*N-1))+sin(240/float(2*N-1)*pi*arange(0,2*N-1))
X=fft(xext)
W=fft(append(wext,zeros(N-1)))
XW=fft(append(xw,zeros(N-1)))
Ex=dt**2*abs(X)**2
Ew=dt**2*abs(W)**2
Exw=dt**2*abs(XW)**2
plot(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1))
show()
semilogy(roll(f,N-1),roll(Ex,N-1),roll(f,N-1),roll(Ew,N-1),roll(f,N-1),roll(Exw,N-1),'--')
ylim(1e-18,1e6)
show()
f=circshift((-(N-1):N-1)/((2*N-1)*dt),N);
xext=sin(90/(2*N-1)*pi*(0:2*N-2))+sin(240/(2*N-1)*pi*(0:2*N-2));
X=fft(xext);
W=fft([wext,zeros(1,N-1)]);
XW=fft([xw,zeros(1,N-1)]);
Ex=dt^2*abs(X).^2;
Ew=dt^2*abs(W).^2;
Exw=dt^2*abs(XW).^2;
plot(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
waitforbuttonpress
semilogy(circshift(f,N-1),circshift(Ex,N-1),circshift(f,N-1),circshift(Ew,N-1),circshift(f,N-1),circshift(Exw,N-1),'--')
ylim([1e-18,1e6])

Missbrauch von Fensterfunktionen zur Unterdrückung von Signalsprüngen

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
Signalausschnitt (ohne Modulation der Amplitude) eines periodischen Signals mit vollen Perioden im Beobachtungsfenster, …
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(2*pi*t*28/float(N))+normal(0,0.1,N)
M=250
w=ones(M)
tw=arange(0,M)*dt
xw=x[0:M]*w
plot(t,x,tw,xw)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(2*pi*t*28/N)+normrnd(0,0.1,[1,N]);
M=250;
w=ones(1,M);
tw=(0:M-1)*dt;
xw=x(1:M).*w;
plot(t,x,tw,xw)
… Leistungsspektrum eines periodischen Signals (ähnlich dem Spektrum des langen Signals, nur geringere spektrale Auflösung) …
f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
X=fft(x)
P=abs(X)**2/N**2
fw=roll(arange(-(M//2),M-(M//2))/float(M*dt),M-(M//2))
Xw=fft(xw)
Pw=abs(Xw)**2/M**2
semilogy(roll(f,N//2),roll(P,N//2),roll(fw,M//2),roll(Pw,M//2))
show()
f=circshift((-fix(N/2):N-fix(N/2)-1)/(N*dt),N-fix(N/2));
X=fft(x);
P=abs(X).^2/N^2;
fw=circshift((-fix(M/2):M-fix(M/2)-1)/(M*dt),M-fix(M/2));
Xw=fft(xw);
Pw=abs(Xw).^2/M^2;
semilogy(circshift(f,fix(N/2)),circshift(P,fix(N/2)),circshift(fw,fix(M/2)),circshift(Pw,fix(M/2)))
… und Korrelationsfunktion (wie die des langen Signals, nur kürzer)
tau=arange(0,N)*dt
R=N*real(ifft(P))
tauw=arange(0,M)*dt
Rw=M*real(ifft(Pw))
plot(tau,R,tauw,Rw)
show()
tau=(0:N-1)*dt;
R=N*real(ifft(P));
tauw=(0:M-1)*dt;
Rw=M*real(ifft(Pw));
plot(tau,R,tauw,Rw)
Signalausschnitt mit Fensterfunktion (mit Modulation der Amplitude) eines periodischen Signals mit vollen Perioden im Beobachtungsfenster, …
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(2*pi*t*28/float(N))+normal(0,0.1,N)
M=250
w=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
tw=arange(0,M)*dt
xw=x[0:M]*w
plot(t,x,tw,xw)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(2*pi*t*28/N)+normrnd(0,0.1,[1,N]);
M=250;
w=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
tw=(0:M-1)*dt;
xw=x(1:M).*w;
plot(t,x,tw,xw)
… Leistungsspektrum eines periodischen Signals (verschmiert) …
f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
X=fft(x)
P=abs(X)**2/N**2
fw=roll(arange(-(M//2),M-(M//2))/float(M*dt),M-(M//2))
Xw=fft(xw)
Pw=abs(Xw)**2/M**2
semilogy(roll(f,N//2),roll(P,N//2),roll(fw,M//2),roll(Pw,M//2))
show()
f=circshift((-fix(N/2):N-fix(N/2)-1)/(N*dt),N-fix(N/2));
X=fft(x);
P=abs(X).^2/N^2;
fw=circshift((-fix(M/2):M-fix(M/2)-1)/(M*dt),M-fix(M/2));
Xw=fft(xw);
Pw=abs(Xw).^2/M^2;
semilogy(circshift(f,fix(N/2)),circshift(P,fix(N/2)),circshift(fw,fix(M/2)),circshift(Pw,fix(M/2)))
… und Korrelationsfunktion (ähnlich der des langen Signals/phasenrichtig, aber zusätzliche Modulation)
tau=arange(0,N)*dt
R=N*real(ifft(P))
tauw=arange(0,M)*dt
Rw=M*real(ifft(Pw))
plot(tau,R,tauw,Rw)
show()
tau=(0:N-1)*dt;
R=N*real(ifft(P));
tauw=(0:M-1)*dt;
Rw=M*real(ifft(Pw));
plot(tau,R,tauw,Rw)
Versuch der Korrektur durch leistungsangepasste Amplitude der Fensterfunktion
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(2*pi*t*28/float(N))+normal(0,0.1,N)
M=250
w=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
w*=sqrt(M/sum(w**2)) # Leistungsanpassung des Fensters
tw=arange(0,M)*dt
xw=x[0:M]*w
plot(t,x,tw,xw)
show()
f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
X=fft(x)
P=abs(X)**2/N**2
fw=roll(arange(-(M//2),M-(M//2))/float(M*dt),M-(M//2))
Xw=fft(xw)
Pw=abs(Xw)**2/M**2
semilogy(roll(f,N//2),roll(P,N//2),roll(fw,M//2),roll(Pw,M//2))
show() # Verschmierung des Spektrums bleibt
tau=arange(0,N)*dt
R=N*real(ifft(P))
tauw=arange(0,M)*dt
Rw=M*real(ifft(Pw))
plot(tau,R,tauw,Rw)
show() # Leistungsanpassung bei tau=0, Modulation bleibt
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(2*pi*t*28/N)+normrnd(0,0.1,[1,N]);
M=250;
w=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
w=w*sqrt(M/sum(w.^2)); % Leistungsanpassung des Fensters
tw=(0:M-1)*dt;
xw=x(1:M).*w;
plot(t,x,tw,xw)
waitforbuttonpress
f=circshift((-fix(N/2):N-fix(N/2)-1)/(N*dt),N-fix(N/2));
X=fft(x);
P=abs(X).^2/N^2;
fw=circshift((-fix(M/2):M-fix(M/2)-1)/(M*dt),M-fix(M/2));
Xw=fft(xw);
Pw=abs(Xw).^2/M^2;
semilogy(circshift(f,fix(N/2)),circshift(P,fix(N/2)),circshift(fw,fix(M/2)),circshift(Pw,fix(M/2))) % Verschmierung des Spektrums bleibt
waitforbuttonpress
tau=(0:N-1)*dt;
R=N*real(ifft(P));
tauw=(0:M-1)*dt;
Rw=M*real(ifft(Pw));
plot(tau,R,tauw,Rw) % Leistungsanpassung bei tau=0, Modulation bleibt
Signalausschnitt (ohne Modulation der Amplitude) eines periodischen Signals, passt nicht mit vollen Perioden in das Beobachtungsfenster, …
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(2*pi*t*28/float(N))+normal(0,0.1,N)
M=300
w=ones(M)
tw=arange(0,M)*dt
xw=x[0:M]*w
plot(t,x,tw,xw)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(2*pi*t*28/N)+normrnd(0,0.1,[1,N]);
M=300;
w=ones(1,M);
tw=(0:M-1)*dt;
xw=x(1:M).*w;
plot(t,x,tw,xw)
… Leistungsspektrum des periodischen Signals (verschmiert) …
f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
X=fft(x)
P=abs(X)**2/N**2
fw=roll(arange(-(M//2),M-(M//2))/float(M*dt),M-(M//2))
Xw=fft(xw)
Pw=abs(Xw)**2/M**2
semilogy(roll(f,N//2),roll(P,N//2),roll(fw,M//2),roll(Pw,M//2))
show()
f=circshift((-fix(N/2):N-fix(N/2)-1)/(N*dt),N-fix(N/2));
X=fft(x);
P=abs(X).^2/N^2;
fw=circshift((-fix(M/2):M-fix(M/2)-1)/(M*dt),M-fix(M/2));
Xw=fft(xw);
Pw=abs(Xw).^2/M^2;
semilogy(circshift(f,fix(N/2)),circshift(P,fix(N/2)),circshift(fw,fix(M/2)),circshift(Pw,fix(M/2)))
… und Korrelationsfunktion (zeigt Interferenz)
tau=arange(0,N)*dt
R=N*real(ifft(P))
tauw=arange(0,M)*dt
Rw=M*real(ifft(Pw))
plot(tau,R,tauw,Rw)
show()
tau=(0:N-1)*dt;
R=N*real(ifft(P));
tauw=(0:M-1)*dt;
Rw=M*real(ifft(Pw));
plot(tau,R,tauw,Rw)
Signalausschnitt (mit Modulation der Amplitude, um Sprünge des periodisch fortgesetzten Signalausschnitts zu unterdrücken, gleich mit Leistungsanpassung) …
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(2*pi*t*28/float(N))+normal(0,0.1,N)
M=300
w=0.42-0.5*cos(pi*(2*arange(0,M)+1)/float(M))+0.08*cos(2*pi*(2*arange(0,M)+1)/float(M)) # Blackman
w*=sqrt(M/sum(w**2)) # Leistungsanpassung des Fensters
tw=arange(0,M)*dt
xw=x[0:M]*w
plot(t,x,tw,xw)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(2*pi*t*28/N)+normrnd(0,0.1,[1,N]);
M=300;
w=0.42-0.5*cos(pi*(2*(0:M-1)+1)/M)+0.08*cos(2*pi*(2*(0:M-1)+1)/M); % Blackman
w=w*sqrt(M/sum(w.^2)); % Leistungsanpassung des Fensters
tw=(0:M-1)*dt;
xw=x(1:M).*w;
plot(t,x,tw,xw)
… Leistungsspektrum des periodischen Signals (bleibt verschmiert) …
f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
X=fft(x)
P=abs(X)**2/N**2
fw=roll(arange(-(M//2),M-(M//2))/float(M*dt),M-(M//2))
Xw=fft(xw)
Pw=abs(Xw)**2/M**2
semilogy(roll(f,N//2),roll(P,N//2),roll(fw,M//2),roll(Pw,M//2))
show()
f=circshift((-fix(N/2):N-fix(N/2)-1)/(N*dt),N-fix(N/2));
X=fft(x);
P=abs(X).^2/N^2;
fw=circshift((-fix(M/2):M-fix(M/2)-1)/(M*dt),M-fix(M/2));
Xw=fft(xw);
Pw=abs(Xw).^2/M^2;
semilogy(circshift(f,fix(N/2)),circshift(P,fix(N/2)),circshift(fw,fix(M/2)),circshift(Pw,fix(M/2)))
… und Korrelationsfunktion (Leistungsanpassung ist gegeben, Interferenz bleibt)
tau=arange(0,N)*dt
R=N*real(ifft(P))
tauw=arange(0,M)*dt
Rw=M*real(ifft(Pw))
plot(tau,R,tauw,Rw)
show()
tau=(0:N-1)*dt;
R=N*real(ifft(P));
tauw=(0:M-1)*dt;
Rw=M*real(ifft(Pw));
plot(tau,R,tauw,Rw)
Lösung nur durch geeignete Normierung (ohne Fensterfunktion), …
N=1000
dt=1.0
t=arange(0,N)*dt
x=sin(2*pi*t*28/float(N))+normal(0,0.1,N)
M=300
w=ones(M)
tw=arange(0,M)*dt
xw=x[0:M]*w
plot(t,x,tw,xw)
show()
N=1000;
dt=1;
t=(0:N-1)*dt;
x=sin(2*pi*t*28/N)+normrnd(0,0.1,[1,N]);
M=300;
w=ones(1,M);
tw=(0:M-1)*dt;
xw=x(1:M).*w;
plot(t,x,tw,xw)
… Leistungsspektrum des Signals unter der Annahme einer Periodenlänge von 250 Werten …
NT=250 #Periodenlägne
f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
X=fft(x)
P=abs(X)**2/N**2
xp=zeros(NT)
np=zeros(NT)
for i in range(0,M):
  xp[i%NT]+=x[i]
  np[i%NT]+=1

X1=fft(xp)
X0=fft(np)
P1=abs(X1)**2/NT**2
P0=abs(X0)**2/NT**2
R=real(ifft(P1))/real(ifft(P0))
fw=roll(arange(-(NT//2),NT-(NT//2))/float(NT*dt),NT-(NT//2))
Pw=real(fft(R))/float(NT)
semilogy(roll(f,N//2),roll(P,N//2),roll(fw,NT//2),roll(Pw,NT//2))
show()
NT=250 %Periodenlänge
f=circshift((-fix(N/2):N-fix(N/2)-1)/(N*dt),N-fix(N/2));
X=fft(x);
P=abs(X).^2/N^2;
xp=zeros(1,NT);
np=zeros(1,NT);
for i=1:M
xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i);
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
X1=fft(xp);
X0=fft(np);
P1=abs(X1).^2/NT^2;
P0=abs(X0).^2/NT^2;
R=real(ifft(P1))./real(ifft(P0));
fw=circshift((-fix(NT/2):NT-fix(NT/2)-1)/(NT*dt),NT-fix(NT/2));
Pw=real(fft(R))/NT; semilogy(circshift(f,fix(N/2)),circshift(P,fix(N/2)),circshift(fw,fix(M/2)),circshift(Pw,fix(M/2)))
… und Korrelationsfunktion
tau=arange(0,N)*dt
R=N*real(ifft(P))
tauw=arange(0,NT)*dt
Rw=NT*real(ifft(Pw))
plot(tau,R,tauw,Rw)
show()
tau=(0:N-1)*dt;
R=N*real(ifft(P));
tauw=(0:NT-1)*dt;
Rw=NT*real(ifft(Pw));
plot(tau,R,tauw,Rw)