Projektowanie filtrów cyfrowych Butterwortha i Czebyszewa

Klasyczne filtry Butterwortha i Czebyszewa są to bądź bierne obwody LC bądź układy aktywne z ujemnym sprzężeniem zwrotnym zbudowane na wzmacniaczach operacyjnych. Cyfrowym odpowiednikiem takich filtrów są filtry również z ujemnym sprzężeniem zwrotnym, zwane filtrami rekurencyjnymi, bądź też filtrami z nieskończoną odpowiedzią impulsową (ang. IIR - Infinite Impulse Response).

Najczęstszym sposobem formułowania problemu zaprojektowania filtru jest podanie jego częstotliwości granicznej, przy której wzmocnienie spada o 3 dB (71% w skali liniowej), oraz rzędu filtru. Jest to jednak sformułowanie nieprecyzyjne, bo nic nie mówi o tłumieniu w paśmie zaporowym i zakłada z góry rząd filtru. Znacznie bardziej precyzyjnym i lepiej dostosowanym do potrzeb praktycznych jest określenie charakterystyki filtru przez podanie dwóch punktów: końca pasma przepustowego (punkt P) i początku pasma zaporowego (punkt R). Punkty te znajdują się na charakterystyce amplitudowej filtru (wzmocnienie w funkcji częstotliwości) i wyznaczają pole, w którym mieścić się musi charakterystyka projektowanego filtru. Na rysunku obok obie charakterystyki, zielona i niebieska spełniają założone wymagania, ponieważ nie przebiegają przez zakreskowane obszary zabronione wyznaczone punktami P i R. Warto zauważyć, że o ile charakterystyka zawsze przechodzi przez punkt P, o tyle nie musi przechodzić przez punkt R, może przechodzić pod nim, wtedy możemy powiedzieć, że filtr z zapasem spełnia wymagania projektowe. Z położenia punktów P i R wynika rząd filtru. Oczywiście im bardziej stromo musi opadać charakterystyka amplitudowa między punktami P i R, tym rząd filtru będzie większy.

Tak więc danymi początkowymi do obliczenia filtru są:

• Częstotliwość końca pasma przepustowego (f_P),
• Największe dopuszczalne tłumienie w paśmie przepustowym (A_P), dla filtru Czebyszewa ten parametr określa jednocześnie dopuszczalne zafalowanie charakterystyki amplitudowej w paśmie przepustowym,
• Częstotliwość początku pasma zaporowego (f_R),
• Wymagane minimalne tłumienie w paśmie zaporowym (A_R).

W celu praktycznego pokazania omawianego sposobu projektowania filtrów, podając wzory i teorię, będę jednocześnie projektował konkretny filtr dolnoprzepustowy o następujących parametrach:

A_P = 0,5 dB, f_P = 3 kHz, A_R = 20 dB, f_R = 7 kHz.

Obliczenie rzędu filtru

Pierwszym krokiem jest obliczenie dwóch współczynników określających nasz filtr. Są to współczynnik selektywności k, oraz współczynnik dyskryminacji d. Współczynnik selektywności określa jak blisko siebie położone są częstotliwości punktów P i R, gdy punkty te zbliżają się w poziomie do siebie, to współczynnik dąży do 1. Współczynnik ten określony jest wzorem:

  (1)

Współczynnik dyskryminacji natomiast określa wzajemną relację między maksymalnym tłumieniem w paśmie przepustowym i minimalnym tłumieniem w paśmie zaporowym. Obliczamy go ze wzoru:

  (2)

Dla podanych przykładowych wartości liczbowych po podstawieniu ich do wzoru otrzymamy: k = 0,42857 i d = 0,035107. Rząd filtru Butterwortha obliczamy ze współczynników k i d z następującego wzoru:

  (3)

Rząd filtru Czebyszewa natomiast określa się z podobnego, ale jednak nieco innego wzoru:

  (4)

