Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
Instrukcja do laboratorium z cyfrowego przetwarzania sygnałów
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”
Do sprawnego wykonania ćwiczenia nie jest konieczna wcześniejsza praktyczna znajomość nie wprowadzonych w
ramach poprzednich ćwiczeń funkcji pakietu MATLAB, jednak niezbędna jest dobra orientacja w materiale
przedstawionym w częściach 1 oraz 2 tej instrukcji oraz w zagadnieniach będących przedmiotem poprzednich ćwiczeń.
Dlatego też wskazane jest dokładne przeczytanie obu wymienionych części instrukcji oraz zanalizowanie podanych
przykładów. UWAGA: znajomość i zrozumienie części 1 i 2 oraz materiału z poprzednich ćwiczeń mogą zostać przez
prowadzącego skontrolowane w trakcie zajęć.
Zapoznanie się z materiałem zamieszczonym z dodatkach teoretycznych może ułatwić zrozumienie ćwiczeń.
UWAGA - dobre opanowanie materiału zawartego w tym ćwiczeniu jest bardzo ważne z punktu widzenia
dalszych ćwiczeń.
1
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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łą
τ
:
(
)
)
(
)
(
2
f
H
e
f
H
f
j
⋅
=
⋅
⋅
⋅
⋅
−
τ
π
(1)
Jest to równoważne stwierdzeniu, że faza określona jest wzorem opisującym funkcję liniową od zmiennej
, o
nachyleniu
f
− ⋅ ⋅
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
K
−
(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:
[ ]
(4)
H f
h n e
j
f n
n
( )
[ ]
=
⋅
− ⋅ ⋅ ⋅ ⋅
=−∞
+∞
∑
2
π
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:
(5)
(
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
1
0
π
π
)
2
π
)
]
co ostatecznie daje:
(6)
(
H f
h
h n
f n
n
M
( )
[ ]
[ ]
cos
=
+
⋅ ⋅
⋅ ⋅ ⋅
=
∑
0
2
2
1
π
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
zależy od wartości
. 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
H f
( )
f
±
π
. 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
2
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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
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
0
( )
(7)
H f
e
H f
j
f M
( )
( )
=
⋅
− ⋅ ⋅ ⋅ ⋅
2
0
π
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:
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)
3
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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
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
. Punkt
określa zatem wartość odniesienia, co oznacza, że amplituda wynosi w nim 0dB.
f
= 0
f
= 0
f
= 0
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
.
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ść
odpowiada częstotliwości próbkowania.
f
f
= 1
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,
4
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
5) okno Kaisera (najbardziej uniwersalne),
6) okno Czebyszewa.
Wzory opisujące powyższe okna w dziedzinie czasu można znaleźć w odpowiedniej literaturze.
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.
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
oznacza zadaną charakterystykę częstotliwościową - transformatę D-TFT - filtru idealnego, natomiast
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:
H
f
Z
( )
H f
( )
J
W f
H f
H
f
Z
1
2
0
0 5
=
⋅
−
∫
( )
( )
( )
.
df
(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)
5
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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
(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.
f
n
Inne bardziej popularne kryterium bazuje na bardzo prostej funkcji błędu zawierającej funkcję wag:
[
]
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
3
0 0 5
=
∈<
>
max
( )
, .
f
0 5
(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
(przedział
F obejmuje pasmo przepustowe i zaporowe). Jeżeli funkcja błędu
zdefiniowana według (12) wykazuje w tym
przedziale F co najmniej
0
≤
≤
f
,
E f
( )
M
+ 2
przerzuty, czyli:
(
)
E f
E f
E
E f
k
k
( )
(
)
max ( )
= −
= ±
=
−1
dla
f
f
f
f
M
k
0
1
1
F
≤
≤ ≤
∧
∈
+
K
,
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
oraz punktów brzegowych pasma przejściowego. Można też
przypomnieć, że rząd filtru
< 0 0 5
, .
>
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
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:
J
3
[
]
W f
H f
H
f
J
k
M
k
k
Z
k
k
( )
( )
( )
( )
:
, ,
,
⋅
−
= −
⋅
=
+
1
0 1
1
3
K
+ 1
(14)
Podstawienie w (14) za
na podstawie równania (6) dla wszystkich wartości k, daje układ równań, który można
zapisać macierzowo:
H f
k
( )
6
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
(
)
(
)
( )
(
)
(
)
( )
(
)
(
)
( )
(
)
(
)
(
)
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
L
L
M
M
O
M
M
L
L
M
M
M
h M
J
H
f
H
f
H
f
H
f
Z
Z
Z
M
Z
M
[
]
( )
( )
(
)
(
)
3
0
1
1
⎡
⎣
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎤
⎦
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
=
⎡
⎣
⎢
⎢
⎢
⎢
⎢
⎢
⎤
⎦
⎥
⎥
⎥
⎥
⎥
⎥
+
M
(15)
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
, do wyznaczenia odpowiedzi impulsowej i wartości
błędu
. Zgodnie z twierdzeniem Remeza, jeżeli punkty
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:
f
k
J
3
f
k
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
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
h n
[ ]
N
2
1
−
otrzyma się po następującej modyfikacji odpowiedzi impulsowej:
h n
n k
N
h n
1
2
2
[ ]
cos
[ ]
= ⋅
⋅ ⋅ ⋅
⎛
⎝⎜
⎞
⎠⎟ ⋅
π
(16)
7
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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 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
, czyli delta Kroneckera. Zatem odpowiedź impulsowa filtru górnoprzepustowego (FGP)
komplementarnego względem “1” do danego filtru dolnoprzepustowego (FDP), to:
h
n
d n
ALL
[ ]
[ ]
=
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
. 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
x
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
stosuje się
następujące zabiegi:
L
x
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
i
drugą, o długości
. Do kolejnych początkowych
L
x
L
h
− 1
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
stanowi odpowiedni fragment splotu liniowego.
L
x
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ść
i jest po prostu odrzucana, natomiast jako odpowiedni fragment splotu
liniowego pozostawia się część drugą, o długości
.
L
h
− 1
L
x
8
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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.
9
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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.
Wybrane funkcje pakietu Matlab - część 2
funkcje realizujące procedury projektowania filtrów FIR
7Nazwa 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)
10
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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.
11
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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:
>> 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;
12
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
>> 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):
>> 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;
13
Projektowanie filtrów FIR, ćwiczenie 2-2 © P.Korohoda
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)
14