Projektowanie filtrów FIR, © P.Korohoda
1
Projektowanie filtrów FIR
(o skończonej odpowiedzi impulsowej)
Przemysław Korohoda, KE, AGH
Zawartość instrukcji:
1 Wybrane zagadnienia z zakresu DSP
1.1 Filtr FIR o liniowej fazie
1.2 Projektowanie filtrów FIR - wstęp
1.3 Projektowanie filtrów FIR metodą okien czasowych
1.4 Projektowanie filtrów typu FIR metodami iteracyjnymi
1.4.1 Próbkowanie w dziedzinie częstotliwości
1.4.2. Błąd aproksymacji
1.4.3 Twierdzenie Remeza
1.4.4 Procedura iteracyjna Parksa-McClellana
1.5 Transformowanie dolnoprzepustowych filtrów FIR do innych postaci
1.6 Filtracja za pomocą filtrów FIR z wykorzystaniem FFT
2 Projektowanie filtrów cyfrowych typu FIR za pomocą pakietu MATLAB
2.2 Przykłady
2.2.1 Ilustracja faktu, iŜ zbyt krótkie widmo DFT moŜe prowadzić do charakterystyki
filtru odbiegającej od załoŜonej
2.2.2 Projektowanie górnoprzepustowego filtru FIR 10. rzędu, o zadanej charakterystyce
2.2.3 Porównanie wyniku projektowania filtru o zadanej charakterystyce amplitudowej z
wykorzystaniem procedur projektowych funkcji Matlab’a oraz „krok po kroku”
2.2.4 Filtracja filtrem typu FIR z wykorzystaniem FFT (splotu kołowego)
Projektowanie filtrów FIR, © P.Korohoda
2
Projektowanie filtrów FIR, © P.Korohoda
3
1 Wybrane zagadnienia z zakresu DSP
1.1 Filtr FIR o liniowej fazie
Jednym z argumentów przemawiających na korzyść filtrów typu FIR jest fakt, Ŝe filtry typu IIR mogą
posiadać fazę jedynie zbliŜoną do liniowej w pasmie przepustowym, a często faza ta dość znacznie
odbiega od fazy liniowej w pasmie przejściowym. Natomiast w przypadku filtrów typu FIR warunek
liniowej fazy w interesującym projektanta zakresie częstotliwości moŜna stosunkowo łatwo
zrealizować.
Filtr posiada liniową fazę, gdy jego charakterystyka częstotliwościowa moŜe być zapisana w
następującej postaci, zawierającej stałą
τ
:
(
)
H f
e
H f
j
f
( )
( )
=
⋅
− ⋅ ⋅ ⋅ ⋅
2
π τ
(1)
Jest to równowaŜne stwierdzeniu, Ŝe faza określona jest wzorem opisującym funkcję liniową od
zmiennej
f
, o nachyleniu
− ⋅ ⋅
2
π τ
:
(
)
ϕ
π
τ
H f
f
( )
= − ⋅ ⋅ ⋅
2
(2)
W systemach z czasem dyskretnym, chcąc uniknąć dodatkowych zniekształceń, naleŜy zadbać o to, by
przesunięcie w dziedzinie czasu wynikające z charakterystyki fazowej filtru wynosiło całkowitą
wielokrotność okresu próbkowania. Dlatego teŜ stała
τ
powinna przybierać tylko wartości całkowite.
Gdy będą to wartości dodatnie, to system będzie opoźniający, gdy ujemne, to odpowiadałoby to
systemowi dającemu odpowiedź z wyprzedzeniem. System (filtr) o zerowej fazie odpowiada wartości
τ
=
0
.
Warto przypomnieć, Ŝe pomiędzy częstotliwością cyfrową (
f
) stosowaną w przypadku D-TFT i
indeksem częstotliwościowym, typowym dla DFT, zachodzi następujący związek:
f
k
N
k
N
=
=
−
:
, ,
,
0 1
1
…
(3)
gdzie
N
oznacza długość ciągu poddawanego transformacji DFT.
Odpowiedź częstotliwościowa - czyli transformata D-TFT - filtru o odpowiedzi impulsowej
h n
[ ]
określona jest wzorem:
H f
h n e
j
f n
n
( )
[ ]
=
⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
+∞
∑
2
π
(4)
Dla skończonej, nieparzystej odpowiedzi impulsowej o długości
2
1
⋅
+
M
, symetrycznej względem
punktu
n
=
0
, wzór ten sprowadza się do postaci:
(
)
H f
h n e
h
h n
e
e
j
f n
n
M
M
j
f n
j
f n
n
M
( )
[ ]
[ ]
[ ]
=
⋅
=
+
⋅
+
− ⋅ ⋅ ⋅ ⋅
=−
− ⋅ ⋅ ⋅ ⋅
+ ⋅ ⋅ ⋅ ⋅
=
∑
∑
2
2
2
1
0
π
π
π
(5)
co ostatecznie daje:
(
)
H f
h
h n
f n
n
M
( )
[ ]
[ ]
cos
=
+
⋅ ⋅
⋅ ⋅ ⋅
=
∑
0
2
2
1
π
(6)
Projektowanie filtrów FIR, © P.Korohoda
4
Otrzymana odpowiedź częstotliwościowa nie posiada części urojonej. Oznacza to, Ŝe faza, po
sprowadzeniu do przedziału
(
,
]
− +
π π
, moŜe przybierać jedynie wartości 0 lub
π
. Wartość fazy
π
oznacza, Ŝe po filtracji nastąpi zmiania znaku składowej sygnału o odpowiedniej częstotliwości. W
przypadku określonej odpowiedzi impulsowej pojawienie się ujemnych wartości
H f
( )
zaleŜy od
wartości
f
. Metody projektowania filtrów dopuszczają pewien ograniczony błąd w charakterystyce
częstotliwościowej w stosunku do filtru idealnego. Nie spotyka się jednak załoŜeń projektowych
dopuszczających, by w pasmie przepustowym i istotnej części pasma przejściowego amplituda
charakterystyki filtru zbliŜała się do wartości zero. Wyklucza to zatem moŜliwość zmiany znaku części
rzeczywistej charakterystyki filtru we wspomnianym zakresie częstotliwości, a więc faza w tym zakresie
nie moŜe zawierać skoków o wartości
±
π
. MoŜliwe jest to tylko w zakresie niewielkich amplitud,
gdzie nie odgrywa to juŜ większego znaczenia. Stąd bierze się stwierdzenie, Ŝe praktycznie filtr taki
posiada zerową fazę, choć ze ścisłego punktu widzenia w pewnym zakresie częstotliwości faza ta moŜe
mieć wartość
π
. Filtr o zerowej - w myśl powyŜszej interpretacji - fazie powinien spełniać warunek
symetrii względem punktu
n
=
0
. Oznacza to jednak, Ŝe odpowiedź impulsowa jest nieparzystej
długości oraz, Ŝe filtr nie jest przyczynowy.
W ramach ćwiczeń rozwaŜane bedą wyłącznie filtry FIR o nieparzystej długości zatem dyskusja
problematyki filtrów o parzystej długości będzie tutaj pominięta. Odpowiednie dla takiego przypadku
rozwaŜania moŜna znaleźć w literaturze uzupełniającej i przestudiowanie tego zagadnienia moŜna
potraktować jako samodzielne ćwiczenie dodatkowe.
Nieprzyczynowy filtr FIR o zerowej fazie moŜna łatwo przekształcić do postaci filtru przyczynowego
FIR o liniowej fazie przez opóźnienie odpowiedzi impulsowej o
M
taktów. JeŜeli przez
H
f
0
( )
oznaczona zostanie odpowiedź częstotliwościowa filtru nieprzyczynowego o zerowej fazie, to
otrzymany filtr będzie miał następującą charakterystykę częstotliwościową:
H f
e
H
f
j
f M
( )
( )
=
⋅
− ⋅ ⋅ ⋅ ⋅
2
0
π
(7)
1.2 Projektowanie filtrów FIR
Projektowanie filtru typu FIR, czyli o skończonej odpowiedzi impulsowej, polega najczęściej na
wyznaczeniu ciągu odpowiedzi impulsowej, który przez operację splatania liniowego z sygnałem
wejściowym da poŜądany rezultat. Zazwyczaj przez określenie “poŜądany rezultat” rozumie się
odpowiednią zmianę w charakterystyce częstotliwościowej sygnału, przy czym powinna to być
charakterystyka w rozumieniu D-TFT (splot liniowy).
Projektowanie filtrów FIR, moŜe wydawać się zadaniem prostym. Wystarczyłoby w dziedzinie
częstotliwości zdefiniować poŜądaną charakterystykę, wyznaczyć odwrotną transformatę Fouriera i w
ten sposób otrzymać odpowiedź impulsową realizującą projektowany filtr. Rozwiązanie takie nie jest
jednak zwykle stosowane z dwóch zasadniczych, powiązanych wzajemnie, powodów:
1) filtr realizowany za pomocą odpowiedzi impulsowej splatanej z nadchodzącym na bieŜąco ciągiem
sygnału nie realizuje splotu kołowego, lecz splot liniowy, natomiast powyŜsza propozycja wiąŜe się z
wykorzystaniem odwrotnej wersji FFT, czyli byłaby odpowiednia dla filtracji właśnie poprzez splot
kołowy;
2) w pewnych przypadkach moŜna by uznać, Ŝe widmo filtru określone w dziedzinie DFT jest
dostatecznie dobrym przybliŜeniem widma w dziedzinie D-TFT - widmo DFT moŜna otrzymać przez
spróbkowanie jednego okresu widma D-TFT - jednak, by tak moŜna było przyjąć, to ilość próbek , czyli
długość ciągu transformaty DFT, a zatem i odpowiedzi impulsowej, musiałaby być zbyt długa dla
większości praktycznych zastosowań.
Najczęściej stosowane rozwiązania moŜna podzielić na dwie grupy:
Projektowanie filtrów FIR, © P.Korohoda
5
a) przez podejście identyczne do opisanego powyŜej dla dostatecznie długiego ciągu transformaty DFT
- by moŜna ją było uznać, za dobre przybliŜenie transformaty D-TFT - i następnie skracanie
wyznaczonej odpowiedzi impulsowej w taki sposób, by w wyniku tego zabiegu widmo filtru ulegało jak
najmniej szkodliwym zmianom;
b) przez wykorzystanie teorii aproksymacji, dzięki czemu zdefiniowanie widma tylko w wybranych
punktach częstotliwości moŜe prowadzić do spełnienia odpowiednich załoŜeń w kaŜdym punkcie
częstotliwości.
Oba rozwiązania będą dokładniej przedstawione w kolejnych rozdziałach.
1.3 Projektowanie filtrów FIR metodą okien czasowych
Jak juŜ wspomniano, długą odpowiedź impulsową wyznaczoną z odwrotnej transformacji FFT
zadanego widma moŜna skracać, starając się przy tym zachować najistotniejsze cechy widma. Do tego
celu słuŜą okna czasowe, czyli ciągi, które są mnoŜone - element po elemencie - przez daną odpowiedź
impulsową. Warto wspomnieć, Ŝe wprowadzanie zmian w odpowiedzi impulsowej owocuje efektem
Gibbsa w dziedzinie częstotliwości.
Najprostszym oknem czasowym jest okno prostokątne, obcinające wprost odpowiedź do długości
2
1
⋅
+
M
:
[
]
w
n
dla
n
M
M
dla
pozost n
Π
[ ]
,
.
=
∈ −
+
1
0
(8)
Przy podawaniu indeksów w definicji okien przyjęto załoŜenie, Ŝe rozwaŜana jest odpowiedź
impulsowa o zerowej fazie. Przejście do odpowiedzi o fazie liniowej będzie się wiązało z
przesunięciem indeksów tak, by po skróceniu odpowiedź impulsowa rozpoczynała się w punkcie
n
=
0
.
Zatem odpowiedź impulsowa po zastosowaniu okna (8), będzie określona następująco:
[
]
h n
h n w
n
h n
dla
n
M
M
dla
pozost n
Π
Π
[ ]
[ ]
[ ]
[ ]
,
.
=
⋅
=
∈ −
+
0
(9)
MnoŜenie w dziedzinie indeksów czasowych odpowiada splataniu w dziedzinie częstotliwości. Zatem
chcąc przewidzieć, jakie zajdą zmiany w widmie filtru, wystarczy zbadać widmo zastosowanego okna i
zastanowić się jak splatanie takiego widma z widmem filtru wpłynie na końcowy rezultat. Z tego punktu
widzenia najkorzystniejszym widmem okna byłaby oczywiście delta Kroneckera z odpowiednią
amplitudą, określona w dziedzinie częstotliwości DFT, co odpowiadałoby delcie Diraca w przypadku
ciągłej dziedziny częstotliwości D-TFT. Wynika to wprost z faktu, Ŝe splatanie z ciągiem delty
Kroneckera nie wnosi Ŝadnych zmian do splatanego ciągu - w tym przypadku jest to ciąg próbek widma
filtru. Jak łatwo sprawdzić, widmo okna (8) dość znacznie róŜni się od delty Kroneckera i dlatego teŜ
okno to nie jest często stosowane. Na rysunku 1 pokazano charakterystykę amplitudową okna o
długości 17, wyliczoną przez FFT dla ciągu wydłuŜonego zerami do długości 1024. Charakterystkę
przedstawiono po przesunięcie punktu ‘dc’ do środka zakresu. Wykres ma kształt zafalowań - kaŜde
takie zafalowanie nazywa się “listkiem”. Listek obejmujący punkt o częstotliwości 0 nazywany jest
listkiem głównym, pozostałe to listki boczne.
W celu stwierdzenia, na ile dane okno spełnia typowe wymagania definiuje się dwa parametry
liczbowe: szerokość listka głównego i tłumienie dla listków bocznych. Szerokość listka głównego to
podwojona wartość częstotliwości cyfrowej, dla której część rzeczywista widma okna przechodzi
oddalając się od punktu
f
=
0
pierwszy raz przez 0 (metoda taka jest poprawna tylko gdy część
urojona widma okna jest zerowa!). Tłumienie dla listków bocznych określa się w dB jako odstęp
pomiędzy maksymalną wartością amplitudy w listkach bocznych i wartością amplitudy dla
f
=
0
.
Punkt
f
=
0
określa zatem wartość odniesienia, co oznacza, Ŝe amplituda wynosi w nim 0dB.
Projektowanie filtrów FIR, © P.Korohoda
6
Rys.1 Część rzeczywista i amplituda widma dla okna prostokątnego o długości 17 elementów (środek
okna w punkcie n=0)
Dla przykładu przedstawionego na rys.1 szerokość listka głównego wynosi około 0.12, natomiast
tłumienie listków bocznych 0.66dB. Pierwszy parametr jest bezwymiarowy, poniewaŜ wyraŜony jest w
częstotliwościach cyfrowych
f
. W przypadku określenia częstotliwości próbkowania w Hz, moŜna ten
parametr wyrazić równieŜ w Hz opierając się na tym, iŜ wartość
f
=
1
odpowiada częstotliwości
próbkowania.
Okno jest tym lepsze im mniejsza jest szerokość listka głównego i większe jest tłumienie dla listków
bocznych - wtedy bowiem widmo jest bardziej zbliŜone do delty Kroneckera z odpowiednią amplitudą
(lub delty Diraca, intepretując widmo jako charakterystykę w dziedzinie D-TFT). Jednak szerokość
listka głównego jest odwrotnie proporcjonalna do długości okna w dziedzinie czasu. PoniewaŜ istotą
stosowania okien jest skracanie odpowiedzi impulsowych, więc prowadzi to do sprzecznych wymagań i
dlatego teŜ istnieje szereg róŜnych propozycji okien czasowych, z których najbardziej znane to:
1) okno Bartletta,
2) okno Blackmanna,
3) okno Hamminga,
4) okno von Hanna, zwane często oknem Hanninga,
5) okno Kaisera (najbardziej uniwersalne),
6) okno Czebyszewa.
Wzory opisujące powyŜsze okna w dziedzinie czasu moŜna znaleźć w literaturze uzupełniającej.
Po wyliczeniu ciągu reprezentującego dane okno, ciąg dalszy jest taki sam jak opisano dla okna
prostokątnego. Za kaŜdym razem chcąc skrócić odpowiedź impulsową do określonej długości naleŜy
pełną odpowiedź impulsową pomnoŜyć przez ciąg okna dla okna o szerokości identycznej z poŜądaną
długością odpowiedzi impulsowej. NaleŜy jedynie pamiętać o zsynchronizowaniu obu mnoŜonych
ciągów z punktu widzenia indeksów czasowych - patrz teŜ uwaga na początku tego rozdziału.
1.4 Projektowanie filtrów typu FIR metodami iteracyjnymi
W rozdziale tym zostanie przedstawiona metoda naleŜąca do grupy metod projektowania filtrów o
charakterystyce równomiernie falistej. Oznacza to, Ŝe projektując filtr dopuszcza się istnienie zafalowań
charakterystyki częstotliwościowej filtru w pasmie przepustowym i zaporowym o załoŜonej, stałej
amplitudzie. Zafalowania te oraz pasmo przejściowe o niezerowej szerokości stanowią róŜnicę takiego
filtru w stosunku do filtru idealnego. Charakterystyka fazowa zaprojektowanego filtru jest liniowa w
znaczeniu opisanym w rozdziale 1.1 tej instrukcji.
Projektowanie filtrów FIR, © P.Korohoda
7
1.4.1 Próbkowanie w dziedzinie częstotliwości
Charakterystyka częstotliwościowa filtru cyfrowego przeznaczonego do filtracji realizującej splot
liniowy powinna być określana jako transformata D-TFT. Często jednak poprzestaje się na
transformacie DFT zawierającej próbki transformaty D-TFT rozłoŜone równomiernie na osi
częstotliwości. Okazuje się jednak, Ŝe równomierne rozmieszczenie próbek znacznie utrudnia właściwe
zdefiniowanie charakterystyki do załoŜeń projektowych. Przykładowo, znaczną poprawę w tłumieniu
dla pasma zaporowego uzyskuje się przez przyjęcie, iŜ jedna lub więcej próbek znajduje się w zakresie
pasma przejściowego. W przypadku równomiernego rozmieszczenia próbek oznacza to znaczne
poszerzenie tego zakresu, co nie jest efektem korzystnym. Dlatego teŜ opracowano szereg metod
umoŜliwiających definiowanie charakterystyki częstotliwościowej filtru w sposób inny niŜ przez
podanie wartości widma rozmieszczonych równomiernie na osi częstotliwości. Część z nich zakłada
wartości częstotliwości ograniczające pasmo przepustowe i zaporowe oraz rząd filtru FIR i wyznacza
współczynniki filtru oraz moŜliwą do osiągnięcia dokładność aproksymacji charakterystyki, natomiast
część przyjmuje jako załoŜenia rząd filtru oraz dopuszczalny błąd aproksymacji, wyliczając oprócz
współczynników graniczne wartości częstotliwości.
1.4.2. Błąd aproksymacji
Niech
H
f
Z
( )
oznacza zadaną charakterystykę częstotliwościową - transformatę D-TFT - filtru
idealnego, natomiast
H f
( )
aproksymację tej idealnej charakterystyki za pomocą realizowalnego filtru
cyfrowego typu FIR. W celu liczbowego określenia na ile otrzymana charakterystyka róŜni się od
zadanej wprowadza się szereg wskaźników. Wskaźniki te zawierają funkcję róŜnicy charakterystyk
zadanej i otrzymanej oraz funkcję wagową umoŜliwiającą skupienie uwagi na wybranych fragmentach
widma. Tam, gdzie funkcja wagowa jest większa, tam minimalizacja danego wskaźnika da lepszą
zgodność obu charakterystyk. Jednym z najprostszych sposobów wykorzystania funkcji wagowej jest
zadanie większej precyzji apkroksymacji w zakresie pasma przepustowego i mniejszej w zakresie
pazma zaporowego. Jednym z bardziej popularnych wskaźników - czyli kryteriów jakości aproksymacji
- jest błąd średniokwadratowy wyznaczony dla ciągłego zakresu częstotliwości:
J
W f
H f
H
f
df
Z
1
2
0
0 5
=
⋅
−
∫
( )
( )
( )
.
(10)
Przedział całkowania we wzorze (10) obejmuje połowę jednego okresu transformaty D-TFT. Wynika to
z załoŜenia, iŜ odpowiedź impulsowa filtru jest rzeczywista.
W niektórych przypadkach rozwaŜa się jedynie aproksymację w wybranych punktach charakterystyki
częstotliwościowej - na przykład w punktach odpowiadających transformacie DFT:
dla parzystej długości ciągu
N
:
J
N
W f
H f
H
f
n
n
Z
n
n
N
2
2
0
2
1
2
1
=
+
⋅
⋅
−
=
∑
(
)
(
)
(
)
(11a)
dla nieparzystej długości ciągu
N
:
J
N
W f
H f
H
f
n
n
Z
n
n
N
2
2
0
1
2
1
1
2
=
− ⋅
⋅
−
=
−
∑
(
)
(
)
(
)
(11b)
Wskaźnik jakości (11) nie uwzględnia w ogóle wartości charakterystyk poza wskazanymi punktami,
dlatego teŜ waŜne jest odpowiednie wybranie tych punktów na osi częstotliwości. Kryterium (11) nosi
nazwę najmniejszych kwadratów i moŜe być stosowane równieŜ przy nierównomiernym rozłoŜeniu
N / 2
1
+
punktów
f
n
(według wzoru (11a)). W takiej sytuacji (11a) nie ma juŜ bezpośredniego
związku z transformatą DFT, a jedynie z nierównomiernie próbkowaną transformatą D-TFT.
Inne bardziej popularne kryterium bazuje na bardzo prostej funkcji błędu zawierającej funkcję wag:
Projektowanie filtrów FIR, © P.Korohoda
8
[
]
E f
W f
H f
H
f
Z
( )
( )
( )
( )
=
⋅
−
(12)
Wskaźnik jakości, zwany teŜ kryterium Czebyszewa, określa się wówczas następująco:
(
)
J
E f
f
3
0 0 5
=
∈<
>
max
( )
, .
(13)
Wzór (13) oznacza maksymalną waŜoną rozbieŜność charakterystyki zadanej i aproksymującej w
połowie okresu charakterystyki częstotliwościowej filtru.
1.4.3 Twierdzenie Remeza
Podstawą dla prezentowanej tutaj metody jest twierdzenie znane między innymi pod nazwami:
twierdzenie Remeza lub twierdzenie o przerzutach.
Twierdzenie Remeza:
Niech F będzie dowolnym domkniętym podprzedziałem przedziału częstotliwości cyfrowych
0
0 5
≤ ≤
f
,
(przedział F obejmuje pasmo przepustowe i zaporowe). JeŜeli funkcja błędu
E f
( )
zdefiniowana według (12) wykazuje w tym przedziale F co najmniej
M
+
2
przerzuty, czyli:
(
)
E f
E f
E
E f
k
k
(
)
(
)
max
( )
= −
= ±
=
−
1
dla
f
f
f
f
F
M
k
0
1
1
≤
≤ ≤
∧
∈
+
…
,
to istnieje dokładnie jeden filtr rzędu
M
typu FIR posiadający charakterystykę częstotliwościową,
będącą optymalną - w sensie minimalnej wartości (13) - aproksymacją załoŜonej charakterystyki filtru,
przy czym powyŜszy warunek istnienia takiego filtru FIR jest konieczny i wystarczający.
Warto zauwaŜyć, Ŝe punktów charakterystyki spełniających warunek twierdzenia moŜe być więcej niŜ
M
+
2
. W szczególności dotyczy to końców przedziału
<
>
0 0 5
, .
oraz punktów brzegowych pasma
przejściowego. MoŜna teŜ przypomnieć, Ŝe rząd filtru
M
oznacza
M
+
1
współczynników filtru,
czyli odpowiedź impulsową o długości
M
+
1
.
1.4.4 Procedura iteracyjna Parksa-McClellana
Opisana poniŜej procedura bywa równieŜ nazywana algorytmem Remeza. Metoda polega na
minimalizowaniu wskaźnika
J
3
przy spełnieniu warunków twierdzenia Remeza. oznacza to, Ŝe w
punktach ekstremalnych funkcji błędu (12), spełniających twierdzenie Remeza, zaleŜność (12) moŜna
przepisać do postaci:
[
]
W f
H f
H
f
J
k
M
k
k
Z
k
k
(
)
(
)
(
)
(
)
:
, ,
,
⋅
−
= −
⋅
=
+
+
1
0 1
1
1
3
…
(14)
Podstawienie w (14) za
H f
k
(
)
na podstawie równania (6) dla wszystkich wartości k, daje układ
równań, który moŜna zapisać macierzowo:
(
)
(
)
( )
(
)
(
)
( )
(
)
(
)
( )
(
)
(
) ( )
1
2
2
1
1
2
2
1
1
2
2
1
1
2
2
1
0
2
1
2
0
0
0
1
1
1
1
1
1
1
cos
cos
cos
cos
cos
cos
(
)
cos
cos
(
)
[ ]
[ ]
⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅
⋅ ⋅ ⋅
−
⋅ ⋅
⋅ ⋅
⋅
−
⋅ ⋅
⋅ ⋅
⋅
−
⋅
⋅
⋅
+
+
+
+
π
π
π
π
π
π
π
π
f
f
M
W f
f
f
M
W f
f
f
M
W f
f
f
M
W f
h
h
M
M
M
M
M
M
M
M
⋯
⋯
⋮
⋮
⋱
⋮
⋮
⋯
⋯
⋮
⋮
⋮
h M
J
H
f
H
f
H
f
H
f
Z
Z
Z
M
Z
M
[
]
(
)
(
)
(
)
(
)
3
0
1
1
=
+
⋮
(15)
Projektowanie filtrów FIR, © P.Korohoda
9
Warto tu przypomnieć, Ŝe jest jedynie kwestią interpretacji, czy dany filtr potraktujemy jako
nieprzyczynowy o zerowej fazie, czy teŜ przyczynowy - opoźniający - o fazie liniowej. We wzorze (15)
przyjęto pierwszą intepretację - róŜnica byłaby widoczna w indeksach czasowych współczynników
odpowiedzi impulsowej.
Rozwiązanie układu (15) prowadzi przy załoŜonych punktach
f
k
, do wyznaczenia odpowiedzi
impulsowej i wartości błędu
J
3
. Zgodnie z twierdzeniem Remeza, jeŜeli punkty
f
k
leŜą w lokalnych
estremach funkcji aproksymującej, to wartość błędu jest mininalna. Jednak na początku poszukiwań
filtru nie są zwykle znane odpowiednie punktu na osi częstotliwości. Dlatego teŜ stosuje się następującą
procedurę iteracyjną, której dowód zbieŜności moŜna znaleźć w literaturze:
0) równanie (15) jest rozwiązywane dla wstępnego układu
M
+
2
punktów
f
k
,
1) dla otrzymanego filtru wyznaczana jest charakterystyka częstotliwościowa, której lokalne ekstrema
określają kolejny układ
M
+
2
punktów
f
k
,
2) ponownie rozwiązywany jest układ równań (15), z którego otrzymuje się współczynniki filtru i
wartość kryterium jakości,
3) jeŜeli w wyniku kolejnej iteracji wartość kryterium jakości nie uległa zmianie, to procedura jest
zakończona, w przeciwnym przypadku naleŜy powrócić do punktu 1).
W praktyce ekstrema nowej charakterystyki częstotliwościowej z kroku 1) wyznacza się przez
zastosowanie interpolacyjnego wzoru Lagrange’a.
1.5 Transformowanie dolnoprzepustowych filtrów FIR do innych postaci
Projektowanie filtru FIR polega na wyznaczeniu odpowiedzi impulsowej odpowiadającej załoŜonej
charakterystyce częstotliwościowej. Wiele z metod projektowania prowadzi do uzyskania odpowiedzi
impulsowej filtru dolnoprzepustowego. Dość waŜnym zatem zagadnieniem jest: jak korzystając z tych
odpowiedzi impulsowych wyznaczyć odpowiedź impulsową filtru innego typu niŜ dolnoprzepustowy?
Wśród moŜliwych rozwiązań moŜna wyróŜnić trzy główne podejścia:
1) przez skorzystanie z właściwości filtrów kwadraturowych - umoŜliwia to zamianę filtru
dolnoprzepustowego na górnoprzepustowy;
2) przez zastosowanie modulacji, czyli przesunięcia widma filtru w dziedzinie częstotliwości, co
otrzymuje się mnoŜąc odpowiedź impulsową fitru dolnoprzepustowego przez ciąg pochodzący z
próbkowanej funkcji kosinus;
3) poprzez bezpośrednie wykorzystanie liniowości transformacji Fouriera i złoŜenie poszukiwanego
widma z widm filtrów, które moŜna łatwo otrzymać - na przykład z widm filtrów
dolnoprzepustowch.
Podejście 1) polega na odwróceniu kolejności elementów odpowiedzi impulsowej oraz zmianie znaku
co drugiego elementu.
Podejście 2) opiera się na pomnoŜeniu ciągu
h n
[ ]
odpowiedzi impulsowej filtru dolnoprzepustowego
przez odpowiedni ciąg, przykładowo przesunięcie środka pasma filtru do punktu “k” naleŜącego do
przedziału od 0 do
N
2
1
−
otrzyma się po następującej modyfikacji odpowiedzi impulsowej:
h n
n k
N
h n
1
2
2
[ ]
cos
[ ]
= ⋅
⋅ ⋅ ⋅
⋅
π
(16)
Warto pamiętać, Ŝe, w przypadku szerokiego pasma przepustowego filtru FDP i niewielkiej wartości
“k”, przesunięte w dziedzinie DFT widma filtru FDP nałoŜą się na siebie, co będzie wynikiem splotu
kołowego w dziedzinie częstotliwości.
Podejście 3) opiera się na liniowości transformacji Fouriera. JeŜeli poŜądane widmo moŜna otrzymać z
kombinacji liniowej widm posiadanych filtrów, to odpowiedź impulsowa szukanego filtru będzie
Projektowanie filtrów FIR, © P.Korohoda
10
wynikiem takiej samej kombinacji liniowej odpowiedzi impulsowych tych samych filtrów.
Szczególnym przypadkiem moŜe być na przykład filtr pasmowoprzepustowy, którego charakterystykę
częstotliwościową moŜna otrzymać przez odjęcie charakterystyk częstotliwościowych dwóch filtrów
dolnoprzepustowych o róŜnych częstotliwościach granicznych. Odpowiedź impulsową szukanego filtru
otrzyma się przez odjęcie odpowiedzi impulsowych wspomnianych filtrów dolnoprzepustowych. Innym
szczególnym
przypadkiem
moŜe
być
odjęcie
charakterystyki
częstotliwościowej
filtru
dolnoprzepustowego od widma składającego się z samych wartości “1”. W ten sposób moŜna otrzymać
na przykład filtr górnoprzepustowy. Oznacza to, Ŝe od odpowiedzi impulsowej filtru
wszechprzepustowego o zerowej fazie naleŜy odjąć odpowiedź impulsową filtru dolnoprzepustowego.
Odpowiedź impulsowa filtru wszechprzepustowego, to
h
n
d n
ALL
[ ]
[ ]
=
, czyli delta Kroneckera.
Zatem odpowiedź impulsowa filtru górnoprzepustowego (FGP) komplementarnego względem “1” do
danego filtru dolnoprzepustowego (FDP), to:
h
n
d n
h
n
FGP
FDP
[ ]
[ ]
[ ]
=
−
(17)
Warto zwrócić uwagę, Ŝe dwa filtry kwadraturowe i oba filtry powiązane zaleŜnością (17) są
komplementarne w myśl innych definicji komplementarności.
1.6 Filtracja za pomocą filtrów FIR z wykorzystaniem FFT
W dotychczasowych rozwaŜaniach przyjmowano, Ŝe przez filtr FIR rozumie się jego odpowiedź
impulsową, bezpośrednio stosowaną do wyznaczania splotu liniowego. Charakterystyka widmowa filtru
stanowiła jedynie wstępne załoŜenie projektowe. Warto przypomnieć, Ŝe dla dostatecznie długiej
odpowiedzi impulsowej znacznie szybsze obliczeniowo moŜe się okazać zastosowaniu algorytmu
overlap-save (lub overal-add). Podejście to polega na realizacji splotu liniowego poprzez filtrację w
dziedzinie DFT. Wspomniane algorytmy były wprawdzie przedmiotem jednego z poprzednich ćwiczeń,
jednak dla przypomnienia na rys.2 i rys.3 przedstawiono istotę obu metod. W obu metodach ciąg
wejściowy dzielony jest w miarę napływania danych na bieŜąco na bloki o długości
L
x
. KaŜdy blok, w
chwili jego filtracji, jest uzupełniony zerami (algorytm overlap-add) lub początkiem kolejnego bloku
(algorytm overlap-save). Długość tego uzupełnienia jest równa długości odpowiedzi impulsowej filtru
FIR pomniejszonej o 1, czyli
L
h
−
1
. W kolejnym kroku blok o łącznej długości
L
L
x
h
+
−
1
poddawany jest FFT, po czym przeprowadzana jest filtracja za pomocą mnoŜenia przez transformatę
DFT odpowiedzi impulsowej filtru i wynik filtracji poddawany jest odwrotnej transformacji - czyli
IFFT. Filtracja taka odpowiada splotowi kołowemu bloku danych o długości
L
L
x
h
+
−
1
. W celu
sprowadzenia wyniku do postaci odpowiadającej fragmentowi splotu liniowego o długości
L
x
stosuje
się następujące zabiegi:
a) w algorytmie overlap-add:
Wynik splotu kołowego, o długości
L
L
x
h
+
−
1
, dzielony jest na dwie części - pierwszą, o
długości
L
x
i drugą, o długości
L
h
−
1
. Do kolejnych początkowych
L
h
−
1
elementów pierwszej części dodawane są kolejne elementy zapamiętane w buforze,
znajdujące się tam w wyniku przetwarzania poprzedniego bloku. Z kolei druga część -
właśnie o długości
L
h
−
1
- jest wpisywana jaka nowa zawartość bufora, do
wykorzystania przy filtracji kolejnego bloku. Cała, tak zmodyfikowana, pierwsza część o
długości
L
x
stanowi odpowiedni fragment splotu liniowego.
b) w algorytmie overlap-save:
Wynik splotu kołowego, o długości
L
L
x
h
+
−
1
, jest dzielony na dwie części. W tym
przypadku pierwsza część ma długość
L
h
−
1
i jest po prostu odrzucana, natomiast jako
odpowiedni fragment splotu liniowego pozostawia się część drugą, o długości
L
x
.
Projektowanie filtrów FIR, © P.Korohoda
11
Rys.2. Algorytm overlap-add
Rys.3. Algorytm overlap-save
W ramach ćwiczenia samodzielnego naleŜy się zastanowić, w jaki sposób, bazując na liniowości splotu
oraz DFT, moŜna uzasadnić poprawność obu algorytmów.
Projektowanie filtrów FIR, © P.Korohoda
12
2 Projektowanie filtrów cyfrowych typu FIR za pomocą pakietu
MATLAB
Wybrane funkcje pakietu Matlab - część 1
funkcje generowania okien
Nazwa funkcji
Opis funkcji
boxcar(N)
N-elementowy ciąg samych jedynek;
triang(N)
N-elementowy ciąg w kształcie trójkąta bez końcowych zer;
bartlett(N)
N-elementowy ciąg teŜ w kształcie trójkąta, ale z końcowymi zerami;
blackman(N, opcja )
N-elementowy ciąg w kształcie dzwonu, schodzący do zera; opcja
‘symmetric’ lub ‘periodic’, bez opcji domyślnie - ‘symmetric’;
hamming(N, opcja )
N-elementowy ciąg w kształcie dzwonu, nie schodzący do zera; opcja
‘symmetric’ lub ‘periodic’, bez opcji domyślnie - ‘symmetric’;
hanning(N, opcja )
N-elementowy ciąg w kształcie dzwonu, schodzący prawie do zera; opcja
‘symmetric’ lub ‘periodic’, bez opcji domyślnie - ‘symmetric’;
kaiser(N, Beta)
N-elementowy ciąg w kształcie dzwonu; drugi, opcjonalny parametr - beta
- od 0 do 700, powoduje, Ŝe okno zmienia się od ciągu samych jedynek,
po ciąg bliski delcie Kroneckera;
chebwin(N, R)
N-elementowy ciąg rozpoczynający się i kończący na elementach
równych “1”, pomiędzy którymi znajduje się ciąg mniejszych wartości o
maksimum w środku ciągu; drugi parametr, R, dotyczy odpowiedzi
częstotliwościowej (D-TFT) i określa w dB tłumienie listków bocznych;
Uwaga - funkcje generowania okien zwracają wektory kolumnowe.
W przypadku
Wybrane funkcje pakietu Matlab - część 2
funkcje realizujące procedury projektowania filtrów FIR
Nazwa funkcji
Opis funkcji
fir1
wyznaczenie metodą okien czasowych współczynników filtru FIR o liniowej
fazie spełniającego wymagania określone w postaci granic pasm
przepustowych i zaporowych,
fir2
wyznaczenie metodą okien czasowych współczynników filtru FIR o liniowej
fazie spełniającego wymagania określone w postaci wartości charakterystyki
we wskazanych punktach częstotliwościowych.
firls
wyznaczenie współczynników filtru FIR o liniowej fazie spełniającego
wymagania określone w postaci wartości charakterystyki we wskazanych
punktach
częstotliwościowych
przez
optymalizację
charakterystyki
projektowanego filtru w sensie najmniejszych kwadratów
remez
wyznaczenie współczynników filtru FIR o liniowej fazie i charakterystyce
równomiernie falistej, spełniającego wymagania określone w postaci wartości
charakterystyki
we
wskazanych
punktach
częstotliwościowych przez
minimalizację odchylenia maksymalnego (13)
Projektowanie filtrów FIR, © P.Korohoda
13
Sposoby wykorzystania funkcji realizujących procedury projektowania filtrów FIR -
część 1
sposoby wywołania funkcji
Instrukcja
Opis
h=fir1(n, Fn,
opcje_pasmowe,
opcjonalne_okno)
Funkcja zwraca wartość współczynników h dla filtru dolnoprzepustowego o
podanym rzędzie n oraz górnej 3-decybelowej znormalizowanej
częstotliwości Fn. Podając odpowiednie opcje moŜna projektować filtry
górno- i pasmowoprzepustowe - patrz kolejna tabela. W roli ostatniego,
opcjonalnego, parametru moŜna zastosować wywołanie funkcji generującej
wybrane okno czasowe, np. bartlett(n+1); w przypadku pominięcia tego
parametru stosowane jest okno Hamminga.
h=fir2(n, F, M)
F to wektor częstotliwości znormalizowanych o początkowej wartości 0 i
końcowej 1. Wartość 1 odpowiada połowie częstotliwości próbkowania.
Kolejne wartości częstotliwości muszą być ułoŜone w parządku
niemalejącym. M to wektor amplitud odpowiadających kolejnym
wartościom częstotliwości F i musi być tej samej długości co wektor F.
h=firls(n, F, M, W)
h=remez(n, F, M, W)
F to wektor częstotliwości znormalizowanych podawany jako ciąg par liczb,
ograniczających pasma częstotliwości zawarte w przedziale częstotliwości
znormalizowanych: 0
÷
1. Wartości wektora M są traktowane jako wartości
amplitudy na krańcach tych pasm częstotliwości. Jedna para wartości z
wektora F odpowiada jednej parze wartości z wektora M. Dlatego długość
wektora M musi być taka sama jak F. Zatem zadana charakterystyka filtru
jest określona jako ciąg par punktów {F(k), M(k)} i {F(k+1), M(k+1)}.
Charakterystykę zadaną otrzymuje się łącząc te pary odcinkiem prostej, lecz
tylko dla nieparzystych wartości k. Dla k parzystych odcinek łączący
{F(k), M(k)} i {F(k+1), M(k+1)} traktuje się jako niezdefiniowany zatem
nie jest on uwzględniany w procesie optymalizacji w trakcie
projektowania filtru. Ostatni, opcjonalny, parametr to wektor wag W. JeŜeli
nie zostanie podany, to przyjęty będzie wektor jedynkowy. Wektor ten
powinien mieć długość dwukrotnie mniejszą niŜ wektory F i M, gdyŜ
dotyczy
tylko
pasm
częstotliwości
uwzględnianych
w
procesie
optymalizacji - zatem “co drugiego” przedziału wynikającego z podziału
odcinka częstotliwości znormalizowanej 0-1. Procedury realizowane przez
funkcje filrls oraz remez róŜnią się optymalizowanym kryterium jakości.
Sposoby wykorzystania funkcji do projektowania filtrów FIR - część 2
wybrane opcje pasmowe funkcji fir1
Opcja
Opis
‘high’
Wymusza projektowanie filtru górnoprzepustowego.
‘stop’
Wymusza projektowanie filtru pasmowozaporowego. Wtedy Fn powinien być wektorem
dwuelementowym określającym graniczne, znormalizowane częstotliwości pasma
zaporowego. Podając Fn jako wektor dwuelementowy bez opcji ‘stop’, zostanie
wymuszone zaprojektowanie filtru pasmowoprzepustowego. Fn będzie wówczas parą
częstotliwości, ograniczającą pasmo przepustowe.
S
zczegóły dotyczące poleceń proszę odczytać przy pomocy polecenia help.
Pojęcie znormalizowanej częstotliwości oznacza w przypadku powyŜszych informacji wartość
częstotliwości cyfrowej - lub pulsacji - przeskalowanej tak, Ŝe połowa częstotliwości - lub pulsacji -
próbkowania odpowiada wartości “1”. Po takiej normalizacji nie ma Ŝadnego znaczenia, czy wartości
znormalizowane nazwiemy znormalizowanymi częstotliwościami cyfrowymi, czy teŜ znormalizowanymi
pulsacjami cyfrowymi.
NaleŜy podkreślić, Ŝe powyŜszy wykaz nie objemuje wszystkich funkcji pakietu Matlab
przeznaczonych do projektowania filtrów typu FIR, ani teŜ wszystkich wariantów uŜycia wymienionych
funkcji.
Projektowanie filtrów FIR, © P.Korohoda
14
2.2 Przykłady
2.2.1 Ilustracja faktu, iŜ zbyt krótkie widmo DFT moŜe prowadzić do charakterystyki filtru
odbiegającej od załoŜonej
Najpierw naleŜy wygenerować w dziedzinie DFT widmo idealnego filtru dolnoprzepustowego o
zerowej fazie:
>> H127=zeros(1,127);
>> H127(1:7)=1;
>> H127(127:-1:127-5)=1;
>> h127=real(ifft(H127)); juŜ bez sprawdzania poprawności symetrii w dziedzinie DFT
generowanie widma w dziedzinie DFT (poprzez FFT) przez przedłuŜenie odpowiedzi impulsowej
zerami (zera wstawiane są do części środkowej ciągu przeznaczonego do DFT, dlatego Ŝe wyjściowy
filtr miał zerową fazę, zatem był nieprzyczynowy):
>> h1016=[h(1:64), zeros(1,1016-127),h(65:127)]; poniewaŜ 1016=8*127;
>>H1016=fft(h1016);
i weryfikacja tezy, Ŝe widmo amplitudowe po przedłuŜeniu odpowiedzi impulsowej zerami zawiera
wprawdzie widmo poprzednie w postaci próbek, ale pomiędzy tymi próbkami przybiera nie zawsze
poŜądany kształt:
>> plot(1:509,abs(H1016(1:509)));
>> hold on
>> plot(1:8:509, abs(H1016(1:8:509)),’r’);
takie samo porównanie dla fazy:
>> figure(2);
>> plot(1:509,angle(H1016(1:509)));
>> hold on
>> plot(1:8:509,angle(H1016(1:8:509)).*abs(H1016(1:8:509)),’r’);
MnoŜenie fazy przez amplitudę ma na celu stłumienie - tylko dla potrzeb tworzenia wykresu - tych
wartości fazy, które odpowiadają niewielkim wartościom amplitudy. MoŜna porównać wykresy faz po
zastosowaniu powyŜszego zabiegu i bez niego:
>> plot(1:8:509,angle(H1016(1:8:509)),’g’);
2.2.2 Projektowanie górnoprzepustowego filtru FIR 10. rzędu, o zadanej charakterystyce
określonej jak na poniŜszym rysunku
Jako pierwsza zostanie wykorzystana funkcja remez. Dlatego teŜ naleŜy rozpocząć od odpowiedniego
zdefniowana zadanej charakterystyki:
Projektowanie filtrów FIR, © P.Korohoda
15
>> F=[0,0.4,0.5,1];
>> M=[0,0,1,1];
przy takim sposobie wprowadzenia zadanej charakterystyki pasmo 0.4
÷
0.5 (przejściowe) nie będzie
brane pod uwagę w procesie optymalizacji; częstotliwość F=1 odpowiada połowie częstotliwości
próbkowania;
>> n=10;
>> b=remez(n,F,M) zwraca współczynniki filtru;
Nie podano średnika, by przyjrzeć się symetrii otrzymanych współczynników.
Filtry cechujące się taką symetrią współczynników mają liniową fazę:
>> freqz(b,1);
Analizując odpowiedź częstotliwościową zaprojektowanego filtru moŜna stwierdzić, Ŝe:
- w paśmie przepustowym ma on liniową fazę
- w paśmie zaporowym ma niewielką tłumienność, co moŜe stanowić dość istotną wadę.
W celu poprawienia tej cechy, kosztem pewnego pogorszenia charakterystyki w paśmie przepustowym,
moŜna wprowadzić odpowiednią funkcję wagową, np.:
>> W=[100,1];
>> b1=remez(n,F,M,W);
>> figure(2)
>> freqz(b1,1)
MoŜna teŜ poprawić dokładność aproksymacji w zakresie pasma przepustowego:
>> W=[1,100];
>> b2=remez(n,F,M,W);
>> figure(3)
>> freqz(b2,1)
Zadany filtr moŜna zaprojektować takŜe innymi metodami, np.:
>> b3=firls(n,F,M);
>> b4=fir2(n,F,M);
i obejrzeć ich odpowiedzi częstotliwościowe:
>> figure(4)
>> freqz(b3,1)
>> figure(5)
>> freqz(b4,1)
JeŜeli na podstawie zadanej charakterystyki przyjmie się dolną częstotliwość graniczną jako Fn=0.5, to
filtr ten moŜna zaprojektować takŜe i w inny sposób:
>> Fn=0.5;
>> b=fir1(n,Fn);
2.2.3 Porównanie wyniku projektowania filtru o zadanej charakterystyce amplitudowej z
wykorzystaniem procedur projektowych funkcji Matlab’a oraz „krok po kroku”
Najpierw pokazany zostanie przykład projektowania filtru „krok po kroku”.
NaleŜy zaprojektować filtr dolnoprzepustowy 10. rzędu, który aproksymowałby idealny filtr
dolnoprzepustowy o charakterystyce wygenerowanej w sposób następujący (z zachowaniem symetrii w
dziedzinie częstotliwości):
Projektowanie filtrów FIR, © P.Korohoda
16
>> H=[ones(1,33),zeros(1,63),ones(1,32)];
Odpowiedź impulsowa tego filtru:
>> h=real(ifft(H));
Odpowiedź impulsowa ma 128 próbek, odpowiadających kolejnym wartościom współczynników filtru.
PoniewaŜ filtr powinien być 10. rzędu, zatem okno skracające odpowiedź impulsową powinno być 11-
elementowe.
>> w=hamming(11)’; transponowanie do postaci wierszowej;
>> w128=[w(6:11),zeros(1,117),w(1:5)]; dostosowanie okna do zakresu indeksów
czasowych;
>> hw=h.*w128; zastosowanie okna do skrócenia
odpowiedzi impulsowej;
Otrzymano w ten sposób odpowiedź impulsową nieprzyczynowego filtru o zerowej fazie. W celu
wyznaczenia współczynników filtru przyczynowego o liniowej fazie przeprowadzane są następujące
przekształcenia:
>> h1=[hw(1:6),hw(124:128)]; wybranie niezerowych wartości próbek
>> h1=fftshift(h1)
>> freqz(h1,1);
Warto zwrócić uwagę, Ŝe zamiast dostosowywać okno do 128-elementowej odpowiedzi impulsowej z
zakresu indeksów od 0 do 127, moŜna było postąpić odwrotnie - najpierw wybrać 11 odpowiednich (to
znaczy jakich?) elementów odpowiedzi impulsowej i następnie przemnoŜyć je przez kolejne elementy
okna. Wynik byłby identyczny.
Dla porównania ten sam filtr zostanie zaprojektowany z uŜyciem procedury umieszczonej w funkcji
fir2. NaleŜy rozpocząć od odpowiedniego opisu charakterystyki. Trzeba zdefiniować wektor
częstotliwości cyfrowych z przedziału 0
÷
1 i odpowiadający mu wektor amplitud (tak jak w przykładzie
2.2.2). Oczywiście teraz charakterystyka będzie miała 2 razy mniej próbek (częstotliwość F=1
odpowiada połowie częstotliwości próbkowania).
>> F=0:1/64:1; 65-elementowy wektor częstotliwości;
>> H2=[ones(1,33),zeros(1,32)]; tylko próbki do częstotliwości znormalizowanej
F=1 (czyli “połowa”);
>> h2=fir2(10,F,H2); filtr 10. rzędu;
MoŜna teraz porównać wektory h1 i h2 oraz odpowiedzi częstotliwościowe:
>> figure
>> freqz(h2,1)
Projektowanie filtrów FIR, © P.Korohoda
17
2.2.4 Filtracja filtrem typu FIR z wykorzystaniem FFT (splotu kołowego)
W przykładzie tym podany zostanie fragment m-pliku, w którym naleŜy samodzielnie uzupełnić część
realizującą dla pojedynczego bloku odpowiedni algorytm (overlap-save lub overlap-add).
function y=over_add(x,h,NB);
% y=over_add(x,h,NB)
%
% x - ciag wejsciowy;
% h - odpowiedz impulsowa (Nh<<NB);
% NB - dlugosc bloku dla FFT - wtedy dlugosc wycinka "x": Nxb;
% y - wynik splatania;
% czesc inicjalizacyjna:
H=fft(h,NB);
Nh=length(h);
Nxb=NB-(Nh-1);
K=floor(length(x)/Nxb);
x=x(:)';
buf=zeros(1,Nh-1);
% algorytm wlasciwy:
for k=0:K-1,
tutaj nale
Ŝ
y wstawi
ć
własny fragment
end;
**********************************************************************************
function y=over_save(x,h,NB);
% y=over_save(x,h,NB)
%
% x - ciag wejsciowy;
% h - odpowiedz impulsowa (Nh<<NB);
% NB - dlugosc bloku dla FFT - wtedy dlugosc wycinka "x": Nxb;
% y - wynik splatania;
% czesc inicjalizacyjna:
H=fft(h,NB);
Nh=length(h);
Nxb=NB-(Nh-1);
K=floor(length(x)/Nxb);
if NB+Nxb*(K-1)>length(x), K=K-1;end; % na wszelki wypadek
x=x(:)';
% algorytm wlasciwy:
for k=0:K-1,
tutaj nale
Ŝ
y wstawi
ć
własny fragment
end;
Projektowanie filtrów FIR, © P.Korohoda
18
3 Zadania do wykonania
UWAGA - w ramach przygotowania się do ćwiczeń naleŜy opracować metodykę postępowania
zmierzającą do realizacji poniŜszych zadań.
1. Zademonstrować rezultat projektowania fitrów typu FIR za pomocą zbyt rzadko próbkowanego
widma (DFT jako przybliŜenie D-TFT).
2. Porównać dwa wybrane okna pod względem szerokości listka głównego i tłumienia dla listków
bocznych - sporządzić porównawcze wykresy okien oraz ich widm.
3. Zaprojektować filtr dolnoprzepustowy o zadanych parametrach za pomocą metody okien czasowych,
nie korzystając z procedur fir1 ani fir2; sprawdzić symulacyjnie działanie filtru; skonfrontować wynik
projektowania z wynikiem otrzymanym za pomocą funkcji fir1 lub fir2.
4. Zademonstrować zastosowanie algorytmu overlap-save (lub overlap-add) do filtracji liniowej za
pomocą pasmowego filtru FIR zdefiniowanego tylko w dziedzinie FFT. Wybrany algorytm naleŜy
zrealizować w języku pakietu Matlab w czasie zajęć, koncentrując się jedynie na poprawności
otrzymanego wyniku, a nie - jak poprzednio - na ilości wykonywanych operacji.
5. Zaprojektować filtr dolnoprzepustowy o zadanych parametrach za pomocą funkcji Matlab’a
realizującej algorytm Parksa-McClellana - sprawdzić działanie filtru.
6. Zaprojektować filtr górnoprzepustowy, korzystając z dowolnej metody dla filtrów
dolnoprzepustowych i wykorzystując symetrie filtrów kwadraturowych.
7. Zaprojektować filtr pasmowo-zaporowy, korzystając z dowolnej metody dla filtrów
dolnoprzepustowych, filtru górnoprzepustowego i liniowości transformacji Fouriera; porównać
otrzymany wynik z wynikiem uzyskanym dzięki procedurze Parksa-McCLellana (funkcja remez).
Projektowanie filtrów FIR, © P.Korohoda
19
Dodatek dla prowadzącego
Kompletne m-pliki do realizacji algorytmów overlap-add i overlap-save:
function y=over_add(x,h,NB);
% y=over_add(x,h,NB)
%
% x - ciag wejsciowy;
% h - odpowiedz impulsowa (Nh<<NB);
% NB - dlugosc bloku dla FFT - wtedy dlugosc wycinka "x": Nxb;
% y - wynik splatania;
% czesc inicjalizacyjna:
H=fft(h,NB);
Nh=length(h);
Nxb=NB-(Nh-1);
K=floor(length(x)/Nxb);
x=x(:)';
buf=zeros(1,Nh-1);
% algorytm wlasciwy:
for k=0:K-1,
n=1+Nxb*k:Nxb+Nxb*k;
xb=[x(n),zeros(1,Nh-1)];
XB=fft(xb);
YB=XB.*H;
yb=real(ifft(YB));
y(n)=[buf+yb(1:Nh-1),yb(Nh:Nxb)];
buf=yb(Nxb+1:NB);
end;
function y=over_save(x,h,NB);
% y=over_save(x,h,NB)
%
% x - ciag wejsciowy;
% h - odpowiedz impulsowa (Nh<<NB);
% NB - dlugosc bloku dla FFT - wtedy dlugosc wycinka "x": Nxb;
% y - wynik splatania;
% czesc inicjalizacyjna:
H=fft(h,NB);
Nh=length(h);
Nxb=NB-(Nh-1);
K=floor(length(x)/Nxb);
if NB+Nxb*(K-1)>length(x), K=K-1;end; % na wszelki wypadek
x=x(:)';
buf=zeros(1,Nh-1);
% algorytm wlasciwy:
for k=0:K-1,
n=1+Nxb*k:NB+Nxb*k;
% ny=Nh+Nxb*k:(Nh-1)+Nxb+Nxb*k; % po to by nie bylo roznicy w przesunieciu
ny=1+Nxb*k:Nxb+Nxb*k; % pojawi sie przesuniecie o "Nh-1"
xb=x(n);
XB=fft(xb);
YB=XB.*H;
yb=real(ifft(YB));
y(ny)=yb(Nh:Nh+Nxb-1);
end;