2IA PS2 2012 2013 05 id 32602 Nieznany (2)

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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

background image

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


Wyszukiwarka

Podobne podstrony:
2IA PS2 2012 2013 04 B id 32601 Nieznany (2)
07 05 2013 odwiert (1)id 6788 Nieznany
plan lekcji 2012 2013 nowy id 3 Nieznany
cw PAiTS 05 id 122324 Nieznany
NAI2006 05 id 313056 Nieznany
matma dyskretna 05 id 287941 Nieznany
cwiczenie 05 id 125057 Nieznany
28 04 2013 cw id 31908 Nieznany
ef 271 4 2012 zal nr 2 id 15072 Nieznany
8 lect8 2013 stud id 46719 Nieznany (2)
lab pwsp 05 id 258618 Nieznany
Zestaw 05 id 587909 Nieznany
26429 05 id 31506 Nieznany
25 3 2013 traduction id 31052 Nieznany (2)

więcej podobnych podstron