Analiza sygnałów z wykorzystaniem DFT
Analiza DFT sygnałów ciągłych (analogowych)
Przy analizie sygnałów ciągłych za pomocą DFT niezbędne są następujące etapy przetwarzania:
1) filtracja antyaliasingowa (przed przetwornikiem AC)
2) okienkowanie sygnału dyskretnego x[n] (tzn. mnożenie przez współczynniki okna)
Ą
1
j� jŚ j(�-Ś)
V (e ) = )dŚ
v[n] = w[n]x[n]
+"X (e )W (e
,
2Ą
-Ą
zawsze obserwujemy splot widma sygnału z widmem okna
(powoduje to wygładzenie widma sygnału)
"
j� - j� n
X (e ) =
3) próbkowanie ciągłego widma sygnału dyskretnego "x[n]e przez DFT
n=-"
N -1
- j(2Ą / N )kn
j�
X[k] =
0 d" n d" N -1, 0 d" k d" N -1 X (e )
"x[n]e , , tj. wyznaczanie wartości dla pulsacji dyskretnych
n=0
� = (2Ą / N )k
).
Etapy 1-3 wpływają na obserwowane widmo DFT (co oznacza, że obserwujemy widmo sygnału
przez pryzmat zastosowanych narzędzi obliczeniowych).
1
 2
Analiza DFT sygnałów sinusoidalnych
x[n] = cos(�0n + Ć)
Transformata Fouriera ciągu harmonicznego składa się z impulsów występujących
"
j� jĆ
X (e ) = [Ąe � (� - �0 + 2Ąk) + Ąe- jĆ� (� + �0 + 2Ąk)].
"
okresowo w pulsacjach �
0
k=-"
1. Wpływ okienkowania:
" Okienkowanie sygnału sinusoidalnego w dziedzinie czasu powoduje rozmycie (zwiększenie
szerokości) reprezentacji Fouriera (widma).
" Długość okna N wpływa na gęstość próbkowania widma sygnału (współczynniki Fouriera
2Ą k / N
wyznaczane są z krokiem ).
Rozważmy sygnał ciągły w postaci sumy dwóch harmonicznych:
sc (t) = A0 cos(&!0t + Ś0 ) + A1 cos(&!1t + Ś1), - " < t < "
W wyniku idealnego próbkowania otrzymujemy sygnał dyskretny:
x[n] = A0 cos(�0n + Ś0 ) + A1 cos(�1n + Ś1), �0 = &!0T , �1 = &!1T , - " < n < "
natomiast po okienkowaniu
v[n] = w[n]x[n] = A0w[n]cos(�0n + Ś0) + A1w[n]cos(�1n + Ś1)
ciąg v[n] można przedstawić w postaci zespolonej:
A0 A0 A1 A1
jŚ0 j�0n jŚ1 j�1n
v[n] = w[n]e e + w[n]e- jŚ0e- j�0n + w[n]e e + w[n]e- jŚ1e- j�1n
2 2 2 2
3
a następnie wyznaczyć jego widmo korzystając z tw. o przesunięciu:
A0 jŚ0 j(�-�0) A0 jŚ0 j(�+�0) A1 jŚ1 j(�-�1) A1 jŚ1 j(�+�1)
j�
V (e ) = e W (e ) + e- W (e ) + e W (e ) + e- W (e )
2 2 2 2
Transformata Fouriera dyskretnego, okienkowanego sygnału sumy harmonicznych składa się z kopii
ą�0, ą �1
widm Fouriera okna powtarzanych w częstotliwościach i przeskalowanych przez
zespolone amplitudy poszczególnych harmonicznych.
Widmo okna Widmo sygnału okienkowanego złożonego z sumy dwóch
sinusoid
4
Listki boczne wpływają na listki główne powodując błędną interpretację (szczególnie
dobrze widoczne dla V2) jest to tzw. przeciek widmowy (leakage).
5
Przeciek widmowy może być przyczyną braku rozróżnienia częstotliwości o dużej
różnicy amplitud.
Częstotliwości położone zbyt blisko siebie mogą być nierozróżnialne.
6
Zwiększenie długość okna w czasie (w tym przypadku 5 razy) zmniejsza szerokość
listka głównego (tyle samo razy) w częstotliwości, co poprawia rozdzielczość
częstotliwościową.
Dwa główne skutki stosowania okna czasowego to:
" zmniejszona rozdzielczość częstotliwościowa i
" przeciek widmowy.
Rozdzielczość częstotliwościowa jest ograniczona przez
szerokość listka głównego widma okna, natomiast przeciek
widmowy przez tłumienie listków bocznych widma okna (w
stosunku do listka głównego).
Szerokość listka głównego widma okna zależy od długości okna
w czasie (im dłuższe okno tym węższy listek główny), natomiast
położenie listków bocznych zależy od rodzaju okna (i nie zależy
od jego długości).
Okno prostokątne jest oknem o najwęższym listku głównym i
najmniejszym tłumieniu listków bocznych (16 dB).
W celu redukcji przecieku widmowego stosuje się okna inne niż
prostokątne (podobnie jak przy projektowaniu filtrów FIR), w
szczególności parametryczne okno Kaisera.
7
Przykład zmniejszenia przecieku widmowego (za cenę zmniejszenia rozdzielczości
częstotliwościowej) przez zastosowanie okna czasowego.
8
Wpływ próbkowania ciągłego widma sygnału dyskretnego
9
Gęstość próbkowania osi częstotliwości przez DFT można zwiększyć przez dodanie zer na końcu sygnału.
10
Krótkoczasowa Transformata Fouriera - spektrogram
Przy analizie DFT zakłada się, że częstotliwości występujące w sygnale są stałe w trakcie
obserwacji (nie zmieniają się przez całą długość okna czasowego), jednak większość sygnałów
posiada widmo częstotliwościowe zmienne w czasie (np. sygnał mowy lub sygnały
telekomunikacyjne), tzn. mogą zmieniać się amplitudy, fazy i częstotliwości składowych
harmonicznych tworzących sygnał. Krótkoczasową Transformatę Fouriera definiuje się następująco:
"
jm
X[n,) =
"x[n + m]w[m]e-
(*)
m=-"
gdzie w[n] oznacza współczynniki okna. Jednowymiarowy ciąg x[n] zostaje przedstawiony jako
X[n,)

