Widmo sygnału
Przypadek ciągły – z rozwinięcia w szereg Fouriera
X =
∫
− ∞
∞
x t e
−
j t
dt gdzie = 2 f a w efekcie X f =
∫
−∞
∞
x t e
−
j 2 f t
dt
transformata odwrotna
x t =
1
2
∫
−∞
∞
X f e
j 2 f t
df
Funkcje e
−
j 2 f t
są do siebie ortogonalne więc stanowią bazę przestrzeni !!!
Przypadek dyskretny – równanie DFT:
X k =
∑
n =0
N −1
x ne
−
j 2 k n / N
gdzie 0≤k ≤ N −1 (dyskretne częstotliwości)
formuła wyznaczania częstotliwości dyskretnych: f
k
=
k F
s
N
(dlaczego akurat takie? !!!)
tr. odwrotna:
x n=
1
N
∑
k=0
N−1
X k e
j 2 k n / N
gdzie
0≤n≤N −1
Jeżeli wynik DFT zapiszemy jako liczbę zespoloną
X k =A k e
j k
=
X
r
k j X
i
k
to:
–
moduł widma
∣
X k
∣
=
Ak =
X
r
2
k X
i
2
k
–
faza widma
arg X k =k =arctan
X
i
k
X
r
k
–
widmowa gęstość mocy
P k =∣X k ∣
2
=
X k conj X k =[ X
r
k j X
i
k ] [ X
r
k − j X
i
k ]= X
r
2
k X
i
2
k
Transformata DFT dla sygnałów 2D
Jeżeli mamy obraz np.
o x , y
x , y ∈ℤ
O k , l=
∑
x=0
N
x
−
1
e
−
j 2 k x / N
x
∑
y=0
N
y
−
1
o x , y e
−
j 2l y/ N
y
=
∑
x=0
N
x
−
1
∑
y=0
N
y
−
1
o x , y e
−
j 2 l y / N
y
e
−
j 2 k x/ N
x
Własności
Symetria
Dla sygnałów rzeczywistych ciągłych !!! i dyskretnych zachodzi:
∣
X k ∣=∣X −k ∣ i k =−−k
lub inaczej X k =conj X −k ,
X
r
k j X
i
k = X
r
−
k − j X
i
−
k (można pokazać
przez wstawienie -k do równania DFT)
Liniowość
Zachodzi dla sygnałów ciągłych i dyskretnych:
dla x n=a y nb z n mamy
X k =a Y k b Z k
Okresowość
Dla dyskretnych
X k =X m∗N k , m∈ℤ (pokazać)
Przesunięcie w czasie
dla
x n= y nn
0
mamy
X k =
∑
k =0
N −1
y ne
−
j
2 nn
0
N
k
=
e
−
j
2 n
0
N
k
Y k =e
j
0
k
Y k
(przesunięcie fazy widma o stały czynnik)
Przesunięcie w częstotliwości
dla X k =Y k k
0
=
1
N
∑
k=0
N −1
Y k k
0
e
−
j
2kk
0
N
n
mamy
x n=e
j
2 k
0
N
n
y n=e
j
0
n
y n
(przemnożenie przez stałą częstotliwość)
Istnienie FT/DFT dla sygnałów okresowych i nieokresowych
DFT istnieje tylko dla sygnałów okresowych !!!
(proszę sobie przypomnieć przykład z aproksymacją linii prostej – Notatki4 BS)
Jeżeli sygnał transformowany nie jest okresowy to konsekwencją jest przeciek widma.
N=512; Fs=512; n=(0:N-1)/Fs; x=sin(2*pi*50*n); y=sin(2*pi*50.17*n);
plot(n,x,n,y);
X = fft(x); Y = fft(y); f = ((0:N-1)/N)*Fs;
plot(f,abs(X),'b*',f,abs(Y),'r*');
Najpowszechniejsze lekarstwo – okienkowanie sygnału
x
w
n=w n x n
N=256; Fs=N; n=(0:N-1)./Fs;
x=sin(2*pi*20*n); y=sin(2*pi*20.17*n);
plot(n,x,';x(n);',n,y,';y(n);');
w = bartlett(N)';plot(n,w,n,y.*w); %% okienko trójkątne
w = hamming(N)';plot(n,w,n,y.*w);
w = hanning(N)';plot(n,w,n,y.*w);
w = gausswin(N,3)';plot(n,w,n,y.*w);
xw= x.*w; yw = y.*w;
X=fft(x);XW=fft(xw);Y=fft(y);YW = fft(yw);f = ((0:N-1)/N)*Fs;
plot(f,abs(Y),'ro',f,abs(YW),'b*'); %% zmniejszony przeciek, ale nic za darmo
plot(f,abs(X),'ro',f,abs(XW),'b*');
Dlaczego tak się dzieje (rozmycie głównego prążka)?
x
w
n=w n x n
X
w
k =X k ∗W k
wb = bartlett(N);
whm = hamming (N);
whn = hanning (N);
wg3 = gausswin (N,3); wg10 = gausswin(N,10);
plot(n,wb,n,whm,n,whn,n,wg3,n,wg10);
Wb = fft(wb);
Whn = fft(whn);
Whm = fft(whm);
Wg3=fft(wg3); Wg10=fft(wg10);
f = (0:N-1)/N;
plot(f,log10(abs(Wb)),f,log10(abs(Whn)),f,log10(abs(Whm)));
plot(f,log10(abs(Wg3)),f,log10(abs(Wg10)));
Zwiększanie rozdzielczości częstotliwościowej
Zwiększenie N powoduje więcej prążków częstotliwości – skąd wziąć dodatkowe próbki ?
N=32;n=(0:N-1);Fs=1;x=sin(2*pi*.27*n);plot(n,x);
X=fft(x);f=(-N/2:N/2-1)/N;plot(f,fftshift(abs(X)));
N=32;n=(0:N-1);Fs=1;x=sin(2*pi*.27*n);
y=[x,zeros(1,N)];M=size(y,2);m=(0:M-1);plot(m,y);
X=fft(x);fx=(-N/2:N/2-1)/N; Y=fft(y);fy=(-M/2:M/2-1)/M;
plot(fx,fftshift(abs(X)),fy,fftshift(abs(Y)));
FFT
Koszt obliczeniowy
DFT = N
2
FFT =
N /2 log
2
N
FFT o podstawie 2
X k =
∑
n=1
N
x ne
−
j 2 k n/ N
podstawmy
W
N
=
e
−
j2/ N
X k =
∑
n=1
N
x nW
N
k n
i podzielmy próbki na parzyste i nieparzyste:
X k =
∑
n=1
N / 2
x 2nW
N
2nk
∑
n=1
N / 2
x 2n1W
N
2n 1 k
=
∑
n=1
N /2
x 2n W
N
2nk
W
N
k
∑
n=1
N / 2
x 2n1W
N
2nk
ponieważ
W
N
2
=
e
−
j 2 2 / N
=
e
−
j 2 /N /2
=
W
N / 2
to
X k =
∑
n=1
N / 2
x 2nW
N /2
nk
W
N
k
∑
n=1
N / 2
x 2n 1 W
N / 2
nk
Podzielmy teraz częstotliwości na dwie połowy l dla 1≤k ≤ N /2 i
lN /2 dla N /21≤k ≤N
X lN /2=
∑
n=1
N /2
x 2n W
N / 2
n l N / 2
W
N
l N / 2
∑
n=1
N / 2
x 2n1W
N / 2
n l N / 2
zobaczmy, że
W
N / 2
n l N / 2
=
W
N / 2
nl
W
N / 2
nN /2
=
W
N /2
nl
e
−
j 2 n N /2 / N / 2
=
W
N / 2
nl
e
−
j 2 n
=
W
N / 2
nl
oraz
W
N
l N /2
=
W
N
l
W
N
N / 2
=
W
N
l
e
−
j 2 n N / 2/ N
=
W
N
l
e
−
j n
=
W
N
l
−
1=−W
N
l
tak więc:
X lN /2=
∑
n=1
N /2
x 2n W
N / 2
nl
−
W
N
l
∑
n=1
N /2
x 2n1W
N /2
nl
i dla dolnych częstotliwości (poniżej N/2)
X l=
∑
n =1
N /2
x 2n W
N / 2
nl
W
N
l
∑
n =1
N / 2
x 2n1W
N / 2
nl
czyli:
X l = AlW
N
l
B l
X lN /2 = Al−W
N
l
B l