DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
1
Próbkowanie i transformacje Fouriera
Zawartość instrukcji:
1 Materiał z zakresu DSP
1.1 Próbkowanie sygnału
e
j
v t
⋅ ⋅ ⋅ ⋅
2
π
1.1.2 Częstotliwość cyfrowa i pulsacja cyfrowa
1.2 Odpowiedź częstotliwościowa systemu
1.2.1 Odpowiedź częstotliwościowa jako fragment transformaty “Z” odpowiedzi
impulsowej
1.3 Transformacja Fouriera z “czasem” dyskretnym
1.3.1 Porównanie Transformacji Fouriera z czasem dyskretnym z Transformacją Fouriera
dla sygnałów z czasem ciągłym
1.4 Dyskretna Transformacja Fouriera
1.5 Transformata DFT jako próbkowanie transformaty D-TFT
2 Korzystanie z pakietu MATLAB
2.1 Opis wybranych funkcji
2.2 Przykłady realizacji wybranych zadań
2.2.1 Demonstracja okresowości sygnału kosinusoidalnego względem “czasu” i
częstotliwości
2.2.2 Wyznaczanie odpowiedzi częstotliwościowej filtru FIR rzędu N-1 dla pojedynczej
wartości częstotliwości cyfrowej
2.2.3 Wyznaczanie transmitancji “Z” filtru typu FIR o zadanej odpowiedzi impulsowej dla
pojedynczej wartości z
2.2.4 Wyznaczanie przybliŜenia transformaty D-TFT ze wzoru definicyjnego (27):
2.2.5 Wyznaczanie transformaty DFT ze wzoru definicyjnego (33):
2.2.6 Wyznaczanie przybliŜenia transformaty D-TFT za pomocą funkcji fft:
2.2.7 Wykres fazy - efekt ograniczania do zakresu od minus pi do pi
3 Dodatek teoretyczny nr 2 - Relacja pomiędzy Transformacją Fouriera z Czasem
Dyskretnym (ang. D-TFT) i Transformacją Fouriera z Czasem Ciągłym (ang.
C-TFT)
4 Dodatek teoretyczny nr 3 - Relacja pomiędzy Transformacją Fouriera z Czasem
Dyskretnym (ang. D-TFT) i Dyskretną Transformacją Fouriera (ang. DFT);
Relacja pomiędzy DFT i C-TFT
1 Materiał z zakresu DSP
1.1 Próbkowanie sygnału
e
j
v t
⋅ ⋅ ⋅ ⋅
2
π
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
2
W podrozdziale tym będzie mowa o próbkowaniu rozumianym jako tworzenie ciągu wartości
chwilowych funkcji pobranych w równych odstępach “czasu”:
s n
s n
t
[ ]
(
)
=
⋅ ∆
, gdzie przez
n
oznaczono dowolną liczbę całkowitą. Odstęp ten będzie dalej oznaczany jako
∆
t
T
=
.
Sygnał
e
j
v t
⋅ ⋅ ⋅ ⋅
2
π
jest zespoloną funkcją określoną na ciągłej zmiennej
t
. Zmienną
v
moŜna
traktować jako parametr, określający częstotliwość sygnału, analogicznie jak w przypadku funkcji
trygonometrycznych, poniewaŜ:
(
)
(
)
s t
e
v t
j
v t
C
j
v t
( )
cos
sin
=
=
⋅ ⋅ ⋅ + ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
2
2
2
π
π
π
(1)
Sygnał cyfrowy pochodzący z próbkowania powyŜszego sygnału moŜna zapisać w następujący sposób:
s n
e
j
v n T
[ ]
=
⋅ ⋅ ⋅ ⋅ ⋅
2
π
(2)
Interpretując (1) moŜna odwrócić role i zmienną
t
przyjąć jako parametr, natomiast (1) określić jako
funkcję zmiennej
v
. Łatwo zauwaŜyć, Ŝe właściwości funkcji będą identyczne, gdyŜ od
t
zaleŜy ona
w identyczny sposób jak od
v
. MoŜna zatem stwierdzić, Ŝe jest to funkcja dwóch zmiennych (
t
oraz
v
) i jest to funkcja symetryczna względem tych zmiennych. Zatem, jeŜeli załoŜy się ustaloną wartość
t
, co z kolei daje
n T
const
⋅ =
, to widać natychmiast, Ŝe funkcja (1) jest okresowa względem
zmiennej
v
.
Funkcja
( )
( )
e
j
j
⋅
=
+ ⋅
ϕ
ϕ
ϕ
cos
sin
(3)
jest okresowa, przy czym okres podstawowy wynosi
∆
ϕ
π
0
2
= ⋅
, z czego wynika, Ŝe kaŜda wartość
∆
ϕ
π
= ⋅ ⋅
n 2
jest takŜe okresem funkcji (3). Z porównania (3) oraz (1) wynika, Ŝe dla kaŜdej
ustalonej wartości
t
n T
n
= ⋅
wyraŜenie okresu sygnału (1) w odniesieniu do zmiennej
v
daje:
∆
v
T
=
1
(4)
W analogiczny sposób moŜna wyobrazić sobie, iŜ ciąg próbek (2) jest funkcją zmiennej całkowitej
n
oraz częstotliwości
v
. Ustalając
n
n
ust
=
otrzymuje się wartość próbki
s n
s v
ust
[
]
( )
=
jako funkcję
ciągłej zmiennej
v
.
PoniewaŜ ciąg
s n
[ ]
powstaje z wartości (1) pobranych właśnie w chwilach
t
n
, więc ciąg ten
wykazuje okresowość względem zmiennej
v
, z okresem (4).
Oznacza to, Ŝe przy ustalonym
T
dla kaŜdej wartości
n
wartość odpowiedniej próbki sygnału
rozpatrywana w funkcji
v
jest okresowa. Okres (4) jest taki sam dla wszystkich wartości
n
. Zatem po
spróbkowaniu (1) i otrzymaniu ciągu w postaci (2) nie da się juŜ w oparciu o ten ciąg rozróŜnić, czy
sygnał sprzed próbkowania (czyli (1)) miał częstotliwość
v
, czy teŜ
v
T
k
+ ⋅
1
(dla dowolnej liczby
całkowitej
k
). Efekt ten jest nazywany
aliasingiem. Jak łatwo zauwaŜyć, identyczny efekt powstaje
takŜe w przypadku próbkowania funkcji sinus lub teŜ kosinus (dla kaŜdej z osobna).
PoniewaŜ kaŜdy sygnał z czasem ciągłym posiadający transformatę Fouriera moŜna zinterpretować jako
kombinację liniową nieskończonej ilości sygnałów typu (1), więc efekt aliasingu dotyczy w taki sam
sposób kaŜdego takiego sygnału. RozwaŜanie efektu aliasingu powstającego w wyniku próbkowania
prowadzi do opisanego w dodatku teoretycznym nr 1 (patrz instrukcja do ćwiczenia 1) twierdzenia o
próbkowaniu.
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
3
1.1.2 Częstotliwość cyfrowa i pulsacja cyfrowa
W dalszych rozwaŜaniach przydatne będzie pojęcie częstotliwości dla sygnałów cyfrowych. PoniewaŜ
okres próbkowania
T
jest wartością stałą, więc opis sygnału w postaci (2) moŜna nieco uprościć
podstawiając:
f
v T
= ⋅
(5)
Pozwala to na przepisanie (2) do postaci:
s n
e
e
j
v n T
j
f n
[ ]
=
=
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
2
2
π
π
(6)
Bezwymiarowa wielkość
f
jest w tym przypadku wielkością ciągłą, jednak ze względu na powiązanie
z sygnałami cyfrowymi będzie dla uproszczenia nazywana skrótowo częstotliwością cyfrową.
ZaleŜność (5) jest dość istotna, poniewaŜ przedstawia, jaka jest zaleŜność pomiędzy częstotliwością dla
sygnałów z czasem ciągłym
v
, wyraŜoną w Hz, oraz bezwymiarową częstotliwością cyfrową.
Niekiedy spotyka się opis sygnału ciągłego (1) zawierający zamiast częstotliwości
v
wyraŜanej w Hz
pulsację
Ω
w radianach na sekundę:
s
t
e
C
j
t
( )
=
⋅Ω⋅
(7)
Postępując podobnie jak w przypadku częstotliwości, moŜna zatem wprowadzić pojęcie pulsacji
cyfrowej, wyraŜanej w radianach:
ω
π
=
⋅
=
⋅ ⋅
Ω
T
f
2
(8)
Sygnał cyfrowy (2) moŜna zatem zapisać równieŜ tak:
s n
e
j
n
[ ]
=
⋅ ⋅
ω
(9)
Warto zauwaŜyć, Ŝe w literaturze anglo-języcznej spotyka się słowo frequency w odniesieniu do:
a) częstotliwości w Hz, b) pulsacji w radianach na sekundę, c) bezwymiarowej częstotliwości cyfrowej
oraz d) pulsacji cyfrowej w radianach. NaleŜy pamiętać, Ŝe wielkości te są wprawdzie powiązane
wzajemnie (zaleŜności (5) i (8)), jednak mimo wszystko są to cztery róŜne wielkości, wyraŜane za
pomocą róŜnych jednostek. Dlatego warto się upewnić jakiego znaczenia słowa frequency uŜywa w
danym przypadku autor.
Efekt aliasingu powoduje, Ŝe sygnał cyfrowy (2) jest okresowy z punktu widzenia parametru
v
i okres
podstawowy wynosi
1
T
. Z zaleŜności (5) oraz (8) wynika, Ŝe sygnał cyfrowy jest okresowy równieŜ ze
względu na
f
oraz
ω
, przy czym w kaŜdym przypadku okres podstawowy jest równy odpowiednio
1
oraz
2
⋅
π
. Okresowość ciągu
s n
[ ]
wyraŜoną względem
f
lub
ω
moŜna zapisać następująco (
k
to dowolna liczba całkowita):
(
)
s n
e
e
f
j
f n
j
f
k n
[ ]
=
=
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ + ⋅
2
2
π
π
(10)
lub teŜ
(
)
s n
e
e
j
n
j
k
n
ω
ω
ω
π
[ ]
=
=
⋅ ⋅
⋅ + ⋅ ⋅ ⋅
2
(11)
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
4
We wzorach (10) i (11) zaznaczono dodatkowo, Ŝe sygnały są uzaleŜnione od wartości parametru, który
moŜe przyjmować dowolne wartości rzeczywiste - odpowiednio:
f
dla (10) oraz
ω
dla (11).
Warto tu wspomnieć o pewnej ciekawostce, iŜ nie zawsze ciąg pochodzący z próbkowania sygnału
okresowego musi być teŜ okresowy. W celu wykazania, Ŝe jest to moŜliwe sformułujemy warunek, jaki
musi być spełniony, by w takim przypadku ciąg próbek był okresowy. JeŜeli warunek ten nie będzie
spełniony, to ciąg próbek nie będzie okresowy.
Rozpatrujemy sygnał jako funkcję czasu ciągłego. W konsekwencji ciąg próbek jest funkcją zmiennej
n
. Niech sygnał z czasem ciągłym
x t
( )
będzie okresowy z okresem
P
:
x t
x t
P
( )
(
)
=
+
(12)
Przyjmując, Ŝe ciąg próbek tego sygnału
x n
x n T
[ ]
(
)
=
⋅
jest takŜe okresowy, otrzymujemy:
x n
x n
N
[ ]
[
]
=
+
(13)
przy czym oczywiście wartość okresu
N
musi być liczbą naturalną. Z porównania (12) i (13) oraz
zaleŜności
t
n T
n
= ⋅
wynika, Ŝe aby zachodziły obie równości ((12) i (13)) muszą istnieć dwie liczby
naturalne
k
oraz
m
, takie Ŝe:
k P
m N T
⋅ = ⋅ ⋅
(14)
Równość (14) moŜna opisać tak, Ŝe musi istnieć taka całkowita wielokrotność okresu
P
, by takiemu
odcinkowi czasu odpowiadała całkowita ilość okresów ciągu (czyli
N T
⋅
).
Z (14) wynika wprost poszukiwany warunek:
P
T
m N
k
=
⋅
(15)
Prawa strona warunku (15) oznacza liczbę wymierną, poniewaŜ wszystkie występujące w niej liczby są
naturalne. Zatem, aby w wyniku próbkowania sygnału okresowego powstał ciąg nieokresowy,
wystarczy by iloraz okresu
P
sygnału do okresu próbkowania
T
był liczbą niewymierną. Nie jest to
takie trudne do spełnienia, przykładowo: dla sygnału sinusoidalnego o częstotliwości
v
=
100
π
Hz
próbkowanego z okresem
T
=
0 01
.
s otrzymuje się, Ŝe:
P
T
T v
=
⋅
=
1
π
(16)
1.2 Odpowiedź częstotliwościowa systemu
Sygnał cyfrowy (2) nie moŜe być naturalnie rezultatem próbkowania sygnału pochodzącego ze świata
rzeczywistego. Jest to sygnał o charakterze teoretycznym, jednak jest on dość waŜny, poniewaŜ
umoŜliwia wygodną interpretację pewnych powszechnie stosowanych pojęć. Przyjmijmy, Ŝe na wejście
systemu cyfrowego, liniowego i stacjonarnego, czyli opisanego przez odpowiedź impulsową
h n
[ ]
(odpowiedź ta moŜe być skończona lub nieskończona), zostanie podany sygnał w postaci ciągu (2).
Okazuje się, Ŝe odpowiedź systemu na ten sygnał będzie takim samym sygnałem, jedynie pomnoŜonym
przez pewną stałą liczbę zespoloną, zaleŜną jedynie od wartości
f
, a niezaleŜnej od indeksu
n
. JeŜeli
odpowiedź systemu
swy
n
f
[ ]
na sygnał typu (10) zapiszemy w postaci splotu tego sygnału
wejściowego i odpowiedzi impulsowej systemu
h n
[ ]
, to:
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
5
swy n
e
H f
h n
e
H f
swy n
s n
f
j
f n
j
f n
f
f
[ ]
( )
[ ]
( )
[ ]
[ ]
=
⋅
=
∗
⇓
=
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
2
2
π
π
(17)
gdzie
H f
( )
jest dla ustalonej wartości
f
wielkością stałą, niezaleŜną od
n
, natomiast ogólnie jest
to zespolona funkcja zmiennej
f
.
Słuszność (17) moŜna wykazać podstawiając do tej zaleŜności pełny wzór na splot:
(
)
h n
e
h k e
e
h k e
j
f n
j
f n k
k
j
f n
j
f k
k
[ ]
[ ]
[ ]
∗
=
⋅
=
⋅
⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ −
=−∞
∞
⋅ ⋅ ⋅ ⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
∑
∑
2
2
2
2
π
π
π
π
(18)
Z porównania (17) i (18) wynika, Ŝe:
H f
h k
e
j
f k
k
( )
[ ]
=
⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
∑
2
π
(19)
Określona przez (19) funkcja zmiennej
f
jest nazywana odpowiedzią częstotliwościową systemu. Jak
wynika z przeprowadzonych rozwaŜań,
H f
( )
jest zespoloną funkcją częstotliwości cyfrowej
i przedstawia zaleŜność pomiędzy sygnałem wyjściowym i wejściowym systemu liniowego
stacjonarnego, gdy sygnał wejściowy jest w postaci (2) (lub teŜ (6)). Odpowiedź częstotliwościowa
moŜe być wyraŜona takŜe jako funkcja pulsacji cyfrowej:
H
h k
e
j
k
k
( )
[ ]
ω
ω
=
⋅
− ⋅ ⋅
=−∞
∞
∑
(20)
Ze wzorów (19) oraz (20) widać, iŜ wszystkie składniki odpowiednich sum posiadają pewien wspólny
okres. Dla (19) jest to
∆
f
=
1
, dla (20)
∆
ω
π
= ⋅
2
. Zatem odpowiedź częstotliwościowa jako
kombinacja liniowa takich składników posiada taki sam okres. MoŜna to wyjaśnić równieŜ w oparciu
o zaleŜność (17) - poniewaŜ sygnał wejściowy jest okresowy ze względu na częstotliwość
f
, zatem
odpowiedź częstotliwościowa jest takŜe okresowa.
PoniewaŜ okresowość odpowiedzi częstotliwościowej wynika z okresowości (względem parametru
f
lub
ω
) sygnału wejściowego, więc dla zaznaczenia tego faktu stosuje się często takŜe następujące
oznaczenie:
H e
h k
e
j
f
j
f k
k
(
)
[ ]
⋅ ⋅ ⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
=
⋅
∑
2
2
π
π
(21)
lub
H e
h k
e
j
j
k
k
(
)
[ ]
⋅
− ⋅ ⋅
=−∞
∞
=
⋅
∑
ω
ω
(22)
1.2.1 Odpowiedź częstotliwościowa jako fragment transformaty “Z” odpowiedzi impulsowej
JeŜeli podstawimy:
e
z
j
f
⋅ ⋅ ⋅
=
2
π
(23)
to wzór (21) przyjmie postać transformaty “Z” odpowiedzi impulsowej
h n
[ ]
:
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
6
H z
h n z
n
n
( )
[ ]
=
⋅
−
=−∞
∞
∑
(24)
Oznacza to, Ŝe odpowiedź częstotliwościowa to fragment transformaty “Z” odpowiedzi impulsowej,
określony dla wartości “z” spełniających zaleŜność (23). ZaleŜność ta wyznacza na płaszczyźnie
zespolonej “z” okrąg jednostkowy o środku w punkcie (0,0). Warto pamiętać, Ŝe gdyby okrąg ten nie
zawierał się w obszarze zbieŜności transformaty “Z”, to dla takiego systemu nie dałoby się wyznaczyć
odpowiedzi częstotliwościowej.
Rys.1. Płaszczyzna zespolona zmiennej “z” z zaznaczonym przykładowym obszarem zbieŜności
transformaty oraz okręgiem określającym fragment, którego dotyczy odpowiedź częstotliwościowa
PoniewaŜ odpowiedź częstotliwościowa jest dla okręgu jednostkowego na płaszczyźnie “z”
równoznaczna z transformatą “Z” odpowiedzi impulsowej, więc dla tego samego okręgu jest ona takŜe
równoznaczna z transmitancją danego systemu. MoŜna zatem, dla tak ograniczonego zbioru wartości
“z”, odpowiedni wzór wielomianowy na transmitancję przepisać do następującej postaci:
H e
b
b
e
b
e
b
e
a
e
a
e
a
e
j
f
j
f
j
f
M
j
f
M
j
f
j
f
K
j
f
K
(
)
(
)
(
)
⋅ ⋅ ⋅
− ⋅ ⋅ ⋅
− ⋅ ⋅ ⋅
− ⋅ ⋅ ⋅ ⋅
−
− ⋅ ⋅ ⋅
− ⋅ ⋅ ⋅
− ⋅ ⋅ ⋅ ⋅ −
=
+ ⋅
+ ⋅
+
+
⋅
+
⋅
+
⋅
+
+
⋅
2
1
2
2
3
4
2
1
2
2
3
4
2
1
1
π
π
π
π
π
π
π
⋯
⋯
(25)
lub teŜ do postaci iloczynowej z zapisanymi wprost zerami i biegunami:
(
) (
) (
)
(
) (
) (
)
H e
K
e
z
e
z
e
z
e
p
e
p
e
p
j
f
m
j
f
j
f
j
f
M
j
f
j
f
j
f
K
(
)
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
−
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
−
=
⋅
−
⋅
−
⋅ ⋅
−
−
⋅
−
⋅ ⋅
−
2
2
1
2
2
2
1
2
1
2
2
2
1
π
π
π
π
π
π
π
⋯
⋯
(26)
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
7
1.3 Transformacja Fouriera z “czasem” dyskretnym
Odpowiedź impulsowa jest sygnałem cyfrowym opisującym system liniowy i stacjonarny. KaŜdy sygnał
cyfrowy moŜna więc zinterpretować jako odpowiedź impulsową odpowiedniego systemu. Dla takiego
systemu moŜliwe jest wyznaczenie - o ile odpowiedni tylko zachodzi zbieŜność transformaty “Z” dla
odpowiedzi impulsowej na okręgu jednostkowym - odpowiedzi częstotliwościowej. MoŜna więc
przyjąć, Ŝe dla niemal dowolnego sygnału cyfrowego
x n
[ ]
- sygnał ten ma umoŜliwiać zbieŜność
transformaty “Z” dla okręgu jednostkowego (23) - określa się funkcję
X e
j
f
(
)
⋅ ⋅ ⋅
2
π
w analogiczny
sposób, jak w poprzednich podrozdziałach wyznaczano odpowiedź częstotliwościową. ZaleŜność
określająca związek pomiędzy
x n
[ ]
oraz
X e
j
f
(
)
⋅ ⋅ ⋅
2
π
nosi nazwę Transformacji Fouriera z Czasem
Dyskretnym i jest często oznaczana skrótowo jako D-TFT (ang. Discrete-Time Fourier Transform).
Nazwa ta bierze się stąd, Ŝe transformacja wykazuje istotne podobieństwo do Transformacji Fouriera
znanej z teorii sygnałów, jednak w tym przypadku czas (umowny) jest reprezentowany w postaci
indeksów
n
, a tylko częstotliwość jest wielkością ciągłą. Wzór definiujący D-TFT jest następujący:
(
)
X f
x n
j
f n
n
( )
[ ] exp
=
⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
∑
2
π
(27)
lub w zapisie bazującym na pulsacji cyfrowej:
(
)
X
x n
j
n
n
( )
[ ] exp
ω
ω
=
⋅
− ⋅ ⋅
=−∞
∞
∑
(28)
W zaleŜnościach (27) i (28) powrócono do zapisu wskazującego na
f
lub
ω
jako na zmienną
niezaleŜną, poniewaŜ transformatę
X f
( )
(lub
X ( )
ω
) przedstawia się zazwyczaj w postaci pary
wykresów amplitudy i fazy lub teŜ części rzeczywistej i urojonej, gdy oś pozioma opisana jest za
pomocą zmiennej
f
(lub odpowiednio
ω
).
Transformacja odwrotna do (27) określona jest wzorem:
(
)
x n
X f
j
f n df
[ ]
( ) exp
=
⋅
⋅ ⋅ ⋅ ⋅ ⋅
−
∫
2
1
2
1
2
π
(29)
natomiast transformacja odwrotna do (28), to:
(
)
x n
X
j
n d
[ ]
( ) exp
=
⋅
⋅
⋅ ⋅ ⋅
−
∫
1
2
π
ω
ω
ω
π
π
(30)
Transformacja odwrotna polega na wyznaczaniu całki w ramach jednego okresu transformaty, który - ze
względu na jej okresowość - zawiera kompletną informację o ciągu pierwotnym
x n
[ ]
. Wyprowadzenie
wzoru (29) moŜna znaleźć w dodatku teoretycznym nr 2.
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
8
1.3.1 Porównanie Transformacji Fouriera z czasem dyskretnym z Transformacją Fouriera dla
sygnałów z czasem ciągłym
W dodatku teoretycznym moŜna równieŜ znaleźć obszerną dyskusję na temat związków D-TFT z
Transformacją Fouriera dla sygnałów z czasem ciągłym, czyli C-TFT (ang. Continuous-Time Fourier
Transform). W skrócie związki te moŜna opisać następująco.
RozwaŜając dany ciąg
x n
[ ]
zawsze moŜna załoŜyć, Ŝe jest to sygnał cyfrowy pochodzący
z próbkowania z okresem
T
odpowiedniego sygnału ciągłego
x t
( )
. Próbkowaniem nazywa się tutaj
pobieranie wartości chwilowych
x t
( )
w odstępach
∆
t
T
=
. Określa się w ten sposób parę - sygnał
x t
( )
oraz ciąg
x n
[ ]
. MoŜna równieŜ załoŜyć, Ŝe sygnał
x t
( )
został spróbkowany funkcją
grzebieniową o takim samym okresie
T
i w wyniku tego powstał sygnał
x t
T
( )
, będący funkcją
grzebieniową zmodulowaną sygnałem
x t
( )
. Zagadnienie takiego próbkowania było obszernie opisane
w dodatku teoretycznym nr 1. Po sporządzeniu wykresów transformat moŜna stwierdzić, Ŝe
Transformata Fouriera typu C-TFT dla sygnału
x
t
T
( )
oraz Transformata Furiera typu D-TFT dla ciągu
x n
[ ]
mają identyczny kształt, pod warunkiem, Ŝe dopasuje się odpowiednio skale częstotliwości:
v
(wyraŜanej w Hz) dla C-TFT i częstotliwości
f
(bezwymiarowej) dla D-TFT. Dopasowanie polega na
tym, Ŝe odcinkowi
v
o długości
v
T
T
=
1
odpowiada odcinek
f
o długości 1. Przez
v
T
oznaczono
częstotliwość próbkowania. Rys. 2 stanowi ilustrację opisanej zaleŜności dla wykresów amplitud
odpowiednich transformat. Warto pamiętać, Ŝe delta Diraca jest pseudo-funkcją określoną dla czasu
ciągłego, zatem i funkcja grzebieniowa, a takŜe sygnał spróbkowany taką funkcją - czyli
x
t
T
( )
- teŜ
są określone dla czasu ciągłego. Ciąg
x n
[ ]
natomiast jest określony wyłącznie dla całkowitych
wartości
n
.
Rys.2. Ilustracja zaleŜności pomiędzy odpowiednimi transformatami - wykresy przedstawiają tylko
amplitudę - dla (kolejno od góry):
a)
x t
( )
- czyli sygnału z czasem ciągłym,
b)
x
t
T
( )
- tego samego sygnału, ale po spróbkowaniu funkcją grzebieniową,
c)
x n
[ ]
- sygnału cyfrowego pochodzącego z pobrania wartości chwilowych sygnału
x t
( )
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
9
Przykładowo, jeŜeli częstotliwość próbkowania wynosi
v
T
=
44
kHz, to wartość transformaty D-TFT
dla częstotliwości cyfrowej
f
=
0 25
,
jest taka sama, jak wartość transformaty C-TFT sygnału
x
t
T
( )
wyznaczona dla częstotliwości
v
=
11
kHz wyznaczona. Z kolei, jeŜeli zostały spełnione
warunki twierdzenia o próbkowaniu, to wartość transfomaty C-TFT od
x
t
T
( )
dla tej częstotliwości jest
T
kHz
=
1
44
mniejsza niŜ wartość transformaty C-TFT od
x t
( )
. Zgodność ta wynika to z faktu, Ŝe
wartość
v
=
11
kHz mieści się w odcinku od
− ⋅
1
2
v
T
do
+ ⋅
1
2
v
T
.
1.4 Dyskretna Transformacja Fouriera
W przypadku D-TFT czas jest reprezentowany przez indeksy przebiegające wartości od minus do plus
nieskończoności. Jednak w zagadnieniach praktycznych zakres indeksów czasowych jest ograniczony
do pewnego skończonego przedziału wartości. W takiej sytuacji moŜliwe jest stosowanie Dyskretnej
Transformacji Fouriera, wskrócie DFT (ang. Discrete Fourier Transform). Wprowadza się w tym celu
nowe indeksy, tak by pierwszy z nich posiadał wartość zero. Nie ogranicza to w Ŝadnej mierze
interpretacji otrzymanych wyników w odniesieniu do oryginalnego zakresu indeksów. JeŜeli,
przykładowo, rozwaŜany jest ciąg
x n
0
[ ]
dla
n
=< −
>
10
100
,
, to dla potrzeb DFT moŜna
wprowadzić przesunięcie indeksów o 10 i otrzymać ciąg
x n
[ ]
dla
n
=<
>
0
110
,
. Związek
pomiędzy
x n
0
[ ]
oraz
x n
[ ]
jest oczywisty, bowiem:
x n
x n
0
10
[
]
[ ]
−
=
dla
n
=<
>
0
110
,
.
DFT wykazuje szereg podobieństw do D-TFT. JeŜeli załoŜymy, Ŝe w ciągu nieskończonym z wyjątkiem
N
kolejnych elementów wszystkie pozostałe mają wartość zero, to po wprowadzeniu nowych
indeksów zakres sumowania we wzorze na D-TFT (22) moŜna następująco ograniczyć:
(
)
X f
x n
j
f n
n
N
( )
[ ] exp
=
⋅
− ⋅ ⋅ ⋅ ⋅
=
−
∑
2
0
1
π
(31)
JeŜeli z otrzymanej według (31) transformaty zostanie pozostawiony tylko jeden okres - na przedziale
f
od zera do 1 - i na tym odcinku zostanie pobranych
N
próbek, w oddalonych o
∆
f
N
=
1
punktach, poczynając od punktu
f
=
0
, to otrzyma się
N
wartości. Punkty próbkowania
f
k
są
rozłoŜone na osi
f
równomiernie, moŜna je więc we wzorze (31) zastąpić ich numerami,
wprowadzając indeksy częstotliwościowe:
k
f
N
k
=
⋅
(32)
Otrzymany według powyŜszego opisu ciąg
N
próbek transformaty D-TFT moŜna zatem zapisać
następująco:
X k
x n
j
k n
N
k
N
n
N
[ ]
[ ] exp
:
, ,...,
=
⋅
− ⋅ ⋅ ⋅ ⋅
=
−
=
−
∑
2
0 1
1
0
1
π
(33)
ZaleŜność (33) jest wzorem Dyskretnej Transformacji Fouriera, czyli DFT. Wzór na odwrotną DFT jest
do (33) bardzo podobny:
x n
N
X k
j
k n
N
n
N
k
N
[ ]
[ ] exp
:
, ,...,
=
⋅
⋅
⋅ ⋅ ⋅ ⋅
=
−
=
−
∑
1
2
0 1
1
0
1
π
(34)
Para wzorów (33) i (34) stanowi przykład wzoru na DFT zaimplementowany w Matlab’ie. Ze względu
na róŜne zastosowania parę zaleŜności na DFT oraz na odwrotną DFT moŜna ogólnie zapisać z
uŜyciem stałej
c
DFT
:
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
10
X k
c
x n
j
k n
N
k
N
DFT
n
N
[ ]
[ ] exp
:
, ,...,
=
⋅
⋅
− ⋅ ⋅ ⋅ ⋅
=
−
=
−
∑
2
0 1
1
0
1
π
(35)
x n
N c
X k
j
k n
N
n
N
DFT k
N
[ ]
[ ] exp
:
, ,...,
=
⋅
⋅
⋅ ⋅ ⋅ ⋅
=
−
=
−
∑
1
2
0 1
1
0
1
π
(36)
Najczęściej stosuje się wartość
c
DFT
ze zbioru
1
1
1
N
N
,
,
. W pierwszym przypadku
X[ ]
0
jest
wartością średnią wszystkich elementów ciągu
x n
[ ]
. W drugim, w wyniku transformacji nie zmienia
się całkowita energia ciągu. Właściwość tą, zwaną równością Parsevala, moŜna zapisać tak:
x n
X k
n
N
k
N
[ ]
[ ]
2
0
1
2
0
1
=
−
=
−
∑
∑
=
(37)
Natomiast, gdy
c
DFT
=
1
, to
X[ ]
0
jest sumą wszystkich elementów ciągu pierwotnego.
Równość Parsevala moŜna zapisać bardziej ogólnie, uwzględniając stałą
c
DFT
:
(
)
c
N
x n
X k
DFT
n
N
k
N
⋅
⋅
=
=
−
=
−
∑
∑
2
2
0
1
2
0
1
[ ]
[ ]
(38)
Transformacja DFT przekształca pierwotny ciąg
N
elementów (zwykle jest to ciąg rzeczywisty)
w ciąg transformaty (zazwyczaj zespolony) składający się równieŜ z
N
elementów. Warto zaznaczyć,
Ŝ
e z punktu widzenia definicji DFT ciąg pierwotny moŜe być takŜe zespolony, co zostało uwzględnione
w zapisie (37) i (38).
1.5 Transformata DFT jako próbkowanie transformaty D-TFT
Jak wynika z poprzedniego podrozdziału, pomiędzy DFT i D-TFT istnieje bezpośredni związek. Rys.3
przedstawia ilustrację tej zaleŜności.
Rys.3. Ilustracja faktu, Ŝe transformata DFT moŜe być interpretowana jako efekt próbkowania jednego
okresu transformaty D-TFT, gdy transformata D-TFT jest wyznaczona dla tego samego skończonego
ciągu jednak przedłuŜonego do nieskończoności elementami zerowymi
Korzystając z relacji między transformatą C-TFT dla sygnału z czasem ciągłym i D-TFT ciągu
pochodzącego z próbkowania tego sygnału oraz informacji podanych powyŜej moŜna określić związek
pomiędzy widmem Fouriera sygnału z czasem ciągłym i DFT odpowiedniego fragmentu ciągu
cyfrowego pochodzącego z próbkowania tego sygnału. W ramach przygotowywania się do ćwiczeń
naleŜy przygotować się do określania, jaka częstotliwość w Hz odpowiada określonemu indeksowi
częstotliwościowemu
k
i na odwrót.
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
11
2 Korzystanie z pakietu MATLAB
2.1 Opis wybranych funkcji
Wybrane polecenia programu MATLAB:
fft( dane wejściowe )
wyznaczenie dyskretnej transformaty Fouriera zadanego
ciągu za pomocą algorytmu szybkiej transformaty
Fouriera
ifft( dane wejściowe )
operacja odwrotna do fft
fftshift( dane wejściowe )
przesunięcie w dziedzinie Fouriera, tak by składowa
stała znalazła się (prawie) w środku ciągu
sum( wektor lub macierz )
wyznaczanie sumy elementów wektora lub, gdy zostanie
podana macierz, wektora sum kolejnych kolumn tej
macierzy
flops
wyświetlenie stanu licznika operacji
zmiennoprzecinkowych, flops(0) zeruje ten licznik
mesh
wykreślanie w aktywnym oknie graficznym wykresu
funkcji dwóch zmiennych podanej w postaci
zdyskretyzowanej jako macierz (wykres trójwymiarowy)
colormap
zmiana palety kolorów uŜywanej do wykreślania
wykresów
view
obrót wykresu trójwymiarowego wokół osi pionowej i
poziomej do pozycji podanej w postaci kątów azymutu i
elewacji
Szczegóły dotyczące formatu danych wejściowych naleŜy odczytać za pomocą polecenia „help”.
W przypadku funkcji fft oraz ifft moŜliwe jest wyznaczanie transformaty DFT dla ilości elementów
innej niŜ długość ciągu podanego jako dane wejściowe. Wystarczy w tym celu podać drugą zmienną,
określającą ilość punktów. Przykładowo dla 8-punktowego ciągu
x n
[ ]
moŜna wyliczyć 32-punktową
transformatę DFT przez podanie polecenia:
>>X=fft(x,32);
W przypadku, gdy ilość punktów transformaty jest większa niŜ długość ciągu pierwotnego, to ciąg ten
jest przedłuŜany do odpowiedniej długości wartościami zerowymi. JeŜeli ilość punktów transformaty
jest mniejsza niŜ długość ciągu, to ciąg ten jest przy wyznaczaniu transformaty obcinany do wskazanej
długości.
Uwaga: jeŜeli zamieniamy zespolony wektor wierszowy na kolumnowy (lub odwrotnie) stosując
operator transpozycji macierzy, naleŜy pamiętać, Ŝe otrzymane wartości będą miały odwrócony znak
części urojonej - wynika to z definicji transpozycji dla macierzy zespolonych. MoŜna to łatwo
„poprawić” stosując polecenie „conj”.
2.2 Przykłady realizacji wybranych zadań
PoniŜej pokazano przykładowe sekwencje poleceń realizujące zadania podobne do tych, które będą
realizowane w ramach laboratorium - w większości poniŜszych przykładów pominięto końcowy etap
wykreślania odpowiednich wykresów.
2.2.1 Demonstracja okresowości sygnału kosinusoidalnego względem “czasu” i częstotliwości:
>> t=0:0.02:3;
>> v=0:0.02:2:
>> s=cos(2*pi*v’*t);
>> mesh(t,v,s);
>> colormap([0,0,0]);
>> axis([min(t),max(t),min(v),max(v),-1,1]);
>> grid
>> view([75,65]);
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
12
2.2.2 Wyznaczanie odpowiedzi częstotliwościowej filtru FIR rzędu N-1 dla pojedynczej wartości
częstotliwości cyfrowej:
>> N=8;
>> h=rand(1,N);
>> f=0.25;
>> n=0:N-1;
>> s_we=exp(j*2*pi*f*n);
>> s_wy=filter(h,1,s_we);
>> H=s_wy/s_we;
W przykładzie wyznaczono odpowiedź systemu tylko dla ciągów 8-punktowych, poniewaŜ w 8.
punkcie odpowiedzi zanika róŜnica pomiędzy odpowiedzią na odpowiedni nieskończony ciąg i na ciąg
skończony zastosowany w przykładzie. Wynika to z faktu, iŜ odpowiedź impulsowa jest właśnie 8-
elementowa.
2.2.3 Wyznaczanie transmitancji “Z” filtru typu FIR o zadanej odpowiedzi impulsowej dla
pojedynczej wartości z:
>> f=0.5;
>> z=exp(j*2*pi*f)
>> h=rand(1,8);
>> p=0:length(h)-1;
>> H=sum(h.*(z.^p));
2.2.4 Wyznaczanie przybliŜenia transformaty D-TFT ze wzoru definicyjnego (27):
>> f=-30:0.01:30;
>> T=1/20;
>> x=rand(1,32);
>> n=-10:21;
i następnie:
>> k=0; for ff=f, k=k+1; X(k)=sum(x.*exp(-j*2*pi*f*n*T)); end;
albo:
>> k=(1:length(f))’;
>> X(k)=exp(-j*2*pi*f’*n*T)*x’;
2.2.5 Wyznaczanie transformaty DFT ze wzoru definicyjnego (33):
>> N=32;
>> n=0:31;
>> x=rand(1,N);
i następnie:
>> for k=n; X(k+1)=sum(x.*exp(-j*2*pi*k*n/N)); end;
albo:
>> k=n’;
>> X(k+1)=exp(-j*2*pi*k*n/N)*x’;
2.2.6 Wyznaczanie przybliŜenia transformaty D-TFT za pomocą funkcji fft:
>> N=8;
>> Nmax=1024;
>> x=rand(1,N);
>> X=fft(x,Nmax);
2.2.7 Wykres fazy - efekt ograniczania do zakresu od minus pi do pi:
>> fi=(-10:0.1:10)*pi;
>> X=exp(j*fi);
>> stem( fi, angle(X));
>> stem(fi, unwrap(angle(X)));
Jak uzasadnić wprowadzenie demonstrowanego powyŜej ograniczenia na zakres wartości fazy?
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
13
D
odatek teoretyczny nr 2
Relacja pomiędzy Transformacją Fouriera z Czasem Dyskretnym (ang.
D-TFT) i Transformacją Fouriera z Czasem Ciągłym (ang. C-TFT)
Definicja D-TFT dla ciągu
x n
[ ]
:
X f
x n e
f
j
f n
n
( )
[ ]
:
=
⋅
∈ℜ
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
∑
2
π
(1)
Definicja C-TFT dla sygnału
x t
( )
(dotyczy takich sygnałów, dla których poniŜsza całka jest zbieŜna -
dla sygnałów okresowych stosuje się pewne poszerzenie definicji, do znaku trzech gwiazdek ***
rozwaŜania dotyczyć będą wyłącznie sygnałów gwarantujących zbieŜność wzoru definicyjnego):
Ψ
( )
( )
:
v
x t
e
dt
v
j
v t
=
⋅
∈ℜ
− ⋅ ⋅ ⋅ ⋅
−∞
∞
∫
2
π
(2)
Definicja transformacji odwrotnej do C-TFT (z uwagami odnośnie zbieŜności analogicznymi do wzoru
(2)):
x t
v
e
dv
t
j
v t
( )
( )
:
=
⋅
∈ℜ
⋅ ⋅ ⋅ ⋅
−∞
∞
∫
Ψ
2
π
(2a)
ZałóŜmy, Ŝe wartości ciągu
x n
[ ]
są równe amplitudom delt Diraca po próbkowaniu sygnału
x t
( )
w
równych odstępach czasu T, czyli po prostu chwilowym wartościom sygnału w równo oddalonych
punktach czasu:
x n
x n T
[ ]
(
)
=
⋅
(3)
Dla przypomnienia, funkcja grzebieniowa to ciąg równo oddalonych delt Diraca (w tym przypadku
oddalonych o T):
g
t
t
n T
T
n
( )
(
)
=
− ⋅
=−∞
∞
∑
δ
(4)
Zatem spróbkowany funkcją grzebieniową sygnał
x t
( )
moŜna zapisać następująco (inaczej jest to
funkcja grzebieniowa zmodulowana amplitudowo sygnałem
x t
( )
):
x
t
x t
g
t
T
T
( )
( )
( )
=
⋅
(5)
Wynik C-TFT dla sygnału (5) określony jest wzorem:
Ψ
T
T
j
v t
v
x t
g t
e
dt
( )
( )
( )
=
⋅
⋅
− ⋅ ⋅ ⋅ ⋅
−∞
∞
∫
2
π
(6)
Skorzystajmy z (4) i podstawmy za funkcję grzebieniową do (6):
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
14
Ψ
T
n
j
v t
v
x t
t
n T
e
dt
( )
( )
(
)
=
⋅
− ⋅
⋅
=−∞
∞
− ⋅ ⋅ ⋅ ⋅
−∞
∞
∑
∫
δ
π
2
(7)
PoniewaŜ w wyraŜeniu podcałkowym występuje iloczyn wyraŜeń, skorzystamy z przemienności
mnoŜenia oraz rozdzielności mnoŜenia względem dodawania i przepiszemy (7) do postaci:
[
]
Ψ
T
j
v t
n
v
x t
e
t
n T
dt
( )
( )
(
)
=
⋅
⋅
− ⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
−∞
∞
∑
∫
2
π
δ
(8)
Zakładając, Ŝe całka (7) jest zbieŜna dla kaŜdego składnika sumy moŜemy całkę sumy wyraŜeń zastąpić
sumą całek:
Ψ
T
j
v t
n
v
x t
e
t
n T dt
( )
( )
(
)
=
⋅
⋅
− ⋅
− ⋅ ⋅ ⋅ ⋅
−∞
∞
=−∞
∞
∫
∑
2
π
δ
(9)
Warto teraz przypomnieć, Ŝe całka z funkcji pomnoŜonej przez deltę Diraca przesuniętą o
∆
T
jest
równa wartości tej funkcji w punkcie
∆
T
(dowód w tym miejscu pomijamy), czyli przykładowo:
p t
t
T dt
p
T
( )
(
)
(
)
−∞
∞
∫
⋅
−
=
δ
∆
∆
(10)
Traktując iloczyn
x t
e
j
v t
( )
⋅
− ⋅ ⋅ ⋅ ⋅
2
π
jako powyŜszą funkcję ciągłą moŜemy przepisać (9) do postaci:
Ψ
T
j
v n T
n
v
x n T
e
( )
(
)
=
⋅ ⋅
− ⋅ ⋅ ⋅ ⋅ ⋅
=−∞
∞
∑
2
π
(11)
WyraŜenie (11) jest juŜ bardzo podobne do (1), naleŜy jeszcze wykonać podstawienie (3) oraz:
v
f
T
=
(12)
Po powyŜszych podstawieniach otrzymujemy, Ŝe:
X f
f
T
T
( )
=
Ψ
(13)
Oznacza to, Ŝe funkcja ciągła częstotliwości f (D-TFT ciągu
x n
[ ]
) jest identyczna (z dokładnością do
skali T określonej zaleŜnością (12) dla zmiennej niezaleŜnej) z funkcją ciągłą zmiennej v (C-TFT
sygnału ciągłego ).
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
15
Wynika z tego szereg wniosków, np. taki:
Ψ
T
v
( )
jest transformatą C-TFT sygnału spróbkowanego z okresem T, zatem ta transformata jest
okresowa względem zmiennej v z okresem
v
T
0
1
=
, z (13) wynika więc, Ŝe transformata D-TFT,
X f
( )
, jest takŜe okresowa względem zmiennej
f
z okresem
f
0
1
=
.
Spotykane są takŜe definicje obu transformat oparte na pulsacji. Otrzymuje się je przez następujące
podstawienia.
Dla D-TFT:
ω
π
= ⋅ ⋅
2
f
(14)
oraz dla C-TFT:
Ω = ⋅ ⋅
2
π
v
(15)
Pomiędzy pulsacjami (14) i (15) zachodzi naturalnie analogiczny do (12) związek:
Ω =
ω
T
(16)
Przy podstawieniach (14) i (15) moŜna wykazać, Ŝe zachodzi takŜe odpowiedni związek pomiędzy
transformatą D-TFT (czyli w tym przypadku:
X ( )
ω
) oraz transformatą C-TFT (czyli teraz:
X
T
( )
Ω
)
:
X
T
T
( )
ω
ω
=
Ψ
(17)
Funkcja
X ( )
ω
jest okresowa względem zmiennej
ω
i jak łatwo z (14) zauwaŜyć okres ten wynosi
2
⋅
π
.
Warto takŜe wskazać inną moŜliwość przejścia od (8) do (11).
[
]
Ψ
T
j
v t
n
v
x t
e
t
n T
dt
( )
( )
(
)
=
⋅
⋅
− ⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
−∞
∞
∑
∫
2
π
δ
(8 - powtórzone)
Zakładamy, Ŝe całka (8) jest zbieŜna w kaŜdym podprzedziale całkowania i następnie dzielimy cały
przedział całkowania na odcinki o długości T i środkach w punktach określonych dla całkowitych k
następująco:
t
k T
k
= ⋅
(18)
PoniewaŜ suma całek na wszystkich podprzedziałach jest równa całce na pełnym przedziale
całkowania, więc (8) moŜna przepisać do postaci:
[
]
Ψ
T
j
v t
n
t
T
t
T
k
v
x t
e
t
n T
dt
k
k
( )
( )
(
)
=
⋅
⋅
− ⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
−
+
=−∞
∞
∑
∫
∑
2
2
2
π
δ
(19)
W danym podprzedziale całkowania mieści się tylko jedna odpowiednio przesunięta delta Diraca,
zatem suma pod całką moŜe być zastąpiona przez pojedynczy składnik:
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
16
Ψ
T
j
v t
t
T
t
T
k
v
x t
e
t
k T dt
k
k
( )
( )
(
)
=
⋅
⋅
− ⋅
− ⋅ ⋅ ⋅ ⋅
−
+
=−∞
∞
∫
∑
2
2
2
π
δ
(20)
Teraz korzystając z (10) i podstawiając n za k moŜna łatwo dokończyć wykazanie równowaŜności (13).
Opierając się na wykazanej zaleŜności pomiędzy D-TFT i C-TFT moŜna wywnioskować, Ŝe
transformacja odwrotna do (1) będzie dana wzorem wynikającym z (2a). PoniewaŜ D-TFT jest
powiązane ściśle z C-TFT sygnału ciągłego w czasie spróbkowanego funkcją grzebieniową, zatem
wypiszemy wzór na transformację odwrotną do C-TFT dla tego właśnie przypadku:
x t
v
e
dv
t
T
T
j
v t
( )
( )
:
=
⋅
∈ℜ
⋅ ⋅ ⋅ ⋅
−∞
∞
∫
Ψ
2
π
(21)
Dla przypomnienia - ogólny wzór na odwrotną transformację jest następujący:
x t
v
e
dv
t
j
v t
( )
( )
:
=
⋅
∈ℜ
⋅ ⋅ ⋅ ⋅
−∞
∞
∫
Ψ
2
π
(2a - powtórzone)
Z transformaty odwrotnej wyznaczonej według (21) wybieramy teraz jedynie te wartości
x t
( )
, które
odpowiadają elementom ciągu
x n
[ ]
, zatem (21) przybiera postać:
x n T
v
e
dv
n
I
T
T
j
v n T
(
)
( )
:
⋅
=
⋅
∈
⋅ ⋅ ⋅ ⋅ ⋅
−∞
∞
∫
Ψ
2
π
(22)
Z kolei transformata odwrotna dla sygnału niespróbkowanego wyznaczona jedynie w tych samych,
wybranych punktach czasu opisana jest w sposób następujący:
x n T
v
e
dv
n
I
j
v n T
(
)
( )
:
⋅
=
⋅
∈
⋅ ⋅ ⋅ ⋅ ⋅
−∞
∞
∫
Ψ
2
π
(23)
Pomiędzy wynikami obu transformat odwrotnych (22) i (23) w tych wybranych punktach czasu
zachodzi następujący związek (wynikający ze sposobu próbkowania):
x
n T
x n T
t
n T
T
(
)
(
)
(
)
⋅
=
⋅
⋅
− ⋅
δ
(24)
DąŜymy zatem do tego, by w oparciu o
Ψ
T
v
( )
(będące w ścisłym związku z
X f
( )
) wyznaczyć
wartości
x n
[ ]
(które są z kolei odpowiednio równe wartościom
x n T
(
)
⋅
). W tym celu naleŜy
doprowadzić do równości zmodyfikowanej odpowiednio prawej strony (22) i lewej strony (23).
Natępnie wystarczy wykorzystać (12) oraz (13) i zastąpić po prawej stronie
Ψ
T
v
( )
przez
X f
( )
.
Przyjmijmy, Ŝe sygnał ciągły przed próbkowaniem ma taką transformatę, Ŝe jego częstotliwość
maksymalna i częstotliwość próbkowania spełniają warunek o próbkowaniu umoŜliwiający wierne
odtworzenie. ZauwaŜmy, Ŝe moŜemy to zawsze załoŜyć, poniewaŜ nie interesuje nas w zasadzie jak
wyglądał sygnał ciągły
x t
( )
- wykonujemy róŜne operacje jedynie na sygnale spróbkowanym
x
t
T
( )
oraz na wybranych wartościach chwilowych
x t
( )
. To jak wygląda sygnał
x t
( )
pomiędzy tymi
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
17
wybranymi punktami nie ma Ŝadnego znaczenia z punktu widzenia przedstawionej dyskusji. W kaŜdym
razie, gdybyśmy chcieli później wyciągnąć w oparciu o odwrotną transformację D-TFT jakieś wnioski
odnośnie całego
x t
( )
, to sygnałów takich moŜe być dla konkretnego ciągu
x n
[ ]
nieskończenie wiele.
Zawsze jednak moŜna zastosować interpolację prowadzącą do takiego
x t
( )
, który spełnia warunek o
próbkowaniu.
JeŜeli zatem
x t
( )
spełnia warunek o próbkowaniu, to jego transformata
Ψ
( )
v
jest związana z
transformatą sygnału spróbkowanego ,
Ψ
T
v
( )
, następująco:
Ψ
Ψ
( )
( )
,
.
v
T
v
dla
v
v
v
dla
pozost v
T
=
⋅
∈ −
0
0
2
2
0
(25)
W celu wyznaczenia wartości
x t
( )
w wybranych punktach
t
n T
= ⋅
wystarczy zatem wyznaczyć
transformatę odwrotną do C-TFT dla
T
v
T
⋅ Ψ
( )
na przedziale zawierającym pojedynczy okres - od
−
v
0
2
do
v
0
2
:
x n T
T
v
e
dv
n
I
T
j
v n T
v
v
(
)
( )
:
⋅
= ⋅
⋅
∈
⋅ ⋅ ⋅ ⋅ ⋅
−
∫
Ψ
2
2
2
0
0
π
(26)
Podstawiając według (12) i (13) oraz (3), doprowadzając między innymi do zmiany zmiennej po której
wyznaczana jest całka a zatem i do zmiany granic całkowania, otrzymujemy:
x n T
T
e
d
f
T
n
I
T
f
T
j
f
T
n T
f
T
f
T
(
)
( )
:
⋅
= ⋅
⋅
∈
⋅ ⋅ ⋅ ⋅ ⋅
−
⋅
⋅
∫
Ψ
2
2
2
0
0
π
(27)
czyli (warto przypomnieć, Ŝe zawsze
f
0
1
=
):
x n
X f
e
df
n
I
j
f n
[ ]
( )
:
=
⋅
∈
⋅ ⋅ ⋅ ⋅
−
∫
2
1
2
1
2
π
(28)
Wzór (28) stanowi więc przepis na transformację odwrotną do D-TFT.
Jak łatwo zauwaŜyć, dla sygnałów o których była dotychczas mowa (czyli spróbkowanych funkcją
grzebieniową) całka (2a) nie będzie zbieŜna dla Ŝadnej wartości
t
poniewaŜ funkcja
Ψ
T
v
( )
jest
okresowa. Podobny problem pojawia się takŜe w przypadku transformacji w przód dla funkcji
okresowej - zarówno C-TFT jak i D-TFT. Konieczne jest zatem poszerzenie definicji (1) i (2) - oraz
automatycznie (2a) i (28).
***
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
18
Uzupełnienie dyskusji o przypadek dla sygnałów okresowych (czy to w dziedzinie
czasu, czy częstotliwości).
Poszerzenie definicji transformacji polega na przyjęciu, Ŝe jeŜeli para funkcji
x t
( )
oraz
Ψ
( )
v
jest
powiązana jednym z dwóch wzorów całkowych (albo na transformatę w przód albo odwrotną), to
zachodzi równieŜ transformacja w drugą stronę, nawet, gdy całka w odpowiadającym wzorze nie jest
zbieŜna. Przykładowo, jeŜeli odwrotna transformata C-TFT (wzór całkowy (2a)) z przesuniętej delty
Diraca jest zbieŜna i daje zespoloną eksponentę, to takŜe transformata C-TFT z zespolonej eksponenty
(będącej sygnałem okresowym) istnieje i jest równa tej delcie Diraca:
[
]
[
]
F
v
v
e
F
e
v
v
C T
j
v t
C T
j
v t
−
−
⋅ ⋅ ⋅ ⋅
−
⋅ ⋅ ⋅ ⋅
−
=
⇔
=
−
1
2
2
δ
δ
π
π
(
)
(
)
∆
∆
∆
∆
(29)
Rozpoczniemy od okresowego sygnału ciągłego w dziedzinie czasu (podstawowy okres sygnału wynosi
P). Sygnał okresowy moŜna rozwinąć w szereg Fouriera, zatem:
x
t
c
e
P
k
j
k
P
t
k
( )
=
⋅
⋅ ⋅ ⋅ ⋅
=−∞
∞
∑
2
π
(30)
Z liniowości transformacji C-TFT (2) oraz (29) wynika, Ŝe transformata sygnału (30) jest ciągiem
poprzesuwanych delt Diraca, skalowanych amplitudowo współczynnikami rozwinięcia w szereg:
[
]
F
x
t
X
v
c
v
k v
C T
P
P
k
P
k
−
=−∞
∞
=
=
⋅
− ⋅
∑
( )
( )
(
)
δ
(31)
gdzie
v
P
P
=
1
(32)
Zatem transformata ta jest funkcją prąŜkową.
Określony powyŜej sygnał o okresie P po spróbkowaniu z krokiem próbkowania T oznaczymy
następująco:
x
t
x
t
g
t
P T
P
T
,
( )
( )
( )
=
⋅
(33)
Ze względu na rozdzielność dodawania względem mnoŜenia wzór (33) moŜna rozpisać na sumę
próbkowanych z osobna wyrazów szeregu (30). Pojedynczy, przykładowo k-ty, wyraz szeregu (30) po
spróbkowaniu opisany jest następująco:
x
t
x
t
g
t
c
e
g
t
P T k
P k
T
k
j
k t
P
T
, ,
,
( )
( )
( )
( )
=
⋅
=
⋅
⋅
⋅
⋅ ⋅ ⋅
2
π
(34)
Korzystając z tego, Ŝe mnoŜenie w dziedzinie czasu odpowiada splotowi w dziedzinie częstotliwości,
moŜna opisać transformatę C-TFT sygnału (34) jako:
[
]
(
)
F
x
t
c
v
k v
T
g
v
C T
P T k
k
P
v
−
=
⋅
− ⋅
∗
⋅
, ,
( )
(
)
( )
δ
1
0
(35)
PoniewaŜ splot z przesuniętą deltą Diraca daje przesunięcie funkcji o tyle samo o ile przesunięta jest
delta (a funkcja grzebieniowa jest sumą poprzesuwanych delt Diraca), więc (35) opisuje deltę Diraca w
dziedzinie częstotliwości
v
, przesuniętą o
k v
P
⋅
w prawo, o amplitudzie
c
T
k
i powielonej co
v
0
wzdłuŜ całej osi częstotliwości.
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
19
Porównując powyŜsze wzory moŜna zatem zauwaŜyć, Ŝe transformata C-TFT sygnału okresowego i
próbkowanego funkcją grzebieniową jest okresową funkcją prąŜkową. PrąŜkowa transformata sygnału
okresowego nie próbkowanego jest tym razem pomnoŜona przez współczynnik 1/T i powielona wzdłuŜ
osi częstotliwości co
v
0
. Jest to dokładna analogia do transformaty spróbkowanego sygnału
nieokresowego, tyle Ŝe tym razem transformata jest prąŜkowa (suma poprzesuwanych delt Diraca).
Warto zauwaŜyć, Ŝe zazwyczaj powinno się zapewnić warunek:
v
v
P
<<
0
. RównieŜ i tego przypadku
dotyczy związek pomiędzy częstotliwością próbkowania i maksymalną częstotliwością sygnału (to, Ŝe
transformata jest prąŜkowa nie oznacza, iŜ sygnał nie moŜe mieć częstotliwości maksymalnej).
Okazuje się, Ŝe w takim przypadku (gdy sygnał jest okresowy), to D-TFT wiąŜe się z C-TFT w
identyczny sposób jak w (13) (lub (17)). Suma (1) nie jest skończona i dlatego do wyznaczenia
transformaty naleŜy posłuŜyć się poszerzoną definicją i wiedzą o tym, Ŝe według odwrotnej D-TFT (28)
z odpowiedniego ciągu delt Diraca w dziedzinie częstotliwości otrzyma się ciąg wartości
odpowiadający wartościom chwilowym ciągłego sygnału okresowego.
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
20
D
odatek teoretyczny nr 3
Relacja pomiędzy Transformacją Fouriera z Czasem Dyskretnym (ang.
D-TFT) i Dyskretną Transformacją Fouriera (ang. DFT)
Relacja pomiędzy DFT i C-TFT
Definicja D-TFT dla ciągu
x n
[ ]
:
X f
x n e
f
j
f n
n
( )
[ ]
:
=
⋅
∈ℜ
− ⋅ ⋅ ⋅ ⋅
=−∞
∞
∑
2
π
(1)
Definicja DFT dla ciągu
x n
[ ]
:
X k
x n e
k
N
j
k n
N
n
N
[ ]
[ ]
:
, ,...,
=
⋅
=
−
− ⋅⋅ ⋅ ⋅ ⋅
=
−
∑
2
0
1
0 1
1
π
(2)
Z definicji od razu wynika, Ŝe transformata D-TFT jest funkcją ciągła ze względu na
f
, natomiast
transformata DFT jest ciągiem o tej samej ilości elementów, co ciąg pierwotny - dyskretny jest i czas i
częstotliwość, a poniewaŜ określony jest odstęp pomiędzy kolejnymi punktami w obu dziedzinach, więc
do określenia zmiennej niezaleŜnej pozostawia się jedynie indeks (czasowy lub częstotliwościowy).
Ponadto transformata D-TFT określona jest dla nieskończonego ciągu
x n
[ ]
, natomiast DFT dla ciągu
skończonego.
Jaka jest zatem relacja pomiędzy obiema powyŜszymi transformatami?
Wartości elementów ciągu z indeksami czasowymi mogą być traktowane jako chwilowe wartości
odpowiedniej funkcji czasu wybrane w chwilach czasu oddalonych o T sekund (patrz relacja pomiędzy
D-TFT i C-TFT):
x n
x n T
[ ]
(
)
=
⋅
(3)
Taka interpretacja umoŜliwia powiązanie Transformacji Fouriera z Czasem Dyskretnym (ang. D-TFT) i
Transformacji Fouriera z Czasem Ciągłym (and. C-TFT). WykaŜemy teraz, Ŝe jeŜeli załoŜymy
analogiczną do (3) zaleŜność w dziedzinie częstotliwości, to transformata DFT stanowi ciąg N
wartości chwilowych ciągłej transformaty D-TFT wybranych z jednego okresu tej transformaty
(a jest ona, jak wiadomo, okresowa), oddalonych na osi częstotliwości
f
o
1
N
. Przepiszmy zatem
(2) wstawiając okres próbkowania w dziedzinie częstotliwości
f
N
N
=
1
(dla przypomnienia: okres
transformaty D-TFT wynosi
f
0
1
=
):
X k
x n e
k
N
j
k f
n
N f
n
N
N
N
[ ]
[ ]
:
, ,...,
=
⋅
=
−
− ⋅⋅ ⋅ ⋅ ⋅ ⋅
⋅
=
−
∑
2
0
1
0 1
1
π
(4)
W (4) moŜna zauwaŜyć, Ŝe:
N f
f
N
⋅
=
=
0
1
(5)
oraz, oznaczając przez
f
k
k-ty punkt na osi częstotliwości, Ŝe:
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
21
k f
f
N
k
⋅
=
(6)
Podstawiając (5) i (6) do (4) oraz zakładając, Ŝe ciąg
x n
[ ]
przyjmuje dla nie uwzględnionych
indeksów czasowych wartości zerowe, moŜna wykazać, Ŝe wzór (4) określa D-TFT (wzór (1)) dla
wybranych wartości
f
:
X k
x n e
X k
f
f
X k f
j
f n
f
n
N
k
k
k
[ ]
[ ]
(
)
=
⋅
=
⋅
=
⋅
− ⋅⋅ ⋅ ⋅ ⋅
=
−
∑
2
0
1
0
0
π
(7)
Prawa strona (7) oznacza wybrane wartości transformaty D-TFT. Pamiętając, Ŝe
k przyjmuje jedynie
całkowite wartości od 0 do
N-1, widzimy, Ŝe faktycznie DFT dotyczy tylko jednego okresu D-TFT.
Dotatkowo krok próbkowania w dziedzinie częstotliwości związany jest z długością interesującego nas
odcinka w dziedzinie czasu (dyskretnego) następująco:
f
N
N
=
1
(8)
lub dzieląc obie strony przez krok próbkowania
T w dziedzinie czasu (otrzymując w ten sposób wymiar
czasu w sekundach, a częstotliwości w Hz):
v
f
T
N T
t
N
N
=
=
⋅
=
1
1
∆
(8a)
Analogiczny związek zachodzi pomiędzy krokiem próbkowania w dziedzinie czasu i okresem w
dziedzinie częstotliwości. Dla czasu dyskretnego krok próbkowania wynosi po prostu „1”:
1
1
0
=
f
(9)
MnoŜąc obie strony (9) przez
T otrzymujemy:
T
T
T
f
v
= ⋅ =
=
1
1
0
0
(9a)
gdzie
T jest krokiem próbkowania osi czasu wyraŜonym w sekundach, a
v
0
okresem transformaty
C-TFT wyraŜonym w Hz.
MoŜna więc teraz takŜe odnieść DFT do C-TFT.
JeŜeli załoŜymy, Ŝe:
x n
dla
n
n
N
[ ]
lub
(
)
=
<
>
−
0
0
1
(10)
to moŜna zastosować następującą interpretację:
x n
x n T
X k
X k f
N
[ ]
(
)
[ ]
(
)
=
⋅
∧
=
⋅
(11)
Oznacza to, Ŝe wartości ciągu
x n
[ ]
są chwilowymi wartościami ciągłego sygnału
x t
( )
, dla chwil
czasu określonych przez zaleŜność:
DSP-MATLAB,
Ć
wiczenie 4,
P.Korohoda, KE AGH
22
t
n T
n
= ⋅
(12)
Natomiast wartości ciągu
X k
[ ]
(czyli w dziedzinie indeksów częstotliwościowych) stanowią wartości
chwilowe ciągłej funkcji
X f
( )
(wyniku D-TFT), a więc zarazem funkcji
Ψ
T
v
( )
(tj. wyniku C-TFT
dla spróbkowanej funkcją grzebieniową - z krokiem próbkowania T - funkcji
x t
( )
):
X k
X k f
k f
T
v
N
T
N
T
k
[ ]
(
)
(
)
=
⋅
=
⋅
=
Ψ
Ψ
(13)
Zatem
X k
[ ]
jest jednocześnie wartością
Ψ
T
v
( )
(czyli wyniku C-TFT) dla:
v
v
k f
T
k v
k
N
N
=
=
⋅
= ⋅
(14)
gdzie
v
t
N
=
1
∆
i oznacza odwrotność przedziału czasu, z którego pochodzą próbki czasowe, czyli
przedziału, w którym badamy funkcję czasu zakładając, Ŝe poza nim jest ona zerowa (a ściślej, Ŝe
próbki są zerowe poza tym przedziałem czasu).