Po podstawieniu danych przykładowych otrzymamy dla filtru Butterwortha rząd 3,952968, dla filtru Czebyszewa rząd 2,711062. W rzeczywistości rząd filtru nie może być ułamkowy, musi być liczbą naturalną. Otrzymane liczby należy więc zaokrąglić w górę. Zaokrąglenie w dół spowodowałoby, że charakterystyka filtru przechodziłaby ponad punktem R, a więc filtr nie spełniałby założeń. Dlatego dla przykładowego filtru Butterwortha przyjmujemy rząd 4, dla filtru Czebyszewa rząd 3. Rząd filtru Czebyszewa dla tych samych danych będzie zawsze mniejszy lub równy rzędowi filtru Butterwortha.

Częstotliwość graniczna filtru

Częstotliwość graniczna inaczej definiowana jest dla filtru Butterwortha i filtru Czebyszewa. W dolnoprzepustowym filtrze Czebyszewa jest to najwyższa częstotliwość, dla której tłumienie filtru nie przekracza założonych zafalowań charakterystyki. Jest to więc dokładnie częstotliwość punktu P na wykresie. Dla filtru Butterwortha natomiast częstotliwość graniczna, to taka, przy której następuje spadek wzmocnienia o 3 dB względem sygnału stałego (o częstotliwości 0). Częstotliwość ta jest większa od f_P jeżeli A_P jest mniejsze niż 3 dB. Jeżeli tłumienie w punkcie Pwynosi 3 dB, to oczywiście częstotliwość graniczna jest częstotliwością punktu P. Częstotliwość graniczną filtru Butterwortha można określić ze wzoru:

   (5)

Dla filtru przykładowego w wersji Czebyszewa częstotliwość graniczna jest równa f_P a więc wynosi 3000 Hz. Dla filtru Butterwortha otrzymamy po podstawieniu do wzoru (5) częstotliwość graniczną równą 3902,27 Hz.

Filtr prototypowy

Po określeniu rzędu i częstotliwości granicznej obliczamy filtr prototypowy, to znaczy taki filtr analogowy którego częstotliwość graniczna wynosi 1 Hz. Transmitancja tego filtru to funkcja operatora różniczkowania s o postaci:

   (6)

Mianowniki transmitancji filtrów prototypowych zostały stablicowane. Postać mianownika dla filtru Butterwortha zależy wyłącznie od rzędu filtru, dla filtru Czebyszewa dodatkowym parametrem jest założone zafalowanie charakterystyki amplitudowej (równe tłumieniu A_P). Dla filtru przykładowego odczytujemy z tablic:

filtr Butterwortha

(s² + 0,765366864s + 1)·(s² + 1,847759065s + 1)

filtr Czebyszewa

(s² + 0,62646s + 1,14245)·(s + 0,62646)

Jak widać mianownik transmitancji filtru M(s) jest wielomianem zmiennej s. Miejsca zerowe tego mianownika są zwanebiegunami filtru. Bieguny filtru są liczbami zespolonymi o ujemnej części rzeczywistej (gdyby któryś biegun miał dodatnią część rzeczywistą, to filtr byłby niestabilny). Filtr ma tyle biegunów ile wynosi jego rząd. Filtry o rzędzie parzystym mają bieguny zespolone sprzężone parami. W filtrze o rzędzie nieparzystym jeden biegun jest ujemny rzeczywisty, reszta jest zespolona, sprzężona parami. Bieguny można bardzo prosto obliczyć wyliczając pierwiastki trójmianów w nawiasach. Dla przykładowego filtru Butterwortha otrzymamy dwie pary biegunów zespolonych sprzężonych:

s₁₂ = −0,382683±j0.923880
s₃₄ = −0,923880±j0.382683

Dla filtru Czebyszewa otrzymamy jeden biegun rzeczywisty i parę biegunów zespolonych sprzężonych:

s₁ = −0,62646
s₂₃ = −0,31323±j1,021928

Skalowanie częstotliwości i przekształcenie nieliniowe

Jak napisałem wcześniej, filtr prototypowy obliczony jest dla częstotliwości granicznej 1 Hz. Należy teraz przeskalować jego bieguny, tak aby odpowiadały założonej częstotliwości granicznej. Oprócz tego należy dokonać przekształcenia nieliniowego. Polega ono na potraktowaniu przeskalowanej częstotliwości funkcją tangens rozciągniętą w taki sposób, że wartość π na osi x przesuwa się do częstotliwości próbkowania sygnału. Przekształcenie nieliniowe jest elementem transformacji z dziedziny analogowego operatora s (różniczkowanie) do dziedziny cyfrowego operatora z (opóźnienie o 1 takt zegara). Przy projektowaniu filtrów cyfrowych łączy się operację skalowania i przekształcenia nieliniowego w jedną, obliczając współczynnik skalujący W z następującego wzoru:

   (7)