dwuwymiarowa funkcja czasu (z n - dyskretnym) i częstotliwości (z ciągłym). są to
transformaty Fouriera fragmentów sygnału x[n] wybieranych przez okno w[m] (sygnał przesuwa się
przez okno).
m' = n + m
Zmieniając kolejność sumowania przez podstawienie otrzymujemy:
"
' j(n-m')
X[n,) =
"x[m ]w[-(n - m')]e-
m'=-"
jn
X[n,) = x[n]" h[n] h[n] = w[-n]e
Powyższe równanie jest splotem , gdzie , a w częstotliwości
j( -�)
H (ei� ) = W(e )
.
11
W definicji (*) położenie czasowe okna jest stałe (zafiksowane), a sygnał jest przesuwany przez
okno. Analogiczny rezultat można otrzymać przez przesuwanie okna po sygnale:
"
(
jm
X[n,) =
"x[m]w[m - n]e-
(**)
m=-"
(
X[n,) = e- jn X[n,)
Zależność pomiędzy definicjami (*) i (**) jest następująca .
Wpływ okienkowania
Stosowanie okna w krótkoczasowej transformacie Fouriera ma na celu wybór fragmentu
sygnału, w którym jego parametry są w przybliżeniu stacjonarne. Analiza sygnałów
szybkozmiennych wymaga stosowania krótkich okien czasowych, co pogarsza rozdzielczość
częstotliwościową. Wybór długości okna jest kompromisem pomiędzy rozdzielczością czasową i
częstotliwościową.
12
Próbkowanie w czasie i częstotliwości
X[n,) w[m] = 0, 0 d" m d" L - 1
N e" L
W praktyce oblicza się wartości za pomocą, DFT dla okna i :
L-1
j(2Ą / N )km
X[n,k] = X[n, = 2Ąk / N) = , 0 d" k d" N -1
"x[n + m]w[m]e-
m=0
zachodzi również:
j(2Ą / N )kn
X[n,k] = x[n]" hk[n], 0 d" k d" N -1 h[n] = w[-n]e
, gdzie
Powyższe równanie można interpretować jako bank N
filtrów, których odpowiedzi częstotliwościowe są
następujące:
j((2Ą k / N )-�)
Hk (ei� ) = W(e )
13
X[n, k]
Wprowadz
my do definicji dodatkowe oznaczenia (r - liczba całkowita, a R naturalna):
L-1
j(2Ą / N )km
X[rR, k] =
"x[rR + m]w[m]e-
m=0
X [k] = X[rR, k]
Upraszczając notację widzimy, że próbkowana, krótkoczasowa transformata
r
Fouriera składa się z N punktowych ciągów DFT liczonych z fragmentów sygnału x wybieranych
oknem w o długości L, które przesuwa się po sygnale, co R próbek. Aby istniała możliwość
N e" L e" R
rekonstrukcji sygnału musi zachodzić .
X[n, k]
Innym punktem spojrzenia na płaszczyznę jest stwierdzenie, że dla danego k otrzymujemy
ciąg będący efektem filtracji filtrem pasmowym. Co oznacza, że częstotliwość próbkowania każdego
2Ą / "ml "ml
z podpasm może być zredukowana , gdzie jest szerokością listka głównego transformaty
Fouriera okna.
14
Na powyższych rysunkach L=N=16. Pasma przepustowe filtrów zachodzą na siebie, również listki
boczne zachodzą na pasma przepustowe, co może być przyczyną aliasingu.
Najbardziej powszechnym zastosowaniem spektrogramu jest analiza dyskretnych sygnałów
niestacjonarnych, czyli takich, których parametry (tj. amplituda, częstotliwość lub faza) zmieniają
się w czasie (np. mowy, sygnały radarowe i inne).
Ograniczenia rozdzielczości częstotliwościowej i amplitudowej DFT dotyczą również spektrogramu
15
Sygnał o liniowej zmianie częstotliwości
16
Sygnał mowy, (kobieta) 'księżyc lśnił mocno noc całą'
Dla sygnałów rzeczywistych wystarczy obserwować
zakres od 0 do fp/2, dodatkowo dla sygnałów o dużej
dynamice (np. mowy) stosuje się skalę decybelową.
Parametrami wyznaczania spektrogramu są długość okna,
rodzaj okna, przesunięcie okna oraz długość FFT (można
stosować większą niż długość okna w celu poprawy
rozdzielczości częstotliwościowej).
L =256; %długość okna Hanninga,
Nfft=512; %długość fft,
M =L/2; %przesunięcie okna
17
Przykład: Detekcja sygnałów radarowych.
System radarowy składa się z: 1) anteny nadawczej i odbiorczej (często tej samej), 2) nadajnika mikrofali (np. impulsów
sinusoidalnych), 3) odbiornika wzmacniającego i detektującego echa impulsów odbitych od oświetlanych obiektów.
cos(&!0t) �(t)
Jeśli wysyłany sygnał radaru jest impulsem sinusoidalnym w postaci , a odległość od anteny wynosi to
s(t) = cos[&!0(t - 2�(t)/ c)]
odebrany (odbity od obiektu) sygnał ma postać (c - prędkość światła): .
�(t) = �0
Jeżeli obiekt się nie porusza względem anteny , to na podstawie pomiaru opóz
nienia można określić odległość
1
& &&
�(t) = �0 + �0t + �0t2 + ...
�(t)
obiektu od anteny. Jeśli jest funkcją czasu to po rozłożeniu na szereg Taylora otrzymujemy , gdzie
2!
& && &&
�0 �0 �0 �0 = 0
jest odległością nominalną, jest prędkością, jest przyspieszeniem itd. Przyjmując stałą prędkość tj.
&
s(t) = cos[(&!0 - 2&!0�0 /c)t - 2&!0�0 /c]
otrzymujemy: . Częstotliwość odebrana różni się od nadanej o częstotliwość Dopplera
&
&!d = -2&!0�0 /c
.
Na podstawie sygnału s(t) można estymować odległość obiektu od anteny (przez opóz
nienie) i prędkość
obiektu (przez częstotliwość Dopplera).
18
Analiza stacjonarnych sygnałów losowych - periodogram
W przypadku, gdy sygnał może być modelowany jako suma sinusoid DFT jest naturalnym
narzędziem analizy, natomiast dla sygnałów mających charakter szumowy stosuje się modele
sygnałów losowych z parametrami takimi jak: wartość średnia (DC), wariancja (moc średnia),
funkcja autokorelacji, gęstość widmowa mocy. Typowa estymata wartości średniej dla
skończonej realizacji stacjonarnego procesu losowego o długości L jest zdefiniowana jako średnia
próbek:
L-1
1
Ć
mx =
"x[n]
L
n=0
podobnie wariancja jest zdefiniowana jako wariancja próbek:
L-1
2
Ćx 1 Ć
� =
"(x[n] - mx )2
L
n=0
Oba powyższe estymatory są nieobciążone (unbiased) i asymptotycznie nieobciążone, tzn. wartość
2
Ć Ćx
mx mx �
oczekiwana osiąga prawdziwą średnią , a wartość oczekiwana wariancji osiągną
2
�
"
prawdziwą wariancję , gdy L osiąga . Dodatkowo oba estymatory są zgodne, tzn. poprawiają
x
"
estymację wraz ze wzrostem L, tzn. ich wariancja maleje do 0, gdy L osiąga .
19
Periodogram - estymacja widmowej gęstości mocy
Rozważmy estymatę widmowej gęstości mocy w postaci:
L-1 L-1
1
1 j� - j�n 2
j�
I(�) = | V (e ) |2
"w[n]x[n]e "(w[n])
, gdzie: V (e ) = , a U = jest stałą normalizującą.
L
LU
n=0 n=0
Dla okna prostokątnego taki estymator nazywa się periodogramem, a dla okna innego niż
I(�)
prostokątne nazywa się zmodyfikowanym periodogramem. Wzór definicyjny można
przedstawić w postaci:
L-1 L-1
1
j�m
I(�) = cvv[m] =
vv
"c [m]e- "x[n]w[n]x[n + m]w[n + m] v[n] = x[n]w[n]
, gdzie jest korelacją ciągu
LU
m=-(L-1) n=0
I(�)
Periodogram jest transformatą Fouriera korelacji okienkowanego ciągu. Wartości są
rzeczywiste i parzyste względem połowy częstotliwości próbkowania.
Periodogram wyznacza się dla skończonego zbioru częstotliwości dyskretnych za pomocą DFT.
Ponieważ sygnał losowy może mieć dużą składową stałą, która ze względu na przecieki i
rozdzielczość amplitudową widma utrudniałaby analizę, wiec przed liczeniem periodogramu
odejmuje się od sygnału wartość średnią.
Periodogram ze stałą normalizującą U jest asymptotycznie nieobciążony, tzn. obciążenie osiąga
zero wraz ze wzrostem długości okna L.
Periodogram nie jest estymatorem zgodnym, ponieważ jego wariancja nie dąży do zera wraz ze
wzrostem długości okna L.
20
Uśrednianie zmodyfikowanego periodogramu -metoda Welcha
Ciąg danych x[n] jest dzielony na fragmenty o długości L -próbek, które następnie są okienkowane:
xr[n] = x[rR + n]w[n], 0 d" n d" L -1
Jeżeli Rposzczególnych segmentów o połowę, gdyż powoduje to najlepsze zmniejszenie wariancji
estymatora. Periodogram r -tego segmentu jest obliczany jako:
1
j�
Ir (�) = | Xr (e ) |2
LU
Następnie K periodogramów poszczególnych segmentów jest uśredniane:
K-1
1
IĆ(�) =
r
"I (�)
K
r=0
Uśrednianie periodogramów daje asymptotycznie nieobciążony i zgodny estymator widmowej
gęstości mocy. W praktyce do liczenia uśrednionych periodogramów stosuje sie DFT, wówczas:
K-1
1
1
IĆ[k] =
Ir[k] = | X [k] |2
r
"I (k)
r
,
K
LU
r=0
21
x=sin(2*pi*100*t)+2.5*randn(1,N); x=x-mean(x);
22
 23