Instrukcja do laboratorium z cyfrowego przetwarzania sygnałów
Ćwiczenie 9
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)
3 Zadania do wykonania
Na instrukcję składają się następujące części:
1 Wybrane zagadnienia z zakresu DSP
2 Projektowanie filtrów cyfrowych typu FIR za pomocą pakietu MATLAB
3 Zadania do wykonania
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 tej instrukcji oraz w materiale poprzednich ćwiczeń.
W razie niejasności należy skonsultować się przed zajęciami ( tzn. na przykład w terminie konsultacji ) z prowadzącym, bezpośrednio lub poprzez e-mail:
korohoda@uci.agh.edu.pl
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łą :
(1)
Jest to równoważne stwierdzeniu, że faza określona jest wzorem opisującym funkcję liniową od zmiennej , o nachyleniu :
(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 .
Warto przypomnieć, że pomiędzy częstotliwością cyfrową () stosowaną w przypadku D-TFT i indeksem częstotliwościowym, typowym dla DFT, zachodzi następujący związek:
(3)
gdzie oznacza długość ciągu poddawanego transformacji DFT.
Odpowiedź częstotliwościowa - czyli transformata D-TFT - filtru o odpowiedzi impulsowej określona jest wzorem:
(4)
Dla skończonej, nieparzystej odpowiedzi impulsowej o długości , symetrycznej względem punktu , wzór ten sprowadza się do postaci:
(5)
co ostatecznie daje:
(6)
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 . 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 . 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 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ą:
(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:
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 :
(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 .
Zatem odpowiedź impulsowa po zastosowaniu okna (8), będzie określona następująco:
(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.
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.
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.
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:
(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 :
(11a)
dla nieparzystej długości ciągu :
(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 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.
Inne bardziej popularne kryterium bazuje na bardzo prostej funkcji błędu zawierającej funkcję wag:
(12)
Wskaźnik jakości, zwany też kryterium Czebyszewa, określa się wówczas następująco:
(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 przerzuty, czyli:
dla ,
to istnieje dokładnie jeden filtr rzędu 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ż . 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 oznacza współczynników filtru, czyli odpowiedź impulsową o długości .
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:
(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:
(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:
0) równanie (15) jest rozwiązywane dla wstępnego układu punktów ,
1) dla otrzymanego filtru wyznaczana jest charakterystyka częstotliwościowa, której lokalne ekstrema określają kolejny układ punktów ,
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 otrzyma się po następującej modyfikacji odpowiedzi impulsowej:
(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 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:
(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 . W kolejnym kroku blok o łącznej długości 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 . W celu sprowadzenia wyniku do postaci odpowiadającej fragmentowi splotu liniowego o długości stosuje się następujące zabiegi:
a) w algorytmie overlap-add:
Wynik splotu kołowego, o długości , dzielony jest na dwie części - pierwszą, o długości i drugą, o długości . Do kolejnych początkowych 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 - 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.
b) w algorytmie overlap-save:
Wynik splotu kołowego, o długości , 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 .
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.
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) |
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. |
Szczegół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.
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;
>> 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;
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)
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;
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).
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;
Projektowanie filtrów FIR, © P.Korohoda
5
Wyszukiwarka