Dla filtru Butterwortha i Czebyszewa współczynnik W liczymy oddzielnie, ponieważ z reguły filtry te mają inne częstotliwości graniczne (chyba, że A_P wynosi 3 dB, wtedy obie częstotliwości będą równe). Dla naszego przykładowego filtru Butterwortha f_g wynosi 3902,27 Hz, a współczynnik W = 0,570758. Dla filtru Czebyszewa f_g = 3000 Hz, W = 0,43405608.

Po obliczeniu współczynnika W transformujemy bieguny filtru prototypowego. Jest to bardzo prosta operacja, dla filtrudolnoprzepustowego po prostu mnożymy W przez biegun. Dla filtru górnoprzepustowego dzielimy W przez biegun. Po przemnożeniu (dla filtru dolnoprzepustowego) otrzymamy następujące bieguny:

Filtr Czebyszewa:

s₁₂ = −0,135959±j0,443573
s₃ = −0,271917

Przejście do transmitancji cyfrowej (operatora z)

Matematycznie ścisłe przejście z operatora analogowego s do cyfrowego operatora z dane jest wzorem:

   (8)

gdzie T_d to okres próbkowania sygnału. Wydaje się to prostą zależnością, dopóki nie zauważymy, że wykładnik potęgi jest liczbą zespoloną, co znacznie komplikuje obliczenia. Dlatego w praktyce stosuje się pewne przybliżenie wzoru teoretycznego zwane przekształceniem dwuliniowym lub bilinearnym. Przybliżenie to wygląda następująco:

   (9)

Gdzie s' jest to przetransformowany wcześniej biegun filtru analogowego, a z biegun filtru cyfrowego. Takie przekształcenie jest znacznie prostsze obliczeniowo i daje dobre rezultaty praktyczne. Przekształcenie zachowuje sprzężenie biegunów, to znaczy że para biegunów zespolonych sprzężonych s pozostaje parą biegunów zespolonych sprzężonych z. Dzięki temu dla pary biegunów sprzężonych można wykonać tylko połowę obliczeń. Podobnie biegun rzeczywisty pozostaje po przekształceniu biegunem rzeczywistym. Dla przykładowego filtru bieguny po przekształceniu dwuliniowym są następujące:

Filtr Czebyszewa:

s₁₂ = 0,795271±j0,372823
s₃ = 0,760628

Wyliczenie transmitancji

Wyliczenie transmitancji polega na wstawieniu obliczonych biegunów do mianownika, oraz wstawieniu odpowiednich zer (miejsc zerowych licznika) do licznika. Do obliczeń w MatLabie wymnażamy wszystkie bieguny (podobnie zera) przez siebie aż do uzyskania w liczniku i mianowniku wielomianów stopnia N równego rzędowi filtru. Do praktycznej realizacji filtru korzystniej jest pozostawić licznik i mianownik jako iloczyn wielomianów stopnia pierwszego i drugiego. Odpowiada to wykonaniu filtru jako kaskadowego połączenia filtrów pierwszego i drugiego rzędu. Takie rozwiązanie jest korzystne ponieważ znacznie zmniejsza wrażliwość filtru na zaokrąglenia współczyników wzmacniaczy. Szczególnie warto więc je zastosować w przypadku pracy w arytmetyce stałoprzecinkowej. Wtedy każdy biegun rzeczywisty odpowiada blokowi pierwszego rzędu, a każda para zespolonych biegunów sprzężonych – blokowi drugiego rzędu.

Blok I rzędu

Blok taki powstaje z bieguna rzeczywistego i ma jedno zero w punkcie −1 (filtr górnoprzepustowy ma zero w punkcie 1). Licznik transmitancji ma zatem postać

Mianownik natomiast (zakładając że biegun wynosi c)

