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;