Od transmitancji można przejść do schematu blokowego, korzystając z definicji transmitancji, jest to bowiem stosunek transformaty sygnału wyjściowego Y(z) do transformaty sygnału wejściowego X(z). Dla bloku pierwszego stopnia mamy zatem:

Z ostatniego wzoru bezpośrednio można narysować schemat blokowy bloku filtru I rzędu.

Przykładowy filtr Butterwortha nie posiada biegunów rzeczywistych, a więc i bloków I rzędu. Filtr Czebyszewa posiada jeden biegun rzeczywisty, odpowiadający mu blok ma transmitancję

Blok II rzędu

Blok ten powstaje z pary biegunów zespolonych sprzężonych i ma podwójne zero w punkcie -1 (filtr górnoprzepustowy ma podwójne zero w punkcie 1). Zatem licznik transmitancji tego bloku ma postać

Mianownik natomiast przyjmuje postać (przy parze biegunów a±jb)

Podobnie jak przy bloku pierwszego rzędu korzystamy z definicji transmitancji, aby uzyskać wzór definiujący postać schematu blokowego

Przykładowy filtr Butterwortha posiada dwie pary biegunów zespolonych sprzężonych, będzie zatem miał dwa bloki drugiego rzędu. Dla bieguna z₁₂ transmitancja wyniesie

Dla bieguna z₃₄ transmitancja wyniesie

Przykładowy filtr Czebyszewa posiada jedną parę biegunów zespolonych sprzężonych. Transmitancja odpowiadającego mu bloku jest następująca

Do weryfikacji obliczeń w MatLabie niezbędna jest kompletna transmitancja filtru jako całości. Uzyskujemy ją wstawiając wszystkie bieguny jako miejsca zerowe mianownika, oraz N-krotne zero w punkcie -1 jako miejsce zerowe licznika, a następnie wymnażając. Równoważnym sposobem jest wymnożenie przez siebie transmitancji wszystkich bloków filtru. Tramsmitancja przykładowego filtru Butterwortha ma postać

Transmitancja filrtu Czebyszewa natomiast

W programie MatLab można sprawdzić charakterystykę amplitudową przykładowego filtru wydając następujące polecenie (wartości liczbowe dla przykładowego filtru Butterwortha):

freqz([1 4 6 4 1],[1 -2.555369 2.622493 -1.245102 0.228714],10:10:20000,44100)

Pierwszy wektor zawiera kolejne współczynniki wielomianu licznika transmitancji, drugi kolejne współczynniki mianownika.

Charakterystyki filtru przykładowego

W celu weryfikacji metody projektowania sprawdziłem charakterystyki otrzymanych filtrów korzystając z funkcji freqz()programu MatLab. Punkty P i R z założeń projektowych są zaznaczone kolorem czerwonym. Pewnym zaskoczeniem może być fakt, że w pasmie przenoszenia wzmocnienie nie wynosi 0 dB, to można jednak łatwo skorygować umieszczając tłumik o odpowiednim tłumieniu (50 dB dla filtru Butterwortha i 45,3 dB dla filtru Czebyszewa). Praktycznie realizując filtr w arytmetyce stałoprzecinkowej warto podzielić tłumik na kilka tłumików rozmieszczonych przed filtrem, między jego blokami oraz za filtrem, tak aby ich sumaryczne tłumienie było równe wymaganemu. W ten sposób zapobiega się powstaniu przepełnienia arytmetycznego w blokach filtru, oraz pogarszaniu stosunku sygnału do szumu kwantyzacji.

Charakterystyka amplitudowa przykładowego filtru Butterwortha w zakresie od 0 do 3,5 kHz.

Charakterystyka amplitudowa przykładowego filtru Butterwortha wokół częstotliwości 7 kHz.

Charakterystyka amplitudowa przykładowego filtru Czebyszewa w zakresie od 0 do 20 kHz (logarytmiczna skala częstotliwości).

Charakterystyka amplitudowa przykładowego filtru Czebyszewa w zakresie od 0 do 3,5 kHz.

Charakterystyka amplitudowa przykładowego filtru Czebyszewa wokół częstotliwości 7 kHz.