1. Wstęp

Filtracja danych cyfrowych jest najwcześniej pojętą dyscypliną w dziedzinie przetwarzania sygnałów cyfrowych. Początki filtracji cyfrowej sięgają wczesnych lat pięćdziesiątych, gdy wzrastająca dostępność komputerów zaowocowała pracami nad wygładzaniem spróbkowanych sygnałów dyskretnych i analizą dyskretnych systemów sterowania. Obecnie filtracja cyfrowa jest stosowana tak szeroko, że liczba pozycji literatury na ten temat jest znacznie liczniejsza niż w innych dziedzinach cyfrowego przetwarzania sygnału.

Mianem filtru określamy, takie urządzenie elektroniczne, które kształtuje widmo sygnału w narzucony przez projektanta sposób tzn. usuwa z widma niepożądane jego składniki pozostawiając pasmo częstotliwości użytecznych nietknięte. W tej pracy zostaną omówione stacjonarne filtry liniowe. Właściwość stacjonarności oznacza niezależność parametrów układu od czasu. Natomiast termin liniowy oznacza, że układ spełnia zasadę superpozycji, czyli odpowiedź tego układu jest liniową kombinacją sygnałów wejściowych. W dziedzinie częstotliwości oznacza to, że układ nie dodaje żadnych składowych od widma sygnału, zmienia jedynie amplitudę i fazę istniejących składowych przebiegu wejściowego. Podobnie jak technice analogowej, układ cyfrowy z powodu niekorzystnych zjawisk, wynikających ze skończonej długości słowa binarnego użytego przy realizacji filtru, przyjęcie, że filtr cyfrowy jest układem liniowym jest pewnym przybliżeniem prawdziwego zachowania się systemu. Istnieje wiele nieliniowych filtrów cyfrowych, które najczęściej dokonują obróbki sygnału na podstawie jego analizy statystycznej. W pracy ten temat zostanie pominięty. Ponadto wszystkie metody projektowania filtrów, zagadnienia implementacji, jak i same zaimplentowane filtry są filtrami jednowymiarowymi, to znaczy, że sygnał wejściowy i wyjściowy jest funkcją jednej zmiennej.

Podstawową przewagą filtrów cyfrowych nad filtrami analogowymi jest możliwość otrzymywania bardzo dobrych parametrów charakterystyk częstotliwościowych. Poziomy tłumienia osiągane w filtrze cyfrowym są dużo większe niż w filtrze analogowym o tym samym rzędzie. Pasmo przejściowe w filtrach cyfrowych jest bardzo wąskie. Ponadto w filtrach cyfrowych można uzyskać pasmo przejściowe pozbawione jakichkolwiek zafalowań. Kolejna zaletą filtrów cyfrowych jest możliwość uzyskania filtru o liniowej fazie w paśmie przepustowym, co powoduje, że taki filtr nie wprowadza zniekształcenia fazy w sygnale użytecznym. Parametry filtrów cyfrowe są stałe w czasie, ponieważ w układach cyfrowych nie istnieje problem starzenia się elementów. Cechuje je także większa elastyczność w procesie projektowania. Konstruktor zmieniając współczynniki filtru i/lub strukturę filtru może uzyskać filtr o lepszych parametrach. Ponadto wszystkie te czynności wykonuje najczęściej bez zmiany strony sprzętowej urządzenia. W wyniku tych własności czas projektowania jest dosyć krótki. Ponadto istnieje możliwość projektowania cyfrowych filtrów adaptacyjnych dostosowujących swoje charakterystyki do aktualnego sygnału na wejściu urządzenia. Filtry cyfrowe mają pewne wady. Ich pasmo przenoszenia jest ograniczone do połowy częstotliwości próbkowania użytej w układzie. Pasmo przenoszenia filtrów analogowych nie zna tego ograniczenia. Jest ono szczególnie widoczne podczas oglądanie wykresów charakterystyk częstotliwościowych. Oś częstotliwości na wykresach filtrów analogowych jest najczęściej podawana w skali logarytmicznej, zaś dla filtrów cyfrowych jest ona zazwyczaj liniowa. W momencie można zauważ konieczność używania filtrów analogowych podczas konstruowania filtru cyfrowego. W torze układu przed przetwornikiem analogowo-cyfrowym musi zostać włączony filtr analogowy, aby usunąć z przetwarzanego przebiegu wszystkie składowe powyżej częstotliwości próbkowania. Kolejną wadą filtrów cyfrowych jest przedział dostępnych amplitud. Stosunek sygnału do szumu w układach analogowych jest znacznie mniejszy niż w układach cyfrowych. Dla przykładu dla przetwornika analogowo cyfrowego 16-bitowego wartość skuteczna szumu wynosi około 9 μV, natomiast dla typowego układu zbudowanego przy użyciu wzmacniaczy operacyjnych 2μV. Ponadto każda operacja mnożenia w układzie cyfrowym powoduje powstawanie dodatkowego szumu, z powodu konieczności obcięcia iloczynu, aby zmieścił się w rejestrze, którego długość jest ustalona. Szybkość działania układów analogowych jest dużo większa. Załóżmy, że na dokonanie filtracji na procesorze sygnałowym taktowanym zegarem 100MHz potrzeba 100 taktów zegara. Daje to maksymalną częstotliwość próbkowania rzędu 1MHz. Tymczasem wzmacniacze operacyjne mogą wzmacniać sygnały o częstotliwościach rzędu setek MHz.

Tradycyjne filtry cyfrowe mogą występować jako jeden z dwóch typów: filtry o skończonej odpowiedzi impulsowej FIR i filtry o nieskończonej odpowiedzi impulsowej IIR. Pierwsze z nich obliczają próbkę wyjściową korzystając tylko z przeszłych próbek wejściowych, z kolei drugie używają także poprzednich wartości wyjścia układu. Transmitancja filtru FIR zawiera tylko zera, tymczasem transmitancja fitru IIR zawiera zarówno zera, jak i bieguny. Dlatego filtry FIR są zawsze stabilne, natomiast zapewnienie stabilności filtru typu IIR należy do projektanta. Ponadto, jeśli projektujemy filtr FIR możemy uzyskać charakterystykę o liniowej fazie poprzez dobranie odpowiedniej symetrii współczynników filtru. Dla filtru IIR nie mamy takiej możliwości. Największą zaletą filtru typu IIR jest ich duża efektywność. Podobne parametry charakterystyki częstotliwościowej można przy zastosowaniu struktury filtru zawierającej dużo mniej mnożników. Stosowane rzędy filtrów typu IIR wynoszą, w zależności od rodzaju, co najwyżej 30, z kolei projektuje się filtry typu IIR rzędu kilkuset. Ponadto filtry FIR mają dużo większe wymagania, jeśli chodzi o ilość pamięci potrzebnej do zaimplementowania filtru. Kiedy projektujemy filtry IIR musi zwrócić szczególną uwagę na dobór struktury filtru, ponieważ filtry tego typu są szczególnie wrażliwe na błędy kwantyzacji oraz prawdopodobieństwo wystąpienia błędów nadmiaru w strukturach tego filtru jest dość duże w stosunku do struktur filtru typu FIR. Ponadto stopień trudności analizy szumu kwantyzacji dla filtrów typu IIR jest bardzo duży. Dla filtrów typu IIR istnieje możliwość wykorzystania metod projektowania filtrów analogowych, które są bardzo rozwinięte.

Przed rozpoczęciem projektowania filtru cyfrowego należy określić parametry, jakie ma spełnić wykonany filtr. Na te parametry składają się: granice pasm filtru, szerokość pasm przejściowych pomiędzy kolejnymi pasmami, wymaganie tłumienie w poszczególnych pasmach zaporowych, dopuszczalne zafalowania w kolejnych pasmach zaporowych oraz kształt charakterystyki fazowej. Mając zdefiniowane wymagania nałożone na charakterystykę częstotliwościową musimy wybrać, który typ filtru będzie lepiej je spełniał. Należy też wziąć pod uwagę rodzaj formy w jakiej będzie realizowany filtr cyfrowy: może być programem komputerowym, procesorem sygnałowym lub specjalizowanym układem scalonym typu ASIC. Kolejnym krokiem jest wybór metody otrzymania współczynników filtru. Posiadanie wiedzy na temat właściwości metod projektowania filtrów umożliwia otrzymanie filtru o jak najlepszych parametrach w jak najkrótszym czasie. Jednak otrzymane w tym kroku współczynniki filtru nie gwarantują otrzymania o dobrych parametrach. Skutki implementacji filtru na sprzęcie o skończonej długości słowa, a zwłaszcza używającego arytmetyki stałoprzecinkowej, mogą spowodować pogorszenie lub wręcz utratę dobrych parametrów filtru. Dlatego konieczny jest dobór odpowiedniej struktury filtru, która powinna zminimalizować występujący szum kwantyzacji, charakteryzować się małą wrażliwością na skutki dyskretyzacji współczynników oraz zapewnić, jak największą efektywność wykonanego filtru poprzez, jak największe zmniejszenie liczby operacji matematycznych potrzebnych do uzyskania wyniku filtrowania. Następnie należy dokonać implementacji filtru w postaci programu lub układu scalonego. Tak otrzymany filtr należy przetestować, czy spełnia nałożone na niego wymagania. Jeśli ich nie spełnia należy wrócić do poprzednich punktów, aby w toku procedury iteracyjnej uzyskać filtr spełniający nałożone na niego wymagania.

W pracy skupiono się na dwóch punktach: metodach projektowania i błędach kwantyzacji. Rozdział 2 dotyczy metod projektowania filtrów FIR, zaś rozdział 3 omawia metody projektowania filtrów IIR. W rozdziale 4 omówiono skutki użycia skończonego słowa do implementacji filtrów cyfrowych. Rozdział 5 omawia stworzone przy pomocy pakietu MATLAB interfejsy GUI, dzięki którym możemy projektować filtry IIR i FIR, oraz dokonać analizy wpływu kwantyzacji. Natomiast w rozdziale 6 przedstawiono realizację filtrów cyfrowych przy pomocy programu Code Composer Studio. Filtry zostały zrealizowane na zmiennoprzecinkowym procesorze sygnałowym TMS 320C6711.

2. Projektowanie filtrów o skończonej odpowiedzi impulsowej

2.1. Wprowadzenie do metod projektowania filtrów FIR

Główna zaletą filtrów FIR jest możliwość realizacji dokładnie liniowej charakterystyki fazowej filtru. Dlatego prawie wszystkie metody projektowania tych filtrów zajmują się filtrami o tej własności. Jeżeli założymy, że projektujemy filtr FIR o liniowej fazie, procedura projektowa sprowadza się do problemów aproksymacji funkcją o wartościach rzeczywistych, gdzie współczynniki są optymalizowane jedynie pod względem uzyskania odpowiedniej charakterystyki amplitudowej. Filtr FIR N-tego rzędu jest określony poprzez N + 1 współczynników i jego charakterystykę częstotliwościową można zapisać w następujący sposób:

(2.1)

Wyłączając przed sumę wyrażenie odpowiadające opóźnieniu o NT/2 (połowa długości rejestru zawierającego poprzednie próbki wejściowe) otrzymujemy:

(2.2)

Dla współczynników b_r istnieją dwa możliwe rodzaje symetrii, które prowadzą do sytuacji w której, poza wyłączonym przed sumę czynnikiem o liniowej fazie, nie ma innych zależnych od częstotliwości czynników wpływających na fazę charakterystyki. Są to:

br = bN-r br = - bN-r

(2.3)

Wprowadzając te warunki do poprzedniego wzoru i odpowiedni wzór Eulera, otrzymujemy zamiast sumy ważonych eksponent zespolonych sumę ważonych rzeczywistych funkcji trygonometrycznych: odpowiednio kosinusoid dla symetrycznych współczynników i sinusoid dla asymetrycznych współczynników.

(2.4)

Sumy w powyższych równaniach są czysto rzeczywiste, więc nie wpływają one na kształt charakterystyki fazowej. Stanowią one o kształcie charakterystyki amplitudowej. Należy zauważyć, że filtry asymetryczne, poza wyrażeniem stanowiącym o liniowej fazie filtru, zawierają także współczynnik j, który powoduje niezależne od częstotliwości przesunięcie fazy o 90 stopni. Takie filtry są więc szczególne często stosowane do realizacji integratorów, układów różniczkujących czy przetworników Hilberta w postaci filtrów FIR.

Oprócz rozróżnienia ze względu na rodzaj symetrii ważnym własnością filtrów jest parzystość rzędu filtru. Ze względu na parzystość i rodzaj symetrii rozróżniamy cztery typy filtrów FIR:
Typ I - symetryczny o parzystej długości,
Typ II - symetryczny o nieparzystej długości,
Typ III - asymetryczny o parzystej długości,
Typ IV - asymetryczny o nieparzystej długości.

Niestety na prawie każdy typ filtru są narzucone pewne ograniczenia na kształt charakterystyki filtru:
Typ I - bez ograniczeń,
Typ II - H₀(π/T) = 0, czyli można projektować filtry dolno- i pasmowoprzepustowe,
Typ III - H₀(0) = 0 i H₀(π/T) = 0, czyli można projektować tylko filtry pasmowoprzepustowe,
Typ IV - H₀(0) = 0, czyli można projektować filtry górno- i pasmowoprzepustowe.

Rys. 2.1. Typy filtru FIR.

Ponadto filtry o liniowej fazie posiadają ciekawą właściwość dotyczącą położenia zer. Transmitancja filtru FIR w dziedzinie z dana jest wzorem:

(2.5)

Jeśli podstawimy z^-1 w miejsce z i pomnożymy obie strony równania przez z ^{-(M - 1)}, to otrzymamy:

(2.6)

Oznacza to, jeśli z₁ jest zerem filtru FIR, to także 1/z₁ jest także zerem tego filtru. Co więcej, jeśli współczynniki filtru są liczbami rzeczywistymi, to każde zero zespolone musi posiadać w transmitancji swój odpowiednik sprzężony. Czyli jeśli z₁ jest zerem filtru, to z₁^* też jest zerem tego filtru. W konsekwencji tego co powiedzieliśmy wcześniej także 1/z₁^* musi być zerem tego filtru. Sytuacje ilustruje poniższy rysunek 2.2.

Rys. 2.2. Położenie zer w filtrze FIR o liniowej fazie.

Dla filtrów FIR nie istnieje żadna metoda projektowania filtrów w dziedzinie czasu, która może stanowić podstawę projektowania odpowiedników cyfrowych. Dlatego przy ich projektowaniu używa się odpowiednich formuł matematycznych. Zależnie od wybranej metody aproksymacji uzyskuje się bardziej lub mniej złożone algorytmy, iteracyjne lub analityczne, które w większości przypadków wymagają użycia komputera.

Celem wszystkich metod projektowania filtrów FIR jest jak najlepsza aproksymacja pożądanej charakterystyki częstotliwościowej skończoną liczbą współczynników filtru. Punktem startowym dla wszystkich tych metod jest założenie idealnej charakterystyki lub specyfikacja dopuszczalnych odchyłek w paśmie przepustowym i zaporowym. Skończona długość odpowiedzi impulsowej, a co za tym idzie skończona liczba współczynników filtru ogranicza dokładność aproksymacji. Małe zmiany amplitudy (zafalowania) w paśmie przejściowym, duże tłumienie w paśmie zaporowym i strome przejście pomiędzy nimi są interesującymi nas parametrami filtru.

Istnieje pewna liczba możliwych kryteriów oceny optymalnego wyboru współczynników, aby leżały możliwie blisko aproksymowanej charakterystyki. Wśród nich można wyróżnić cztery najważniejsze:

1. Współczynniki filtru są dobrane w ten sposób, aby minimalizować średniokwadratową odchyłkę pomiędzy założoną, a aktualnie otrzymaną charakterystyką.

2. Współczynniki filtru są tak dobrane, aby minimalizować maksymalna różnicę pomiędzy charakterystyką pożądaną, a projektowaną.

3. Współczynniki są tak wyznaczane, żeby charakterystyka zadana i aktualnie otrzymana zgadzały się w kilku punktach z zachowaniem wzmocnienia i liczby pochodnych.

n współczynników filtru jest tak ustalonych, by pożądana i projektowana charakterystyka pokrywały się w n punktach. W tym przypadku należy rozwiązać układ n równań.

2.2. Projektowanie filtrów FIR metodą próbkowania w dziedzinie częstotliwości

Metoda ta polega na otrzymaniu współczynników filtru ze zbioru n próbek charakterystyki częstotliwościowej. Otrzymana charakterystyka pokrywa się z zadaną w punktach, określonych przez położenie próbek amplitudy. Istnieją trzy sposoby podejścia do tego problemu:

1. użycie wprost odwrotnej dyskretnej transformaty Fouriera dla ciągu równo rozłożonych próbek charakterystyki częstotliwościowej,

2. ułożenie wzoru na podstawie własności filtru FIR o liniowej fazie i odwrotnej transformaty Fouriera, dzięki któremu możemy obliczyć wartości poszczególnych współczynników filtru FIR. To podejście także wymaga równo rozłożonych próbek charakterystyki częstotliwościowej,

3. ułożenie i rozwiązanie układu równań. Ten sposób nie wymaga próbek równo oddalonych od siebie.

W pracy zostaną omówione pewne własności tej metody projektowania, w oparciu na podejściu z punktu 2. W tej metodzie projektowania filtrów cyfrowych, na wstępie musimy podać wartości charakterystyki częstotliwościowej w równo od siebie oddalonych częstotliwościach:

(2.7)

Na tak zdefiniowanym zbiorze próbek charakterystyki częstotliwościowej dokonuje się odwrotnego dyskretnego przekształcenia Fouriera w wyniku otrzymujemy współczynniki filtru. Można też korzystając z właściwości filtru o rzeczywistych współczynnikach i liniowej fazie napisać wzór określający odpowiedź impulsową filtru:

(2.8)

gdzie A_k są rzeczywiste, a [M/2] oznacza całkowitą część M/2. Wzór (2.8) można rozpisać do następującej postaci:

(2.9)

gdzie:

(2.10)

Po dokonaniu transformaty z wyrażenia określonego wzorem (2.10):

(2.11)

Z powyższego wzoru możemy wyznaczyć kształt pojedynczego elementu składowego charakterystyki częstotliwościowej filtru:

(2.12)

Kształt pojedynczej próbki charakterystyki częstotliwościowej jest przedstawiony na rysunku 2.3. Charakterystyka amplitudowa pojedynczego elementu ma maksymalną wartość równą N|A_k| przy częstotliwości równej ω_k=kω_s/N. W idealnym przypadku taki element przenosiłby tylko jedną częstotliwość ω_k. Widać, że rozważany element przenosi, choć stłumione, prawie całe pasmo częstotliwości.

Rys.2.3. Kształt pojedynczej próbki w dziedzinie częstotliwości.

Wiedząc, że charakterystyka filtru jest sumą charakterystyk poszczególnych elementów można spodziewać się w otrzymanej dużych zafalowań w paśmie przejściowym oraz małe tłumienie w paśmie zaporowym. Te niekorzystne cechy nasilają się zwłaszcza przy granicy pomiędzy pasmem zaporowym, a pasmem przepustowym. Istnieją dwie metody usunięcia tej niekorzystnej cechy. Pierwszą z nich jest użycie w spróbowanej charakterystyce pasma przejściowego, które wstawiamy w miejsce nieciągłości w charakterystyce filtru. Przy takim podejściu mamy dwie drogi możliwe do zastosowania. Pierwsza z nich jest użycie funkcji, która będzie opisywała zachowanie się funkcji w paśmie przejściowym. Będę ją nazywał funkcją przejścia (ang. transition function). Najprostszą funkcją przejścia jest funkcja liniowa. Takie zastosowanie funkcji przejścia umożliwia zastosowanie każdej z możliwych implementacji metody próbkowania charakterystyki częstotliwościowej. Drugim zaś podejściem jest zastosowanie nierównomiernego próbkowania, ale możemy go używać tylko wtedy, kiedy dla uzyskania współczynników napiszemy odpowiedni układ równań (metoda 3). Zasada poprawy charakterystyki jest prosta: nie umieszczać próbek charakterystyki paśmie przejściowym. W ten sposób nie nakładamy na pasmo przejściowe żadnych ograniczeń i dlatego w tym paśmie mogą pojawić się pewne niespodziewane zniekształcenia charakterystyki. Takie pasmo przejściowe możemy nazwać pasmem do zaniedbania (ang. don't care). Drugą możliwością polepszenia właściwości filtru projektowanego omawianą metodą jest użycie funkcji okna. Otrzymane w standardowy sposób współczynniki filtru są przemnażane przez wartości okna. Dzięki temu zmniejszane zafalowania w paśmie przepustowym, a tłumienie w paśmie zaporowym rośnie. Okna są szerzej opisane w rozdziale dotyczącym projektowania filtrów FIR metodą okien. Obie metody poprawiania charakterystyki maja tą samą wadę: ich użycie zmniejsza stromość przejścia pomiędzy pasmem przepustowym, a pasmem zaporowym.

2.3. Projektowanie filtrów FIR metodą okien

Najprostszą i najbardziej oczywistą metodą uzyskania współczynników filtrów FIR jest obcięcie przesuniętej odpowiedzi impulsowej idealnego filtru. Idealny filtr dolnoprzepustowy dany posiada następującą charakterystykę częstotliwościową:

(2.13)

gdzie ω₀ jest częstotliwością graniczną pomiędzy pasmem przepustowym, a pasmem zaporowym. Odwrotna dyskretna transformata Fouriera A_d(ω) wynosi:

(2.14)

Uzyskana odpowiedź impulsowa h_d jest nieprzyczynowa i nieskończona. Chcemy otrzymać filtr rzędu M o następującej odpowiedzi impulsowej:

(2.15)

Aby filtr był przyczynowy musimy przesunąć go w prawo o liczbę L = (M-1)/2 i dokonać obcięcia mnożąc przesuniętą nieskończoną odpowiedź impulsową przez funkcję okna prostokątnego, która dane jest wzorem:

(2.16)

W wyniku przesunięcia i przemnożenia przez okno otrzymujemy następujący wzór określający współczynniki filtru:

(2.17)

Jeśli M jest liczbą nieparzystą, to wartość h(n) dla n równego L wynosi:

(2.18)

Jest dobrze znanym faktem, że minimalny błąd średniokwadratowy przy aproksymacji funkcji okresowej przez skończony szereg Fouriera jest otrzymywany poprzez obcięcie nieskończonego szeregu Fouriera dla tej funkcji. Dlatego projekt zachowujący wyżej podany warunek, jest w rzeczywistości, metodą minimalizującą błąd średniokwadratowy aproksymacji filtru idealnego. Jednak obcięcie szeregu Fouriera powoduje pojawienie się zjawiska Gibbs'a w otrzymanym filtrze, zwłaszcza że charakterystyka filtru idealnego jest nieciągła. Z uwagi na fakt, że wszystkie idealne filtry, które mają za zadanie selekcje pewnego pasma częstotliwości, są nieciągłe w krawędziach pasm, proste obcięcie często prowadzi do niezadowalających wyników, to znaczy za małego tłumienia w paśmie zaporowym oraz stosunkowo dużych zafalowań charakterystyki w paśmie przepustowym.

Aby pokazać bardziej precyzyjnie efekt wymnożenia nieskończonej odpowiedzi impulsowej filtru przez okno prostokątne musimy przejść do dziedziny częstotliwości, w której mnożenie zamienia się w splot idealnej charakterystyki filtru z transformatą Fouriera okna prostokątnego:

(2.19)

Charakterystyka amplitudowa |W_r(ω)| dana jest wzorem:

(2.20)

Jest ona przedstawiona na poniższym rysunku:

Rys 2.4. Kształt charakterystyki amplitudowej okna prostokątnego.

Jak widać na ilustracji, okno posiada listki boczne o stosunkowo dużej wartości. Uwzględniając fakt, że A_d(ω) składa się jedynie z jedności i zer, mnożenie splotowe upraszcza się do sumowania poprzesuwanych charakterystyk W_r(ω). Powoduje to powstanie wyraźnych odchyłek od charakterystyki idealnej, zarówno w paśmie przepustowym, jak i w zaporowym. Zafalowania te mają wartość wynoszącą około 9 procent. Daje to maksymalną odchyłkę w paśmie przepustowym równą 0.75dB i minimalne tłumienie około 21dB. Ponadto, wraz ze wzrostem rzędu okna zmniejsza się szerokość prążka głównego oraz maleją wszystkie oprócz pierwszego prążki boczne. Dlatego zafalowania zawężają się wokół granicy pasm, tłumienie w paśmie zaporowym rośnie oraz przejście pomiędzy pasmami staje się bardziej strome. Niestety maksymalna wartość odchyłki, pomimo wzrostu długości filtru pozostaje prawie taka sama. Jak widać, powstanie zafalowań związane jest z kształtem okna obcinającego nieskończoną odpowiedź impulsową filtru idealnego, a ściślej mówiąc ze stromością zboczy okna. Okazuje się, że można uzyskać mniejsze zafalowania poprzez stosowanie funkcji okien o mniej stromych zboczach. Jednak ceną za zminimalizowanie zjawiska Gibbs'a jest powiększenie obszaru przejściowego pomiędzy pasmami. Właściwości okien związane są z ich charakterystykami częstotliwościowymi. Generalnie nieprostokątne okna mają szerszy prążek główny, natomiast ich prążki boczne są znacznie mniejsze niż okna prostokątnego. Tablica (2.1) podaje najważniejsze parametry najpopularniejszych okien.

Tablica 2.1, Najważniejsze parametry najpopularniejszych okien

Typ okna	Poziom największego prążka bocznego [dB]	Przybliżona szerokość pasma przejściowego	Minimalne tłumienie w paśmie zaporowym [dB]
Prostokątne	-13	4π/M	21
Hanninga	-31	8π/M	44
Hamminga	-41	8π/M	53
Blackmana	-57	12π/M	74

Okno Hanninga ma wszystkie prążki boczne mniej więcej tego samego poziomu. Okno Hamming'a ma najmniejsze listki boczne dla okien o szerokości pasma przejściowego 8π/M. Natomiast okno Blackman'a ma największe tłumienie, że wszystkich stałych okien, ale ceną za to jest relatywnie szerokie pasmo przejściowe, które jest 3 razy większe od pasma przejściowego okna prostokątnego i 1.5 razy większe od szerokości pasma dla okna Hamminga. Poniższy rysunek pozwala unaocznić powyżej omówione właściwości okien:

Rys. 2.5. Charakterystyki amplitudowej różnych okien w skali logarytmicznej.

Bardzo elastyczną rodziną okien są okna o zmiennych parametrach. Do tej grupy należy okna Kaiser'a dane następującym wzorem:

(2.21)

gdzie I₀ jest zmodyfikowaną funkcją Bessela pierwszego rodzaju, a β jest parametrem determinującym stosunek pomiędzy szerokością głównego listka, a największym prążkiem bocznym. Innymi słowy określa relacje pomiędzy szerokością pasma przejściowego, a tłumieniem w paśmie zaporowym. Te okna są niemal optymalne pod względem posiadania największej energii w listku głównym w stosunku do największego prążka bocznego. Oprócz tych dobrych właściwości okno Kaiser'a pozwala na oszacowanie rzędu filtru o danej szerokości pasma przejściowego i tłumieniu w paśmie przepustowym. Co więcej pozwala wyznaczyć współczynnik β, dzięki któremu zaprojektowany filtr spełni narzucone na niego wymagania. Estymacja rzędu M wymaga podania dwóch parametrów. Pierwszym z nich jest znormalizowana szerokość pasma przejściowego:

(2.22)

gdzie ω_r jest częstotliwością graniczną pasma zaporowego, ω_p częstotliwością graniczną pasma przepustowego, a ω_s częstotliwością próbkowania. Drugim parametrem jest tłumienie w paśmie przepustowym dane wzorem:

(2.23)

gdzie δ₂ jest poziomem najwyższego prążka w paśmie zaporowym. Rząd filtru może zostać oszacowany za pomocą wzoru:

(2.24)

Natomiast współczynnik β może być otrzymany ze wzoru:

(2.25)

W poniższej tabeli zawarto zmiany parametrów okna Kaiser'a, jak i odpowiednich własności filtru dla różnych współczynników β:

Tablica 2.2.Zmiany parametrów filtru zaprojektowanego przy pomocy okna Kaiser'a dla różnych β

Parametr β	Poziom największego prążka bocznego [dB]	Szerokość pasma przejściowego	Minimalne tłumienie w paśmie zaporowym [dB]
2,0	-19	3π/M	29
3,0	-24	4π/M	37
4,0	-30	5,2π/M.	45
5,0	-37	6,4π/M.	54
6,0	-44	7,6π/M.	63
7,0	-51	9π/M.	72
8,0	-59	10,1π/M.	81
9,0	-67	11,4π/M.	90
10,0	-74	12,8π/M.	99

Na poniższym rysunku można porównać charakterystyki filtrów FIR otrzymane przy pomocy okna Kaiser'a przy różnych wartościach współczynnika β:

Rys 2.6. Filtry FIR zaprojektowane przy pomocy okna Kaiser'a dla różnych współczynników β.

2.4. Projektowanie filtrów przy pomocy minimalizacji błędu średniokwadratowego

Punktem wyjścia dla tej metody jest wymaganie, aby średniokwadratowa odchyłka pomiędzy otrzymaną w wyniku charakterystyką, a charakterystyką idealną, była minimalna. Przy rozważaniach jest przyjmowane założenie, że charakterystyka idealna H_D(ω) posiada fazę równą zeru, co usuwa wszelkie założenia dotyczące charakterystyki fazowej z naszego projektu. Błąd średniokwadratowy E_LS może być wyrażony wzorem:

(2.26)

Zakładając symetrie odpowiedzi impulsowej (czyli liniowość) oraz jej parzystość lub nieparzystość, możemy w miejsce H₀(ω) podstawić odpowiednie równanie opisujące współczynniki tego filtru. Aby uzyskać minimalizację błędu E_LS należy najpierw zróżniczkować równanie opisujące błąd względem kolejnych współczynników filtru, a następnie przyrównać powstałe równania do zera. Przedstawię ten proces na przykładzie filtru typu I. Błąd średniokwadratowy dla tego filtru można wyrazić wzorem:

(2.27)

Różniczkując powyższe równanie względem k-tego współczynnika filtru i przyrównując wynik do zera otrzymujemy:

(2.28)

Na lewej stronie powyższego równania otrzymujemy niezerowy składnik sumy tylko dla n = k. Wszystkie pozostałe składniki sumy znikają, ponieważ całkujemy na przedziale będącym całkowitą wielokrotnością okresu kosinusów zawartych pod całką.

(2.29)

Z tego wzoru można wyznaczyć poszczególne współczynniki filtru:

(2.30)

Podobnie można wyznaczyć wzory na współczynniki dla innych typów filtru. Przy takim podejściu błąd średniokwadratowy jest minimalizowany w całym przedziale częstotliwości ω = [0,π]. Jednak dla filtru fizycznie realizowanego, przebieg amplitudy charakterystyki częstotliwościowej filtru w paśmie przejściowym musi mieć skończone nachylenie. Dlatego w tej metodzie wyłącza się obszar przejściowy z przedziału, na którym dokonujemy minimalizacji błędu średniokwadratowego. Wyraża to wzór:

(2.31)

Tak więc, optymalizujemy kształt charakterystyki w pasmach przejściowych i zaporowych przy pomocy kryterium błędu średniokwadratowego, ale nie przykładamy żadnej uwagi do wymuszenia przebiegu charakterystyki w pasmach przejściowych, które pozostają nieokreślone. Teraz nie całkujemy już funkcji błędu po całym przedziale częstotliwości, lecz dokonujemy tej operacji na określonych pasmach. Poniższy wzór opisuje przypadek filtru dolnoprzepustowego:

(2.32)

Powyższy wzór może zostać łatwo zapisany w postaci macierzowej:

C * b = h

(2.33)

Elementy macierzy C mogą zostać obliczone ze wzoru:

(2.34)

Macierz C jest niezależną funkcją przejścia, którą należy aproksymować. Zawiera ona tylko informacje o przedziałach częstotliwości, w których ma zostać dokonana minimalizacja błędu średniokwadratowego. Wektor h zawiera informacje na temat kształtu pożądanej charakterystyki amplitudowej i jest określony następującym wyrażeniem:

(2.35)

Całki w powyższych wyrażenia są elementarne. W literaturze można znaleźć odpowiednie wzory determinujące elementy macierzy C i wektora h dla różnych typów filtrów. Aby wyznaczyć współczynniki filtru wystarczy rozwiązać równanie macierzowe. W prezentowanej metodzie, która minimalizuje średniokwadratową odchyłkę pomiędzy założoną, a otrzymaną charakterystyką filtru, istnieją pewne stopnie swobody w projektowaniu filtru. Poza jego rzędem i częstotliwością odcięcia, które możemy arbitralnie wybrać, możemy regulować zafalowania w paśmie przejściowym i tłumienie w paśmie zaporowym poprzez wybór szerokości pasma przejściowego. Zarówno zafalowania, jak i tłumienie nie może zostać wybrane niezależnie od siebie. W omawianym przypadku maksymalny błąd średniokwadratowy δ_max w całym interesującym nas przedziale jest stały.

Zafalowania w paśmie przejściowym D_P i tłumienie w paśmie zaporowym D_S można podać następującymi wzorami wyrażającymi te parametry filtru w skali liniowej:

(2.36)

Widać, zatem oba parametry zależą do stopnia minimalizacji błędu średniokwadratowego i nie można rozpatrywać ich oddzielnie. Przy projektowaniu filtru byłoby cecha pożądaną możliwość poprawy parametrów jednego pasma kosztem drugiego, bez zmiany długości filtru.

Poprzez wprowadzenie funkcji wagi do wzoru określającego miarę błędu otrzymujemy taką możliwość.

(2.37)

W(ω) określa jaki jest stopień redukcji błędu średniokwadratowego w określonym paśmie częstotliwości. Im dane pasmo posiada wyższą wagę, tym błąd w tym paśmie jest mniejszy, a co za tym idzie tłumienie lub zafalowania w tym paśmie są mniejsze. W tablicy pokazano wpływ wag dla filtru rzędu 34:

Tablica 2.3. Wpływ funkcji wagi na parametry filtru.

Waga pasma przepustowego	Waga pasma zaporowego	Maksymalne zafalowania w paśmie przepustowym [dB]	Minimalne tłumienie w paśmie zaporowym [dB]
1	1	0.13	38
100	1	0.03	28
1	100	0.38	48

Nie istnieje matematyczna formuła pozwalająca na estymację rzędu filtru, szerokości pasma przepustowego i funkcji wagi, aby otrzymać filtr o określonych parametrach. Aby otrzymać filtr o określonej charakterystyce należy dokonywać ręcznej zmiany parametrów.

2.5. Projektowanie filtrów FIR przy pomocy minimalizacji błędu maksymalnego

W tej metodzie projektowania kryterium otrzymania filtru optymalnego jest minimalizacja maksymalnej odchyłki filtru zaprojektowanego od filtru zadanego. Minimalizacja błędu nie jest przeprowadzana w pasmach przejściowych. Filtr zaprojektowany tą metodą posiada zafalowania zarówno w paśmie przepustowym, jak i w paśmie zaporowym. Aby opisać procedurę projektowania załóżmy, że projektujemy filtr dolnoprzepustowy o granicy pasma przepustowego przy częstotliwości ω_p i granicy pasma zaporowego w częstotliwości ω_s. Charakterystyka filtru w paśmie przejściowym musi spełniać warunek:

(2.38)

Podobnie w paśmie zaporowym, przebieg charakterystyki musi mieścić się w granicach ±δ_s, co wyrażone jest wzorem:

(2.39)

Współczynnik δ_p. odpowiada więc zafalowaniom w paśmie przejściowym, a parametr δ_s tłumieniu w paśmie zaporowym. Jak we wstępie stwierdziliśmy charakterystykę częstotliwościową H_r każdego typu filtru FIR o liniowej fazie i rzeczywistych współczynnikach można rozpisać jako sumę funkcji trygonometrycznych. Ogólna postać charakterystyki wygląda następująco:

Hr(ω)=Q(ω) P(ω)

(2.40)

Funkcje opisujące Q(ω) i P(ω) dla wszystkich czterech typów filtrów zawarto w tablica:

Tablica 2.4. Funkcje opisujące typy filtrów FIR.

Typ filtru	Q(ω)	P(ω)
I	1
II
III
IV

Współczynniki α(k) reprezentujące charakterystykę filtru są liniowo zależne od odpowiedzi impulsowej filtru. Dodatkowo oprócz funkcji H_r(ω) zdefiniujmy rzeczywistą funkcje H_d(ω) opisującą pożądaną charakterystykę filtru i funkcję wagi błędu aproksymacji W(ω). Zadana charakterystyka częstotliwościowa w najprostszym przypadku definiuje się jako 1 dla częstotliwości leżących w paśmie przepustowym, a jako 0 dla częstotliwości w paśmie zaporowym. Funkcja wagi W(ω) pozwala na wybór stosunku błędów w różnych pasmach. Jest wygodne, aby znormalizować W(ω) do jedności w paśmie zaporowym i używać W(ω) = δ_s/δ_p. w paśmie przejściowym. Wtedy wystarczy zmieniać W(ω) w paśmie przejściowym, bo oznacza ona stosunek odchyłki w paście zaporowym do odchyłki w paśmie przepustowym.

Teraz możemy zdefiniować ważony błąd aproksymacji jako:

(2.41)

Dla wygody zdefiniujmy zmodyfikowaną funkcję wagi W^'(ω) i zmodyfikowaną funkcję pożądanej charakterystyki częstotliwościowej H^'_d jako:

(2.42)

Wtedy ważony błąd aproksymacji można zapisać jako:

(2.43)

Mając daną funkcję błędu E(ω), minimalizacja maksymalnego błędu, zwana także aproksymacją Czebyszewa, sprowadza się do znalezienia zbioru współczynników filtru b(k), który zminimalizuje maksimum wartości bezwzględnej E(ω) w pasmach dla których optymalizacja jest przeprowadzana. Matematyczna formuła opisująca ten problem jest następująca:

(2.44)

gdzie S reprezentuje zbiór pasm częstotliwości, w których ma zostać dokonana optymalizacja. W przypadku filtru dolnoprzepustowego zbiór S zawiera pasmo przepustowe i pasmo zaporowe projektowanego filtru. Rozwiązanie tego problemu podali Parks i McClellan w roku 1972. Jest ono znane pod nazwą twierdzenia o kolejnych zmian (ang. alternation theorem).

Oto jego treść: Niech S będzie zwartym podzbiorem zbioru [0.π). Wystarczającym i koniecznym warunkiem, aby

(2.45)

było jedynym, najlepszą ważoną aproksymacją funkcji H^'_d(ω) w S, jest zajście faktu, że funkcja błędu zawiera najmniej L+2 ekstremów w S. Oznacza to, że musi istnieć L+2 częstotliwości ω_i w zbiorze S, takich że, ω₁<ω₂<...<ω_L+2 dla których E(ω_i) = -E(ω_i+1), a ponadto

(2.46)

Zauważmy, że funkcja błędu zmienia znak pomiędzy dwoma kolejnymi ekstremami.
Użycie twierdzenia o kolejnych zmianach zostanie przedstawione na przykładzie filtru dolnoprzepustowego. Wiedząc, że funkcja wagi i pożądana charakterystyka filtru są funkcjami stałymi można otrzymać:

(2.47)

Skutkiem tego, częstotliwości ω_i odpowiadające szczytom E(ω), odpowiadają szczytom, w których zaprojektowana H_r(ω) styka się z granicami tolerancji błędu.

Ponieważ H_r(ω) można przedstawić jako wielomianem rzędu L, którego wyrazy są cosinusami podniesionymi do potęgi, dla pierwszego typu filtru FIR, otrzymuje się:

(2.48)

Z tego równania można wyciągnąć wniosek, że H_r(ω) co najwięcej L - 1 ekstremów w przedziale (0,π). Dodatkowo ω = 0 i ω = π są zazwyczaj ekstremami H_r(ω), a więc i E(ω). Dlatego H_r(ω) może mieć co najwyżej L+1 ekstremów. Ponadto, częstotliwości graniczne pasm ω_p i ω_s są także ekstremami funkcji E(ω), ponieważ E(ω) posiada maksimum w tych częstotliwościach. W wyniku tych rozważań dochodzi się do wniosku, że najlepsza aproksymacja filtru FIR może posiadać co najwyżej L+3 ekstrema. Z kolei na mocy twierdzenia kolejnych zmian wiemy, że optymalny filtr w sensie Czebyszewa zawiera co najmniej L + 2 ekstrema. Dlatego filtr rozważany filtr dolnoprzepustowy może zawierać L + 2 lub L + 3 ekstremów. Generalnie filtr zawierające więcej niż L + 2 zafalowań nazywamy filtrami z dodatkowymi zafalowaniami (ang. extra ripple filters). Kiedy filtr zawiera maksymalną możliwą liczbę ekstremów nazywany jest filtrem o maksymalnych zafalowaniach (ang. maximal ripple filter).

Twierdzenie o kolejnych zmianach gwarantuje nam uzyskanie jednoznacznego rozwiązania dla problemu Czebyszewa. Dla każdej częstotliwości ekstremum ω_n, możemy napisać równanie:

(2.49)

gdzie δ jest maksymalną wartością funkcji błędu E(ω).Równanie to można przekształcić do bardziej użytecznej postaci poprzez przeniesienie na prawą stronę H^'_d(ω), następnie przeniesienie na lewą stronę równania δ i podstawienie za P(ω). W wyniku otrzymujemy równanie, dzięki któremu możemy wyznaczyć współczynniki filtru b(k) oraz maksymalny błąd δ:

(2.50)

Początkowo nie jest znany ani zbiór częstotliwości ω_n, przy których występują ekstrema, ani zbiór współczynników filtru i maksymalnej wartości błędu δ. Aby rozwiązać ten problem używamy iteracyjnego algorytmu nazwany algorytmem Remeza, w który zaczyna się od zgadnięcia zbioru ekstremalnych częstotliwości, obliczeniu P(ω) i δ oraz policzenia funkcji błędu E(ω). Z E(ω) wyliczamy kolejny zbiór L + 2 częstotliwości ekstremów i powtarzamy proces iteracyjnie, dopóki nie otrzymamy optymalnego zbioru częstotliwości ekstremów.

Bardziej efektywna procedura oblicza δ analitycznie posługując się następującą formułą:

(2.51)

gdzie

(2.52)

Wykorzystując to, że P(ω) jest trygonometrycznym wielomianem postaci:

(2.53)

oraz fakt, że ten wielomian w punktach ω_n, n=0,1,...,L+1 ma następujące wartości:

(2.54)

to istnieje możliwość użycia wzoru interpolacyjnego Lagrange'a do obliczenia P(ω):

(2.55)

gdzie x = cos(ω), x_k = cos(ω_k) i β_k dane jest wzorem:

(2.56)

Mając formułę na P(ω), możemy policzyć funkcje błędu z

(2.57)

na gęstym zbiorze punktów częstotliwości. Zwykle liczba punktów częstotliwości równa 16M, gdzie M jest długością filtru, wystarcza do uzyskania odpowiedniej dokładności. Użycie większej liczby punktów zwiększa dokładność, ale przedłuża czas działania algorytmu. Jeżeli |E(ω)| ≥ δ dla jakiś częstotliwości w zbiorze próbkującym funkcje błędu, wtedy tworzony jest nowy zbiór częstotliwości odpowiadających L+2 największym szczytom |E(ω)| i procedura powtarza się do obliczenia nowego δ. Nowy zbiór L+2 częstotliwości jest tak wybrany, aby odpowiadał ekstremom funkcji błędu |E(ω)|, algorytm wymusza wzrost parametru δ w każdej iteracji, dopóki nie zbiegnie się on do górnej granicy i wtedy otrzymany filtr będzie optymalny w sensie Czebyszewa. Inaczej mówiąc, kiedy |E(ω)| ≤ δ dla wszystkich częstotliwości w zbiorze próbkującym funkcję błędu, znaleźliśmy optymalne rozwiązanie. Kiedy tylko optymalne rozwiązanie w sensie optymalności P(ω) jest uzyskane, odpowiedź impulsowa filtru h(n) może zostać uzyskane wprost, bez obliczania współczynników α(k). Wystarczy wyznaczyć H_r(ω) dla ω =2πk/M, k = 0,1,...,(M - 1)/2 dla M nieparzystego lub M/2 dla M parzystego. Poniższy rysunek przedstawia graf przepływu algorytmu Remeza:

Rys 2.7. Algorytm Remeza.

Częstotliwości graniczne pasm filtru, które określają obszar A, w którym jest dokonywana aproksymacji, rząd filtru M, pożądany kształt charakterystyki H_d(ω) oraz wektor wag W(ω) są naturalnymi parametrami algorytmu Remeza. Liczby δ_s i δ_p., które determinują zafalowania w paśmie przejściowym i tłumienie filtru w paśmie zaporowym, są rezultatem iteracji przeprowadzanych przez algorytm. Ponieważ specyfikacje filtru podawane są zazwyczaj w postaci tolerancji odchyłek w poszczególnych pasmach, byłoby pożądane posiadanie relacji, która określałaby rząd filtru, dla którego filtr mieściłby się w ramach projektowych. Niestety nie istnieje żadna analityczna formuła pozwalająca określić związek między tymi parametrami.

Istnieje natomiast pół-empiryczny wzór, który pozwala dość precyzyjnie oszacować rząd filtru dla przypadku nieważonego W(ω) = 1 (δ = δ_s = δ_p.).

(2.58)

gdzie Δ to szerokość pasma przejściowego. W literaturze są podane wzory na oszacowanie rzędu filtru ważonego.

2.6. Projektowanie filtrów FIR metodą wymuszonej minimalizacji błędu średniokwadratowego bez określonych pasm przejściowych.

Autorzy tej metody rozważali problem definicji optymalności projektowanego filtru cyfrowego i stwierdzili, że kryterium wymuszonej minimalizacji błędu kwadratowego bez specyfikacji pasm przejściowych jest użytecznym uzupełnieniem dla istniejących kryteriów aproksymacji dla projektowania filtrów cyfrowych.

Zauważmy, przykładowo, podstawowy scenariusz projektowania filtru dolnoprzepustowego w którym sygnał użyteczny zawiera się w paśmie od 0 do ω₀ jest zakłócony szumem addytywnym, którego spektrum pokrywa cały zakres częstotliwości. W tym przypadku, bez dalszych założeń, brak pasma przejściowego naturalnie narasta z problemu usuwania szumu z sygnału nas interesującego: żadna część pasma przepustowego nie jest mniej ważna niż każda inna część tego pasma. Podobnie każda część pasma zaporowego nie jest bardziej ważna niż pozostałe. W wielu praktycznych przypadkach, nie rozdzielenia pomiędzy pasmem przepustowym, a pasmem zaporowym poprzez wstawienie pasma przejściowego pomiędzy nie. Rzeczywiście, często widmo sygnału pożądanego i zakłóceń nakładają się na siebie.

Ważnym wyjątkiem od przedstawionej powyżej sytuacji jest projektowanie filtrów, które oddzielają sygnały położone w dobrze odseparowanych pasmach. W tym przypadku i innych, w których filtrowane sygnały nie przenoszą energii w pasmach przejściowych, użycie pasma przejściowego jest dobrze uzasadnione - pasmo przejściowe nie stanowi krytycznej części projektu. Jednak nawet w tym przypadku w tym paśmie pojawia się szum lub inny nie pożądany sygnał, który chcemy usunąć. Z tego powodu anomalie zawarte w paśmie przejściowym są zjawiskiem nie pożądanym. W wielu przypadkach, pasmo przejściowe jest wprowadzane, aby usunąć lub zmniejszyć oscylacje w charakterystyce częstotliwościowej znajdujące się blisko krawędzi pasm spowodowane zjawiskiem Gibbs'a, więc pasma przejściowe nie powstają z fizyki problemu. Dla znaczących metod projektowania filtru, jest koniecznością znalezienie kryterium błędu, które nałoży na sygnał jakiś nierealistycznych wymagań, takich jak istnienie pasm oddzielających pożądany sygnał od szumu.

Załóżmy, że oznaczymy D(ω) oznacza pożądany kształt amplitudy. Przykładowo niech będzie to charakterystyka idealnego filtru dolnoprzepustowego. Aproksymacja tą nieciągłej funkcji wielomianem złożonym z kosinusów A(ω) jest podstawowym problemem projektowania filtrów. Wiele różnych metod opisanych w literaturze można wyróżnić, że względu w jaki sposób radzą sobie z tą nieciągłością. Rzeczywiście, kontrolowanie zachowania A(ω) w rejonie wokół nieciągłości ma silny wpływ na rozwój metod projektowania filtrów.

Dwie główne miary błędu aproksymacji są używane przy projektowaniu filtrów. Oznaczając więc błąd aproksymacji jako E(ω)=D(ω) - A(ω) mamy:

1. ważony, całkowy błąd kwadratowy ("błąd L₂") dany wzorem:

(2.58)

2. Ważony błąd Czebyszewa dany wzorem:

(2.59)

W obu przypadkach W{ω) jest nieujemną funkcją ważenia błędu. Kiedy W(ω) jest ustawiona jako jedynka na całym przedziale [0,π], miara błędu aproksymacji są nazywane nie ważonym całkowym błędem kwadratowym i nie ważonym błędem Czebyszewa. Najprostsza metoda projektowania optymalnych filtrów FIR minimalizuje ||A(ω) - D(ω)||². W wyniku otrzymujemy filtr, który nazywamy najlepszym filtrem L₂. Jak dobrze wiemy, jeśli funkcja wagi jest ustawiona na jeden na całym przedziale [0, π], wtedy najlepszy filtr w sensie L₂ jest otrzymywany poprzez obcięcie odwrotnej transformacji Fouriera D(ω) (metoda okien z użyciem okna prostokątnego). Wtedy, dla prosto opisanej funkcji D(ω) wyrażenie na współczynniki filtru może zostać łatwo znalezione. Jednak, W(ω) = 1 nie jest generalnie używane w praktyce, ponieważ najlepszy filtr w sensie L₂ będzie posiadał duże ekstrema błędu w pobliżu granic pasm. Co więcej szczyty wartości tych zafalowań nie będą się zmniejszać wraz ze wzrostem rzędu filtru.

Aby pokonać tę niedogodność, znaną jako zjawisko Gibbs'a, trzy główne podejścia mogą zostać użyte:

1. użycie okna różnego od prostokątnego,

2. użycie funkcji transmitancji z ciągłym przejściem pomiędzy pasmami,

3. użycie pasm przejściowych umieszczonych pomiędzy pasmami, w których waga błędu wynosi zero.

Użycie tych pomysłów zostało zastosowane w wielu metodach projektowania filtrów posiadających dwie pożądane właściwości:

1. tworzą filtr, na który nie wpływa zjawisko Gibbs'a,

2. procedura projektowa jest zaimplementowana przy użyciu obliczeniowo efektywnego algorytmu,

Tak się nieszczęśliwie składa, że każda z wyżej wymienionych metod posiada własne ograniczenia:

1. Okna. Mimo, że mnożenie współczynników uzyskanych z transformaty Fouriera jest proste, otrzymany filtr nie jest rozwiązaniem optymalnym, ponieważ ma zbyt wielką odchyłkę od idealnej charakterystyki.

2. Funkcja przejścia. Dzięki zmodyfikowaniu pożądanej funkcji amplitudy (tak, aby nie posiadała nieciągłości ) daje aproksymację, która nie jest zniekształcona efektem Gibbs'a, to metoda ta nie potrafi dobrze oddać kształtu zadanej charakterystyki.

3. Pasma przejściowe o wadze równej zero. Poprzez to pojęcie rozumiemy, że pomiędzy dwa pasma wstawiamy obszar, którego błąd jest zawsze równy zeru. Takie rejony w terminologii angielskiej nazywa się obszarami "don't care". Użycie takich obszarów przejściowych ułatwia rozwiązanie problemu aproksymacji. Jednakże, jeśli ich używamy, a ponadto jeśli sygnał wejściowy posiada energię w paśmie przejściowym, to optymalność otrzymanego filtru jest problematyczna.

W świetle powyższej dyskusji autorzy sugerują, że użycie dokładnie wyznaczonych pasm przejściowych w projektowaniu filtrów FIR związane jest z potrzebą redukcji szczytów błędu oraz, w pewnych przypadkach, dla ułatwienia obliczeń. W związku z tym faktem, stwierdzają, że ich istnienie nie wynika z istoty problemu. Co więcej, ponieważ minimalizacja normy Czebyszewa dla filtru dolnoprzepustowego wymaga użycia pasma przejściowego, twórcy algorytmu zaproponowali normę L₂ jako wskaźnik optymalizacji i fakt, że ekstrema błędu będą kontrolowane poprzez narzucone ograniczenie. Poprzez wyrażenie ekstrema błędu należy rozumieć lokalne minima i maksima A(ω). Przyjmijmy, że częstotliwość odcięcia oznaczymy jako ω_o. Wymagania filtru narzucone na filtr, który chcemy zaprojektować używając ważonej normy L₂, są dane poniższymi równaniami:

1. W(ω) = 1 dla każdego ω ∈ [0, π].

2. U(ω) = 1 + δ_p, L(ω) = 1 - δ_p. dla każdego ω ∈[0, ω₀].

3. U(ω) = δ_s, L(ω) = - δ_s. dla każdego ω ∈[ω₀, π].

Minimalizacja nie ważonego całkowego błędu kwadratowego polega na wymuszeniu, żeby lokalne minima i maksima funkcji A(ω) leżało wewnątrz górnej i dolnej funkcji ograniczającej L(ω) i U(ω). Najważniejsze trzy cechy tego podejścia do problemu projektowania filtru:

1. Nie ma rejonu wokół nieciągłości, który byłby wyłączony z miary błędu aproksymacji. Algorytm minimalizuje błąd kwadratowy w całym zakresie częstotliwości ze szczególnym uwzględnieniem ograniczenia maksymalnych wartości błędu.

2. Rejon przejściowy jest wprost zdefiniowany poprzez procedurę minimalizacji normy L₂

3. Nie ma minimalnego osiągalnego błędu. Oznacza to, że można projektować filtry o dowolnie małych zafalowaniach (oczywiście rząd takiego filtru będzie duży)

Algorytm minimalizuje błąd średniokwadratowy amplitudy A(ω), tak aby lokalne ekstrema funkcji A(ω) dotykały górnej i dolnej funkcji ograniczającej w pewnej liczbie częstotliwości. Jeżeli te częstotliwości są znane wcześniej, wtedy współczynniki filtru mogłyby być znalezione poprzez minimalizowanie |E(ω)|₂², w taki sposób aby spełnienie narzuconych wymagań na charakterystykę częstotliwościową w częstotliwościach ekstremów. Procedura opisana poniżej dobiera właściwie zbiór częstotliwości poprzez rozwiązanie liniowego układu równań. Ograniczenia są nałożone na wartości A(ω_i), gdzie ω_i są zawarte w tak zwanym zbiorze ograniczeń. W czasie każdej iteracji zbiór ograniczeń jest zmieniany, tak że, przy osiągnięciu przez algorytm zbieżności, zbiór ograniczeń zawierał tylko częstotliwości dla których charakterystyka dotyka narzuconych na dane pasmo ograniczeń. Problem równego ograniczenia jest rozwiązywany za pomocą mnożników Lagrange'a. Przy pomocy warunków Kuhna-Tuckera, mając dany równo ograniczony problem można rozwiązać odpowiadający mu problem nierówno ograniczony, jeżeli wszystkie mnożniki Lagrange są nieujemne. Jeśli w dowolnej iteracji zdarzy się, że mnożnik będzie ujemny, wtedy rozwiązanie równo ograniczonego problemu nie rozwiązuje odpowiadającego mu problemu nierówno ograniczonego. Z tego powodu algorytm przed zmianą zbioru ograniczeń, wyklucza z tego zbioru częstotliwości dla których mnożniki Lagrange są ujemne. W ten sposób zmniejszana jest czas wykonywania procedury projektowej i wymusza jej zbieżność.

Załóżmy, że zbiór ograniczeń S jest zbiorem częstotliwości S = {ω₁, ...,ω_r}, gdzie każda częstotliwość leży wewnątrz przedziału [0,π]. Załóżmy, że S jest podzielone na dwa podzbiory S_l i S_u, gdzie S_l jest zbiorem częstotliwości, na których chcemy wymusić ograniczenie:

A(ω)=L(ω)

natomiast S_u jest zbiorem częstotliwości, na którym chcemy wymusić ograniczenie:

A(ω)=U(ω).

Przyjmijmy, że S_l = {ω₁,...,ω_q} i że S_u = {ω_q+1,...,ω_r}. Aby zminimalizować błąd |E(ω)|₂ tworzymy lagrangian L:

(2.60)

Koniecznym warunkiem, aby zminimalizować ||E(ω)||₂, jest przyrównanie pochodnych funkcji L względem współczynników filtru a_k i mnożników Lagrange'a μ_i do zera. Prowadzi to do poniższych równań:

(2.61)

Zgodnie z warunkami Kuhna-Tuckera, kiedy wszystkie mnożniki Lagrange'a są nieujemne minimalizuje ||E(ω)||₂ wymuszając ograniczenie nierówności:

(2.62)

Korzystając ze znajomości wzorów określających charakterystykę amplitudową filtru FIR o liniowej fazie możemy układ równań przekształcić do postaci:

a + G^tμ = c Ga = d

(2.63)

gdzie

• A wektor o długości M+1 zawierający nieznane współczynniki projektowanego filtru,

• μ wektor zawierający r nieznanych mnożników Lagrange'a,

• G macierz o rozmiarach r na M+1 o elementach równych wartości charakterystyki amplitudowej filtru w częstotliwościach zawartych w zbiorze ograniczeń,

• c wektor, którego poszczególne wyrazy są współczynnikami wolnego od ograniczeń filtru optymalnego w sensie błędu średniokwadratowego o tej samej charakterystyce amplitudowej. W implementacji tej metody jest to filtr zaprojektowany metodą okien,

• d ten wektor zawiera ograniczenia nałożone na charakterystykę.

Po dokonaniu łatwych przekształceń można uzyskać wzory na współczynniki filtru a i mnożniki Lagrange'a μ:

μ = (GG^t)^-1(Gc - d) a = c - G^tμ

(2.64)

Jeśli liczba ograniczeń r jest mała w stosunku do rzędu filtru M, to wtedy układ równań jest obliczeniowo bardzo oszczędny. Interesującą sprawą jest to, ze współczynniki kosinusów a są obliczane przez dodanie wyrażenia korygującego do współczynników Fouriera. Widać różnicę w stosunku do metody okien, gdzie współczynniki optymalnego filtru są mnożone przez funkcje okna.

Algorytm rozpoczyna się z pustym zbiorem ograniczeń, więc pierwszy filtr wynikowy będzie najlepszym filtrem uzyskanym metodą minimalizacji błędu średniokwadratowego bez wprowadzenia ograniczeń. Wtedy ograniczenia są stopniowo wprowadzane do A(ω) w wybranych częstotliwościach, dopóki nie zostanie otrzymany optymalny filtr. Zbiór ograniczeń jest zmieniany po pierwsze przez znajdowanie lokalnych maksimów A(ω), które przekraczają górne ograniczenie U(ω), a po drugie przez znajdowanie lokalnych minimów A(ω), które przekraczają dolne ograniczenie L(ω). Tworzony jest wtedy zbiór ograniczeń, który tworzą częstotliwości ekstremów przekraczających nałożone na charakterystykę ramy projektowe. Następnie ponownie obliczone są mnożniki Lagrange'a dla częstotliwości zawartych w zbiorze ograniczeń. Jeżeli którykolwiek z mnożników Lagrange'a jest ujemny, to usuwamy ze zbioru ograniczeń częstotliwość odpowiadającą mnożnikowi o najbardziej ujemnej wartości i ponownie liczymy mnożniki Lagrange'a. Jeśli nie ma ujemnych współczynników obliczamy współczynniki filtru, zmieniamy zbiór S znajdując ekstrema funkcji A(ω)i sprawdzamy zbieżność algorytmu. Zbieżność jest osiągana, gdy różnica pomiędzy A(ω) a ograniczeniami nałożonymi na filtr jest mniejsza do pewnej arbitralnie wybranej małej liczby e. Jeśli algorytm nie osiągnął zbieżności, procedura jest powtarzana od policzenia mnożników Lagrange'a.

Rys. 2.8. Przebieg algorytmu projektowania filtru.

Algorytm został omówiony na przykładzie projektowania filtru dolnoprzepustowego, ale istnieje możliwość projektowania filtru o dowolnej charakterystyce. Ponadto jest możliwe wprowadzenie funkcji wagi, określającej stopień minimalizacji błędu w poszczególnych pasmach. Niestety nie istnieje żadna formuła pozwalająca estymować rząd filtru, który spełniłby założone przez projektanta wymagania.

2.7. Projektowanie filtrów FIR uogólnioną metodą Butterworth'a

Charakterystykę amplitudową filtru FIR o liniowej fazie o parzystym rzędzie N można opisać wzorem:

(2.65)

Można otrzymać zbiór współczynników B(n) w takiej formie, aby dawały one charakterystykę amplitudową monotonicznie opadającą w przedziale częstotliwości [0,π]. W celu obliczenia N/2 + 1 współczynników B(ω) musi ułożyć taką samą liczbę równań. Część tych równań, których liczba wynosi L, jest używana do optymalizacji charakterystyki częstotliwościowej w sąsiedztwie częstotliwości ω = 0; pozostałe równanie w liczbie K określają przebieg charakterystyki w punkcje ω = π. Parametry K, L i N powiązane są następującą zależnością:

K + L = N/2 +1

(2.66)

Jedno z równań L jest używane do znormalizowania wzmocnienia filtru w punkcie częstotliwości równym ω=0. Pozostałe L - 1 równań wymusza zerowanie się pochodnych H₀(ω) dla ω = 0.

(2.67)

Natomiast K równań ma za zadanie wymusić zerowanie się funkcji H₀(ω) i K - 1 pochodnych w punkcie ω= π:

(2.68)

Chwilowo, skupmy się na rozwiązaniu problemu za pomocą wielomianu stopnia N/2 zmiennej x:

(2.69)

Współczynniki wielomianu a_v są dobrane w ten sposób, aby P_N/2,K(0) = 1, a P_N/2,K(1) = 0. Ponadto L -1 pochodnych funkcji P_N/2,K(x) w punkcie x = 0 i K - 1 pochodnych w punkcie x=1 zanikało do zera. Rozwiązanie tego problemu dane jest wzorem:

(2.70)

.
Wyrażenie w nawiasach jest symbolem Newtona, który jest zdefiniowany jako:

(2.71)

Rys.2.9. Wykres rodziny wielomianów P.N/2,L(x) dla N=18 i L =1...9.

Rysunek 2.9 pokazuje wykresy rodziny wielomianów dla rzędu filtru N=18 w przedziale x=[0,1]. Dla tego filtru istnieje 9 różnych kształtów krzywej i 9 różnych częstotliwości odcięcia. W celu otrzymanie realizowanego filtru FIR, należy znaleźć transformację pomiędzy przestrzenią zmiennej x, a przestrzenią zmienną ω. Transformacja ta musi spełniać następujące warunki:

• transformacja musi być funkcją geometryczną,

• przedział x=[0,1] musi być odwzorowana na przedział [0,π],

• relacja pomiędzy x i ω musi być monotoniczną funkcją w przedziale [0,1], ponieważ monotoniczność charakterystyki amplitudowej musi zostać zachowana.

Wyrażenie spełniające te warunki wygląda następująco:

(2.72)

Podstawienie wzoru do wzoru prowadzi do następującego wyniku:

(2.73)

Powyższy wzór pozwala na wyznaczanie wzoru obliczającego wprost współczynniki filtru, ale uzyskana formuła jest nieelegancka. Lepszym sposobem jest obliczenie wartości P_N/2,K(ω_i) w N + 1 równo oddalonych punktach częstotliwości:

(2.74)

używając wzoru, a następnie obliczenie współczynników filtru poprzez dokonanie odwrotnej transformaty Fouriera spróbkowanej charakterystyki częstotliwościowej

Opisany proces projektowania filtru potrzebuje podania następującego zbioru parametrów: N, L i K, które w wyniku opisanych wyżej działań prowadzą do otrzymania optymalnej aproksymacji filtru o podanej specyfikacji.

Rząd filtru może być oszacowany z szerokości pasma przejściowego, przy czym jego granice sztywno ustalają częstotliwości, przy których amplituda osiąga, odpowiednio dla pasma przejściowego wartość 0.95, a dla pasma zaporowego 0.05.

(2.75)

Ponadto parametr L można dla częstotliwości odcięcia filtru f_0.5 wyznaczyć za przybliżonego wzoru:

(2.76)

Wszystkie częstotliwości w wzorach są unormowane względem połowy częstotliwości próbkowania. Praktyczna realizacja tej metody projektowania odbywa się w czterech krokach:

1. Estymacja rzędu filtru z szerokości pasma przejściowego przy użyciu wzoru (2.75),

2. Oszacowanie parametrów K i L przy użyciu wzorów. W tym kroku potrzebna jest informacja o położeniu częstotliwości odcięcia,

3. Policzenie równo rozłożonych próbek częstotliwości w N + 1 przy użyciu wyrażenia,

4. Wyznaczenie współczynników filtru poprzez odwrotną transformatę Fouriera próbek otrzymanych w poprzednim kroku.

Zaprezentowana procedura nie daje najlepszych wyników. W literaturze opisane są metody dają lepsze wyniki. Wrodzoną wadą tej metody jest duża wrażliwość współczynników na błędy kwantyzacji. Dotyczy to szczególnie filtrów używających arytmetyki stałoprzecinkowej, a już szczególnie używających krótkiego słowa do zapisywania współczynników filtru.

2.8. Porównanie metod projektowania filtrów FIR

Poniższa tabela zawiera porównanie najważniejszy cech algorytmów projektowania filtrów o liniowej fazie:

Tablica 2.5 Porównanie właściwości metod projektowania filtrów FIR.

	Próbkowania charaktery-styki	Okien	Minimalizacji błędu średniokwa-dratowego	Remeza	Z wymuszeniem ograniczeń	Uogólniona metoda Butterworth'a
Kontrola granic pasm	słaba	Słaba zależy od okna	dobra	Bardzo dobra	Bardzo dobra	Słaba
Pasmo przejściowe	Wąskie powoduje duże zafalowania dla polepszenia parametrów konieczne użycie funkcji Przejścia	Zależy od kształtu okna przy okien Kaisera możliwa jest regulacja jego szerokości	W tym paśmie nie jest dokonywana optymalizacja, więc mogą występować duże odchyłki	Algorytm optymalizuje całkowity błąd w całym przedziale ω więc i kształt tego pasma	Kształt pasma jest optymalny dla nałożonych ograniczeń	Bardzo szerokie i opadające powoli,
Tłumienie	Małe zwłaszcza w pobliży granicy pasma	Zależy od kształtu okna dla okna Kaisera możliwe jego regulacja	Dobre zmniejsza się wraz z oddalaniem od granicy pasma, lecz otrzymanie określonego tłumienia ręczne	Dobre stałe w całym paśmie zaporowym możliwe do ustalenia podczas estymacji rzędu	Algorytm wprost wymaga podania tłumienia i otrzymany filtr doskonale oddaje podaną wartość	Początkowo małe potem osiąga duże wartości
Zafalowania	Duże niekontrolowa-ne	Istnieją tylko przy użyciu okna prostokątnego	Małe rosną przy zbliżaniu do granicy, brak kontroli ich wartości	Stałe w całym paśmie, ich wartość można ustalić kosztem rzędu filtru	Algorytm zachowuje się tak samo jak w przypadku tłumienia	Nie ma zafalowań w paśmie przepusto-wym
Funkcja wagi	Nie ma	Nie ma	Jest	Jest	Jest	Nie ma
Estymacja rzędu	Nie ma	Jest tylko dla okna Kaisera	Nie ma	Jest	Nie ma	Jest
Implementacja	prosta	prosta	Średnio	Trudna	Trudna	Średnio trudna
Szacunkowy unormowany czas projektowania	1.3	1	1.94	1.95	11.5	38.5

Metoda Remeza i metoda z wymuszaniem ograniczeń ma najlepsze właściwości. Ponadto metoda Remeza daje dodatkowe możliwości, co prezentuje program gremez wchodzący w skład pakietu Matlab.

3. Projektowanie filtrów o nieskończonej odpowiedzi impulsowej

3.1. Wstęp do projektowania filtrów IIR

Przy projektowania filtrów IIR szeroko wykorzystuje się dokonania w dziedzinie projektowania filtrów analogowych. Dzięki istnieniu szeregu transformat istnieje możliwość przeniesienia charakterystyki filtru analogowego, z zachowaniem jej kształtu i stabilności filtru, na odpowiadający prototypowi filtr cyfrowy. Oprócz metody prototypu analogowego istnieją także inne:

1. Metoda projektowania z użyciem równań Yule-Walkera bazująca na analizie statystycznej i właściwościach filtru pobudzane szumem,

2. Metoda projektowania filtru w dziedzinie częstotliwości, opierająca się na poszukiwaniu minimum funkcji kryterialnej błędu danej w postaci normy rzędu p,

3. Metody wyznaczania współczynników poprzez modelowanie zadanych charakterystyk, których podstawą są najczęściej właściwości splotu. Metody te projektują filtr na podstawie wymagań opisanych w dziedzinie czasu, najczęściej w postaci odpowiedzi impulsowej. W rozwiązywaniu tego typu problemów często stosuje się kryterium błędu średniokwadratowego, jako drogę do uzyskania jak najlepszej aproksymacji.
   

3.2. Metoda prototypu analogowego

Omówiony poniżej sposób projektowania opiera się na osiągnięciach w dziedzinie projektowania filtrów analogowych. Aby uzyskać współczynniki cyfrowego filtru IIR projektant musi dokonać następujących czynności :

1. Zaprojektować prototyp analogowy,

2. Użyć odpowiedniej transformaty, aby uzyskać niezniekształcony filtr cyfrowy.

W projektowaniu filtrów IIR podstawowe znaczenie mają takie klasyczne metody projektowania filtrów analogowych jak metody: Butterworth'a, Czebyszewa i eliptyczna (zwana też metodą Cauera)

Filtry Butterworth'a posiadają charakterystykę amplitudową opisaną następującym wzorem:

(3.1)

,gdzie N jest rzędem filtru, Ω_c jest częstotliwością 3 decybelową filtru, Ω_p jest częstotliwością graniczną pasma przepustowego, natomiast 1/(1 + ε)² jest wartością charakterystyki amplitudowej dla Ω = Ω_p. Charakterystykę amplitudową opisaną wzorem (3.1) nazywamy maksymalnie płaską, ponieważ spełnia warunek:

(3.2)

Charakterystyki filtru Butterworth'a są monotoniczne i dla wartości Ω = Ω_c przybierają wartość:

(3.3)

a więc maleją w porównaniu z wartością początkową A(0) o 3 dB. Ponadto wartość pierwszej pochodnej w punkcje Ω = Ω_c wynosi:

(3.4)

Nachylenie charakterystyki amplitudowej dla pulsacji 3-decybelowej, jest stosunkowo niewielkie i liniowo rośnie wraz ze wzrostem rzędu filtru. Rysunek 3.1 przedstawia przebieg charakterystyk filtrów Butterworth'a dla różnych rzędów filtru.

Rys. 3.1. Charakterystyki filtrów Butterworth'a o różnych rzędach.

Ponieważ H(s)H(-s) obliczona dla s = jΩ jest wprost równa A²(Ω), możemy równanie (3.1) zapisać wprost w postaci:

(3.5)

Jeśli dokonany rozkładu tego stosunku wielomianów na pierwiastki to okaże się, że:

1. filtr Butterworth'a posiada jedynie bieguny,

2. bieguny są równo rozłożone na kole o promieniu Ω_C,

3. ze wzoru (3.5) można wyliczyć ich położenie:

(3.6)

4. bieguny w lewej części koła należą do H(s), zaś bieguny w prawej części koła należą do H(s).

Rząd filtru potrzebny do uzyskania tłumienia δ w paśmie zaporowym rozpoczynającym się od częstotliwości Ω_S jest dany następującym wzorem:

(3.7)

gdzie δ^' jest dane wzorem:

(3.8)

Filtr Butterworth'a jest kompletnie zdefiniowany podanie parametrów: rzędu filtru N, tłumienia δ, współczynnika ε oraz stosunku granic pasm Ω_S/Ω_P.

Filtry Czebyszewa dzielą się na dwie grupy. Filtry typu I są filtrami posiadającymi tylko bieguny, które kształtują równomierne zafalowania filtru w paśmie przepustowym i monotoniczne opadanie charakterystyki w paśmie zaporowym. Natomiast filtry Czebyszewa drugiego rodzaju posiadają równą liczbę zer i biegunów, tak położonych w przestrzeni zmiennej s, aby otrzymać monotoniczne zachowanie w części przepustowej i równo zafalowaną charakterystykę w części zaporowej. Zera tej klasy filtrów leżą na osi urojonej przestrzeni s. Charakterystyka amplitudowa filtru Czebyszewa typu I jest dana wzorem:

(3.9)

Gdzie parametr ε odnosi się do wielkości zafalowań w paśmie przejściowym, a T_N(x) jest wielomianem Czebyszewa stopnia N zdefiniowanym następująco:

(3.10)

Wielomiany Czebyszewa mogą także być tworzone poprzez następujący wzór rekurencyjny:

TN+1(x)=2xTN(x) - TN-1(x) dla N=1,2,... T0(x) = 1 T1(x) = x

(3.11)

Ponadto wielomiany te spełniają następujące warunki:

1. wartość bezwzględna T_N(x) jest mniejsza lub równa jedności dla x z przedziału [-1,1],

2. T_N(1) = 1 dla wszystkich N,

3. wszystkie pierwiastki wielomianu T_N(x) mieszczą się w przedziale [-1,1],

4. w przedziale (-1,1) wielomian T_N(x) ma N -1 ekstremów o wartościach ±1,

5. spośród wszystkich wielomianów stopnia N o współczynniku a_n = 1 najmniejsze odchylenie od zera w przedziale [-1,1] ma wielomian:

(3.12)

Parametr ε jest zależny od zafalowań δ w paśmie przepustowym. Tą zależność określa wzór:

(3.13)

Bieguny filtru Czebyszewa rodzaju I leżą na elipsie w przestrzeni zmiennej s o długości półosi wielkiej danej wzorem:

(3.14)

i półosi małej:

(3.15)

gdzie β jest współczynnikiem zależnym od ε równaniem:

(3.16)

Bieguny filtru leżą na elipsie w punktach określonych przez parę współrzędnych (x_k,y_k) określonych następującą zależnością:

(3.17)

gdzie φ_k jest katem fazowym biegunów filtru Butterworth'a:

(3.18)

Filtr Czebyszewa typu II zawiera tyle samo zer, co biegunów. Charakterystyka amplitudowa dana jest wzorem:

(3.19)

gdzie T_N jest wielomianem Czebyszewa określonym wzorem (3.10), natomiast Ω_S jest częstotliwością graniczną pasma zaporowego. Zera filtru Czebyszewa rodzaju drugiego leżą na osi urojonej w punktach określonych poniższym wzorem:

(3.20)

Bieguny filtru znajdują się w punktach danych współrzędnymi (v_k,w_k),gdzie

(3.21)

Współrzędne {x_k} i {y_k} są dane wzorem (3.17) ze zmienionym współczynnikiem β, teraz odnoszącym się do tłumienia δ w paśmie zaporowym poprzez równanie:

(3.22)

Z powyższego opisu można wyciągnąć wniosek, ze każdy typ filtru Czebyszewa jest w pełni opisany poprzez parametry N,δ oraz stosunek Ω_P/Ω_S szacujący szerokość pasma przejściowego. Tak więc dla filtrów Czebyszewa typu I optymalizowane jest pasmo przepustowe, zaś dla typu II optymalizowane jest pasmo zaporowe. W literaturze dostępne są formuły na oszacowanie rzędu filtru Czebyszewa spełniającego określone przez projektanta wymagania. W ogólności filtry Czebyszewa spełniają narzucone na nie wymagania dla mniejszego rzędu filtru niż filtry Butterworth'a. Dzieje się tak, ponieważ w filtrach Chebyszewa charakterystyka w paśmie przejściowym opada znacznie gwałtowniej niż w filtrach Butterworth'a. Szybkość opadania charakterystyki w filtrach Czebyszewa zależy od kwadratu rzędu, a w filtrach Butterworth'a zależność ta jest liniowa. Na poniższym rysunku przedstawiono typowe kształty filtrów Czebyszewa:

Rys. 3.2. Charakterystyki filtrów Czebyszewa.

Filtry eliptyczne, zwane też filtrami Cauera, posiadają równomierne zafalowania, zarówno w paśmie przepustowym, jak i w zaporowym. Ta klasa filtrów zawiera równą liczbę zer i biegunów, a jej charakterystyka amplitudowa opisana jest wzorem:

(3.23)

gdzie U_N(x) jest eliptyczną funkcją Jacobiego rzędu N. Zera leżą na osi urojonej przestrzeni s. Błąd aproksymacji w filtrach eliptycznych jest równy rozciągnięty na pasmo zaporowe, jak i pasmo przepustowe. Dlatego są najbardziej efektywne z punktu widzenia osiągnięcia parametrów filtru przy najmniejszym możliwym rzędzie. Równoważnie możemy powiedzieć, że dla danego rzędu filtru i dla danego zbioru specyfikacji, filtr eliptyczny ma najwęższe pasmo przejściowe. Rząd filtru potrzebny do potrzymania filtru spełniającego następujące wymagania: tłumienie w paśmie zaporowym δ_S, zafalowania w paśmie przepustowym δ_P i stosunek częstotliwości granicznej pasma przepustowego do częstotliwości pasma zaporowego Ω_P/Ω_S szacujący szerokość pasma przejściowego, dany jest wzorem:

(3.24)

gdzie K jest całką eliptyczną pierwszego rodzaju, współczynnik δ odnosi się do tłumienia w paśmie zaporowym δ_S, a parametr ε do zafalowań w paśmie przepustowym.

Istnieje jeszcze czwarta ważna grupa filtrów analogowych: filtry Bessela. Posiadają one w paśmie przepustowym liniową fazę, zaś ich charakterystyka amplitudowa posiada bardzo szerokie pasmo przejściowe. Niestety użycie transformaty w celu uzyskania filtru cyfrowego zniekształca charakterystykę fazową, pozbawiając filtry Bessela ich waloru: liniowej fazy. Dlatego filtry Bessela nie są używane jako prototypy filtrów cyfrowych.

Krótkie porównanie właściwości analogowych prototypów filtrów cyfrowych

Najpierw porównajmy charakterystyki amplitudowe wszystkich typów prototypów dla tego samego rzędu.

Rys. 3.3. Porównanie charakterystyk amplitudowych prototypów filtrów rzędu 5 o częstotliwości granicznej 0.5.

Charakterystyka filtrów Czebyszewa typ pierwszy i filtrów eliptycznych zagina się dokładnie w częstotliwości odcięcia. W tych filtrach zafalowania w paśmie przepustowym są stałe w całym paśmie. Można powiedzieć, że filtr Czebyszewa typ I jest efektywny przy filtrach w których zależy nam dobrych parametrach filtru w paśmie przepustowym. Podobna uwaga dotyczy filtrów Czebyszewa typu II, ale w tym przypadku odnosi się do pasma zaporowego. Filtr eliptyczny ma najwęższe pasmo przejściowe, a ponadto umożliwia uzyskanie określonych parametrów zarówno w paśmie zaporowym, jak i przepustowym. Dlatego można uznać go za filtr optymalny, jeżeli chodzi o możliwość kształtowania charakterystyki amplitudowej. Z kolei filtry Butterworth'a mają najszersze pasmo przejściowe, a moment zagięcia charakterystyki jest rozmyty.

Rys. 3.4. Charakterystyki fazowe prototypów analogowych.

Charakterystyka fazowa filtrów eliptycznych jest nieliniowa, szczególnie w pobliżu częstotliwości odcięcia. Skutkiem tego częstotliwości w pobliżu częstotliwości granicznej będą przenoszone wolniej niż pozostałe częstotliwości w paśmie przepustowym. Wprowadza to duże zniekształcenia fazowe w sygnale. Dlatego w przypadkach, gdy zależy projektantowi na jak najmniejszym zniekształceniu fazy przez filtr IIR należy użyć analogowego prototypu filtru Butterworth'a, którego charakterystyka fazowa jest najbardziej zbliżona do liniowej.

Wszystkie powyższe rozważania dotyczą filtru dolnoprzepustowego. W literaturze podano odpowiednie przekształcenia pozwalające uzyskać inne rodzaje filtrów. Przekształcenia te mogą zostać dokonane zarówno na prototypie analogowym, jak i na filtrze cyfrowym. Posiadając wyliczony prototyp analogowy w dziedzinie s, należy przenieść do dziedziny z bez zniekształceń. Filtr analogowy można opisać poprzez jego transmitancję:

(3.25)

gdzie {α_k} i {β_k} są współczynnikami filtru, lub poprzez odpowiedź impulsową filtru, która jest powiązana z H₀(s) poprzez transformatę Laplace'a:

(3.26)

Równoważnie, filtr analogowy posiadający transmitancję H₀(s) daną wzorem (3.25), można korzystając z właściwości transformaty Laplace'a zapisać w postaci liniowego równania różniczkowego o stałych współczynnikach:

(3.27)

Każdy z tych trzech równoważnych zapisów prowadzi do innej metody przekształcenia filtru analogowego w jego cyfrowy odpowiednik. Przypomnijmy, że analogowy układ liniowy, stacjonarny jest stabilny, jeśli wszystkie jego bieguny leżą na lewej półpłaszczyźnie przestrzeni s. Z powyższych rozważań wynikają następujące własności, które musi posiadać transformata, aby dobrze przekształcać filtr analogowy:

1. oś urojona przestrzeni s powinna po transformacji odpowiadać okręgowi jednostkowemu w przestrzeni z. Oznacza to, że musi istnieć jednoznaczna relacja pomiędzy zmiennymi mającymi sens częstotliwości w obu dziedzinach.

2. lewa półpłaszczyzna przestrzeni s musi być odwzorowywana na wnętrze okręgu jednostkowego w przestrzeni z. Wymaganie to jest spowodowane faktem, że stabilny filtr analogowy powinien być przekształcany na stabilny filtr cyfrowy.

Jedną z prostszych metod konwertowania filtru cyfrowego na filtr cyfrowy jest transformata BD. Jej nazwa pochodzi od angielskiej nazwy tej metody (ang. backward difference) i odnosi się do numerycznej metody obliczania pochodnej. Pochodną w punkcie t = nT, możemy w przybliżeniu policzyć jako:

(3.28)

gdzie T reprezentuje okres próbkowania. Aby znaleźć związek pomiędzy zmiennymi zespolonymi s i z należy użyć przekształcenia Laplace'a do lewej strony równania (3.28), natomiast do prawej przekształcenia z. W wyniku przekształceń otrzymujemy związek pomiędzy zmiennymi:

(3.29)

Podobne rozważania dla pochodnych wyższy rzędów prowadzą do tego samego wyniku, co otrzymane powyżej. W konsekwencji transmitancja filtru cyfrowego jako wynik aproksymacji pochodnych poprzez skończone różnice daje zapisać się zależnością:

(3.30)

gdzie H_o(s) jest transmitancją filtru określonego równanie różniczkowym podanym we wzorze (3.27).Niestety transformata BD przekształca lewą półpłaszczyznę przestrzeni s we wnętrze koła o promieniu 0.5 i środku w punkcie z=0.5,co jest pokazane na rysunku 3.5.

Rys. 3.5. Przekształcenie płaszczyzny s na płaszczyznę z przy pomocy transformaty BD.

Przekształcenie przekształca stabilny prototyp w stabilny filtr cyfrowy, ale silnie zniekształca charakterystykę. Możliwe do uzyskania położenie biegunów filtru cyfrowego ogranicza przydatność tej transformaty do filtrów dolnoprzepustowych lub pasmowoprzepustowych o stosunkowo niskiej częstotliwościach odcięcia.

Inną metodą przejścia z dziedziny s na dziedzinę z jest metoda niezmienności odpowiedzi impulsowej. Odpowiedź impulsowa filtru cyfrowego składa się z próbek ciągłej odpowiedzi impulsowej filtru analogowego.

(3.31)

Omawiana transformata ma oba warunki prawidłowego odwzorowania tzn. odwzoruje oś urojoną na okręg jednostkowy i transformuje lewą półpłaszczyznę przestrzeni s we wnętrze okręgu jednostkowego. Niestety w wyniku spróbkowania odpowiedzi impulsowej, uzyskana charakterystyka filtru cyfrowego składa się z wielu poprzesuwanych o 2π/T oryginalnych charakterystyk prototypu analogowego. Wskutek nakładania się widm omawiana transformata nie nadaje się do projektowanie filtrów górnoprzepustowych oraz pasmowozaporowych. Ponadto w przypadku filtrów Czebyszewa rodzaju drugiego i filtrów eliptycznych powstałe zniekształcenia w paśmie zaporowym powodują utratę dobrych cech tych filtrów. Przed dokonaniem transformacji dokonuje się rozkładu transmitancji filtru analogowego określonego wzorem (3.25) na ułamki proste:

(3.32)

Dokonując przekształcenia Laplace'a otrzymujemy odpowiedź impulsową postaci:

(3.33)

Próbkując otrzymaną odpowiedź impulsową czyli stosując wzór (3.31) otrzymujemy odpowiedź impulsową filtru cyfrowego:

(3.34)

Następnie przechodzimy do dziedziny z stosując transformatę z:

(3.35)

Warto zauważyć, ze transformata zamienia bieguny filtru analogowego położone w s_k na bieguny filtru cyfrowego p_k położone w p_k = e^s_k^T.

Omówione dotychczas transformaty działają prawidłowo jedynie dla filtrów dolnoprzepustowych i pewnej klasy filtrów pasmowoprzepustowych. Ograniczenia tego nie posiada transformacja biliniowa. Przekształca ona oś urojoną przestrzeni s w koło jednostkowe w przestrzeni z jednoznacznie, unikając nakładania się charakterystyk. Co więcej, lewa półpłaszczyzna przestrzeni s jest przekształcana do wewnątrz koła jednostkowego, zaś prawa półpłaszczyzna jest odwzorowywana na zewnątrz koła. Transformata biliniowa jest powiązana z numeryczną aproksymacją całkowania przy pomocy metody trapezów. Dla przykładu załóżmy, że filtr analogowy ma transmitancję postaci:

(3.36)

Ten sam układ można zapisać w postaci liniowego równania różniczkowego o stałych współczynnikach:

(3.37)

Zamiast jak w metodzie BD aproksymować pochodną, dokonajmy całkowania pochodnej poprzez zastosowanie metody trapezów. Gdy rozpiszemy y(t) jako:

(3.38)

gdzie y'(t) oznacz pierwszą pochodną y(t). Aproksymując całkę we wzorze (3.38) metodą trapezów dla t = nT i t₀ = nT - 1 otrzymujemy:

(3.39)

Teraz wyznaczmy pochodną z równania (3.37) i podstawiając t = nT:

(3.40)

Podstawiając równanie (3.40) w miejsce pochodnych w wyrażeniu (3.39) i w ten sposób otrzymując równanie różnicowe dla równoważnego układu dyskretnego czasu. Dla y(n) = y(nT) i y(n-1) = y(nT - T) otrzymujemy następujący wynik:

(3.41)

Dokonując przekształcenia z równania różnicowego (3.41) i dokonując odpowiednich przekształceń otrzymujemy transmitancję filtru cyfrowego:

(3.42)

Poprzez porównanie (3.42) z (3.36) można wyznaczyć zależność pomiędzy zmienną s, a zmienną z:

(3.43)

Wadą przekształcenia biliniowego jest nieliniowe przekształcanie osi częstotliwości przestrzeni s w oś częstotliwości przestrzeni z:

(3.44)

Powoduje to, ze otrzymana charakterystyka jest zniekształcona, co objawia się poprzez zagęszczenie osi częstotliwości dla dużych częstotliwości. Aby tego uniknąć należy zadaną częstotliwość graniczną filtru cyfrowego ω wstawić do wzoru (3.44), a następnie otrzymaną wartość użyć w procedurze projektującej prototyp analogowy filtru cyfrowego.

3.3. Metoda Yule-Walkera

Metoda Yule-Walkera opiera się na relacji pomiędzy parametrami filtru, a ciągiem autokorelacji. Jeżeli widmowa gęstość mocy jest stacjonarnym procesem scholastycznym , to istnieje zależność pomiędzy ciągiem autokorelacji {γ_xx}, a współczynnikami liniowego filtru H(z) wyrażonymi następującym wzorem:

(3.45)

,przy czym filtr jest pobudzany białym szumem. Ponadto filtr opisany H(z) powinien być przyczynowy, stabilny i minimalnofazowy. Związek ten można rozpisać jako:

(3.46)

gdzie h(n) jest odpowiedzią impulsową filtru, a σ_w jest wariacją wejściowego szumu białego. Zauważmy, że dla m > q równanie (3.46) upraszcza się do postaci:

(3.47)

Powyżej podane równanie w literaturze przedmiotu jest znane pod nazwą zmodyfikowanego równania Yule-Walkera. Pozwala ono obliczyć współczynniki {a_k} znając ciąg autokorelacji. W przypadku projektowania filtru cyfrowego zazwyczaj jest dana informację o charakterystyce częstotliwościowej filtru H(z). Wiedząc, że widmowa gęstość mocy wiąże się z charakterystyką częstotliwościową za pomocą następującej zależności:

(3.48)

oraz, że ciąg autokorelacji można wyliczyć z widmowej gęstości mocy dokonując odwrotnej transformaty Fouriera:

(3.49)

Na podstawie wyrażenia (3.47) można ukłożyć p - 1 równań, wstawiając do nich odpowiednie wyrazy ciągu korelacji otrzymujemy:

(3.50)

Kolejnym krokiem będzie wyznaczenie współczynników {b_k}. W literaturze podano kilka sposobów na ich wyznaczenie. W tej pracy omówimy jedno z nich. Zauważmy, że przyczynowa część ciągu autokorelacji można zapisać jako:

(3.51)

Innymi słowy Γ_p. jest liniowym systemem, którego odpowiedź impulsowa jest jednostronnym ciągiem autokorelacji. Aby wyznaczyć n(z) używamy przybliżonej zależności:

(3.52)

gdzie π(z) jest z transformatą okna prostokątnego wycinającego przedział [1,N], zaś gwiazdka oznacza mnożenie splotowe. Równanie (3.52) przeniesione do dziedziny czasu pozwala na utworzenie układu równań z którego można wyznaczyć n(z):

(3.53)

Współczynniki h(n) są odpowiedzią impulsową filtru posiadającego tylko bieguny a_k. Wiedząc, że widmowa gęstość mocy jest funkcją parzystą, to znając część przyczynową ciągu autokorelacji możemy policzyć całkowitą gęstość widmową. Kolejnym krokiem będzie obliczenie odpowiedzi impulsowej filtru. W tym celu wykorzystamy cepstrum ciągu autokorelacji:

(3.54)

Samą gęstość widmową mocy można zapisać przy pomocy cepstrum jako:

(3.55)

Z równania 3.55 możemy wyznaczyć H(z)

(3.56)

Wykorzystując komputer do projektowania filtru wygodniej wykorzystywać zamiast transformaty z transformatę Fouriera. Wyznaczenie odpowiedzi impulsowej filtru składa się z następujących kroków:

1. obliczenie cepstrum poprzez wyznaczenie logarytmu widmowej gęstości mocy,

2. policzenie cepstrum ciągu autokorelacji jako odwrotnej transformaty Fouriera cepstrum,

3. usunięcie z cepstrum ciągu autokorelacji części nieprzyczynowej,

4. przejście do dziedziny częstotliwości poprzez transformatę Fouriera i policzenie funkcji wykładniczej o podstawie e z uzyskanego cepstrum,

5. ponowne policzenie odpowiedzi impulsowej, aby otrzymać odpowiedź impulsową filtru.

Ostatnim krokiem jest wyznaczenie współczynników {b_k} znając odpowiedź impulsową filtru h(n) i współczynniki {a_k}. Dokonujemy tego poprzez dekompozycje filtru na filtr zawierający same bieguny i filtr posiadający jedynie zera. Sygnałem wejściowym dla filtru o samych zerach będzie odpowiedz impulsowa h_a(n) filtru o współczynnikach mianownika transmitancji {a_k}, zaś jego wyjściem jest odpowiedź impulsowa całego filtru:

(3.57)

Na mocy układu równań (3.55) możemy obliczyć współczynniki {b_k}. Ten sposób otrzymania współczynników licznika transmitancji nazywamy dopasowaniem optymalnym pod względem minimalizacji błędu średniokwadratowego, ponieważ dokonujemy obcięcia nieskończonej odpowiedzi impulsowej.

3.4. Metoda projektowania filtrów IIR przy użyciu kryterium minimalizacji p-tej normy

Problem projektowania filtru IIR, która ma arbitralnie określoną charakterystykę częstotliwościową, może być rozważone jako klasyczny problem aproksymacji. Zaletą tego podejścia jest fakt, że istnieje wiele metod rozwiązania takich problemów. Klasyczny problem aproksymacji może być zadany w następujący sposób. Niech f(x) będzie daną funkcją rzeczywistą zdefiniowaną na zbiorze X, a F(A,x) jest rzeczywistą funkcją aproksymującą zależną od x∈X i od n parametrów A. Posiadając funkcję odległości p, należy dobrać n parametrów A^*∈P w taki sposób, aby:

(3.58)

Rozwiązanie tego problemu nazywamy najlepszą aproksymacją w zależności od wybranej funkcji odległości p. Dla projektowania filtrów cyfrowych użyte powyżej symbole mają swoją identyfikację w parametrach procesu projektowania. Niech transmitancja filtru H(z) będzie funkcją n zmiennych (współczynników filtru), które umieścimy w wektorze A. Niezależna zmienna x jest teraz częstotliwością ω, a zbiór X jest całym przedziałem częstotliwości X=[0,π]. Charakterystyka częstotliwościowa filtru może być wyrażona jako funkcja rzeczywista F(A,ω), podczas gdy zadana charakterystyka częstotliwościowa jest dana przez funkcję f(ω).

Mając tak zdefiniowany problem, musimy kolejno:

1. zdefiniować funkcję odległości p,

2. wybrać sposób zapisu charakterystyki częstotliwościowej H(z),

3. stworzyć wektor parametrów filtru A,

4. wybrać metodę rozwiązania problemu (3.58).

W tej metodzie jako funkcję odległości wybrano normę nazywaną w literaturze przedmiotu ważoną normą L_p.

(3.59)

gdzie w(ω) jest dodatnią funkcją wagi. Powodem wyboru tego wyboru jest fakt, że jest ona często używana, a jej właściwości są dobrze znane. Co więcej, problem znalezienia najlepszej aproksymacji w sensie normy L_p może zostać przekształcony w problem poszukiwania minimum funkcji n zmiennych, dla którego istnieją bardzo efektywne algorytmy numeryczne. Ponadto dwa najczęściej używane kryteria minimalizacji błędu, to jest minimalizacji błędu średniokwadratowego oraz minimalizacji błędu całkowitego są specjalnymi przypadkami omawianej normy dla p = 2 i dla p = ∞.

Przy implementacji numerycznej tej jest wygodne użycie prostszej funkcji odległości:

(3.60)

przy czym minimum L_2p(A) jest równe minimum ||L(A)||_2p^2p, czyli pokrywa się z minimum ||L(A)||_2p. Potęga 2p została wybrana w celu ułatwienia obliczania pochodnych. Należy też zauważyć, że punkty częstotliwości ω_k nie muszą być równo oddalone od siebie. W rzeczywistości muszą być tak dobrane, aby zmniejszana funkcja błędu okazywała sinusoidalny charakter. Liczba punktów K powinna wynosić około 10 na okres. Funkcja f(ω) opisująca wymagania filtru może zawiera informacje na temat charakterystyki amplitudowej i/lub opóźnienia grupowego. Przy projektowaniu filtru o określonej charakterystyce opóźnienia grupowego, jego charakterystyka amplitudowa ma charakter wszech przepustowy.

Forma przedstawienia transmitancji projektowanego filtru jest ważna z kilku powodów. Fizyczna struktura filtru ma duży wpływ na błędy wywołane poprzez efekty kwantyzacji. Dlatego struktury bezpośrednie, które mają najgorsze właściwości pod tym względem, nie są dobrym wyborem do przedstawienia transmitancji. Korzystniejsze właściwości mają struktury równoległa lub kaskadowa. Zwłaszcza ta ostatnia wydaje się lepszym wyborem, szczególnie jeśli chcemy projektować filtry posiadające zera na okręgu jednostkowym (tzn. dobrze zdefiniowane pasmo zaporowe). Drugim powodem ,dla którego dokonano wyboru struktury kaskadowej jest łatwość sprawdzenia stabilności filtru, podczas gdy dla struktury bezpośredniej wymaga to dłuższych obliczeń. Podczas procesu optymalizacji musimy mieć pewność, że filtr pozostawał stabilny, co znowu potwierdza przewagę struktury kaskadowej nad pozostałymi. Ponadto charakterystyka częstotliwościowa filtru przedstawiona w postaci kaskadowej składa się z prostych wyrażeń, umożliwiając łatwe policzenie pierwszych pochodnych i wprowadzając jasny wgląd na wpływ zer i biegunów na funkcję odległości. Z powyższych powodów na formę przedstawienia transmitancji filtru wybieramy kaskadę sekcji drugiego rzędu:

(3.61)

gdzie k₀ jest dodatnią stałą normalizującą wzmocnienie filtru.

Najbardziej oczywistym sposobem stworzenia wektora parametrów filtru jest ułożenie go wprost ze współczynników struktury kaskadowej:

A=[b11, b12, a11, a12, . . ., b1i, b1i, a1i, a1i, . . .,k0]

(3.62)

Jednak taki wybór sposobu ułożenia wektora A ma poważne wady. Wyprowadzone z niego wyrażenia na charakterystykę amplitudową oraz opóźnienie grupowe filtru są dość skomplikowanymi wyrażeniami zespolonymi, tak że wpływ współczynników a i b na kształt charakterystyk jest trudny do uchwycenia. Ponadto policzenie pochodnych tych wyrażeń jest dosyć trudne. Lepszym wyborem parametrów są współrzędne polarne położenie zer i biegunów:

A=[r01, φ02, rp1, φp2, . . ., r0i, φ0i, rpi, φpi, . . .,k0]

(3.63)

Przy takim wyborze wyrażenia opisujące charakterystykę częstotliwościową i opóźnienie grupowe są funkcjami rzeczywistymi, przy czym wpływ pojedynczego bieguna lub zera jest może być łatwo wyznaczony. Ponadto równania opisujące pochodne cząstkowe charakterystyki amplitudowej i opóźnienie grupowego mogą zostać w łatwy sposób wyprowadzone. Przy użyciu tych pochodnych można uzyskać pochodne funkcji odległości względem współrzędnych położenia biegunów, których obliczenie jest wymagane przez algorytm dokonujący minimalizacji funkcji błędu. Wzory na pochodne cząstkowe można znaleźć w pracy Deczky'ego [].

Jak wynika z powyższej dyskusji nasz problem redukuje się do minimalizacji nieliniowej funkcji n zmiennych. Istnieje wiele metod rozwiązania tego problemu tej klasy, jednak jedną z najbardziej efektywnych jest metoda Flecher-Powell'a. Opiera się ona na wykorzystaniu właściwości hesjanu funkcji n zmiennych. Hesjanem funkcji n zmiennych nazywany macierz utworzoną z pochodnych cząstkowych funkcji:

(3.64)

Zgodnie z twierdzeniem Schwarza o równości pochodnych cząstkowych, hesjan jest macierzą symetryczną. Ponadto badając określoność hesjanu wyliczonego dla pewnego wektora x_k, możemy określić czy funkcja f(x) osiąga dla wektora x_k ekstremum. Algorytm Flecher-Powell'a opiera się na korzystnych właściwościach hesjanu,jeżeli dokonujemy poszukiwania minimum funkcji kwadratowej n zmiennych. Jeżeli do punktu minimum zbliżamy się ciągiem kroków w kierunkach poszukiwań odpowiednio miedzy sobą uzależnionych, zwanych kierunkami sprzężonymi, to minimum funkcji osiągniemy po co najwyżej n krokach. Kierunki sprzężone są układem n niezależnych wektorów ξ względem dodatnio określonego hesjanu A, jeżeli spełniają warunek:

(3.65)

Przy znajdowaniu minimum funkcji f(x) algorytmem Flecher-Powell'a posługujemy się następującym algorytmem iteracyjnym przy przejściu od punktu x(k) od punktu x(k+1):

(3.66)

w którym α_k^* wybieramy tak, aby funkcja f(x) osiągała minimum na kierunku D_kg_k. Macierz D_k jest przybliżeniem macierzy odwrotnej hesjanu policzonym w punkcje x(k), zaś g_k jest gradientem funkcji f(x) również policzonym w punkcje x(k). Po wyznaczeniu nowego punktu x(k+1) dokonujemy modyfikacji macierzy D_k zgodnie ze wzorami podanymi w []. W przypadku projektowania filtru wektor g_k zawiera pochodne cząstkowe funkcji odległości względem współrzędnych polarnych położeni zer i biegunów filtru. Wartość początkową położenia biegunów możemy wybrać dowolnie lub użyć prototypu filtru uzyskanego np. za pomocą metody Yule-Walkera. Proces iteracyjny skończy się, jeżeli funkcja odległości osiągnie określoną małą wielkość.

Jak wspomnieliśmy wcześniej poszukujemy minimum funkcji nieliniowej, więc właściwość osiągnięcia minimum w co najwyżej n krokach nie jest już prawdziwa. Jednak, jeżeli funkcja jest kwadratowa w pobliżu minimum, metoda posiada informację o tym poprzez obliczenie hesjanu w tym punkcie i zbieżność zostanie osiągnięta dość szybko. Ważnym problemem, który pojawia przy projektowaniu jest stabilność uzyskanego filtru. Oznacza to, że interesują nas tylko takie wartości zawarte w wektorze A, dla których odległość bieguna do środka układu współrzędnych jest mniejsza od jedności. Jest to osiągnięte poprzez takie ograniczenie kroku algorytmu α_k^*, aby przy przesunięciu w każdym kierunku nowy punkt znajdował się w stabilnym obszarze P. W takim wypadku algorytm osiągnie zbieżność, jeżeli istnieje lokalne minimum wewnątrz obszaru stabilnego. Deczky w swojej pracy dowiódł, że wewnątrz obszaru stabilnego musi istnieć takie ekstremum lokalne. Algorytm zaproponowany przez Deczky'ego proponuje następujący sposób zagwarantowania stabilności projektowanego filtru. Otóż, jeśli nowy krok iteracji prowadzi do przekroczenia granicy stabilności, to jeżeli taki krok podzielony przez dwa, nowy punkt będzie leżał na pewna wewnątrz obszaru stabilnego. Gwarantuje to nam stabilność projektowanego filtru oraz zbieżność algorytmu.

3.5. Metody wyznaczania współczynników poprzez modelowanie zadanych charakterystyk

Omawiana dotychczas metody opierały się na projektowaniu filtru w dziedzinie częstotliwości. Poniżej omówione metody opierają się na projektowaniu filtru w dziedzinie czasu, na podstawie jego odpowiedzi impulsowej (metoda Padé oraz metoda Prony'ego) lub odpowiedzi filtru na dowolny sygnał rzeczywisty (metoda Steiglitz-McBride'a). Wszystkie te metody jako kryterium błędu używają miary średniokwadratowej.

3.5.1 Metoda Padé

Załóżmy, że dana jest odpowiedź impulsowa filtru h_d(n) dla n ≥ 0, zaś sam filtr jest opisany transmitancją:

(3.67)

gdzie h(n) jest odpowiedzią impulsową filtru. Filtr ma L = M + N +1 parametrów, mianowicie współczynników filtru {a_k} oraz {b_k}, które muszą być tak dobrane, aby minimalizować wybrane kryterium błędu. Jak wspomnieliśmy we wstępie w takim przypadku często używa się kryterium błędu średniokwadratowego. Załóżmy, że minimalizujemy sumę błędu podniesionego do kwadratu

(3.68)

pod względem współczynników filtru, gdzie U jest górną granicą sumowania.

Generalnie h(n) jest nieliniową funkcją współczynników filtru, więc minimalizacja błędu powoduje konieczność rozwiązania układu równań nieliniowych. Jednakże, jeżeli ustalimy górną granicę sumowania U = L - 1, jest możliwe dopasowanie pożądanej odpowiedzi impulsowej h_d(n) z odpowiedzią projektowanego filtru h(n) dla 0 ≤ n ≤ M + N. Jest otrzymywane w następujący sposób.

Równanie różnicowe projektowanego filtru jest dane wzorem:

(3.69)

Załóżmy, że pobudzeniem filtru jest impuls jednostkowy. Wtedy odpowiedź filtru y(n) = h(n) i dlatego równanie (3.69) zmienia się na:

(3.70)

Ponieważ δ(n - k) = 0 z wyjątkiem n = k (3.70) redukuje się do:

(3.71)

Natomiast dla n > M otrzymujemy:

(3.72)

Zbiór liniowych równań zawartych w (3.71) i (3.72) może być użyty w celu otrzymania współczynników filtru {a_k} i {b_k}. Jeżeli podstawimy h(n) = h_d(n) do równania (3.72) możemy wyznaczyć współczynniki {a_k} filtru. Korzystając z otrzymanego wyniku podstawiamy {a_k} do równania (3.71), aby otrzymać współczynniki mianownika filtru. W wyniku tych działań odpowiedź impulsowa h(n) zaprojektowanego filtru doskonale pokrywa się ze zadaną odpowiedzią impulsową dla pierwszych L wartości.

Stopień, w którym ta technika projektuje zadowalające filtry, zależy w części od wybranej liczby współczynników filtru. Ponieważ dokładność odtworzenia zadanej odpowiedzi impulsowej tylko rośnie wraz ze wzrostem liczby współczynników filtru, im bardziej złożony filtr, tym lepsza aproksymacja h_d(n). Jednakże jest to także główną wadą tego, ponieważ niewielka poprawa właściwości filtru powoduje spory wzrost liczby zer i biegunów, co po pierwsze jest nieekonomiczne, a po drugie im wyższy rząd tym większa wrażliwość na skutki kwantyzacji.

3.5.2. Metoda Prony'ego

Ta metoda, podobnie jak metoda Padé, również opiera się na równaniach (3.71) i (3.72). Nie zakłada ona jednak, że przez podstawienie h_d(n) do równania (3.72) otrzymany idealną odpowiedź impulsową, lecz przybliżoną:

(3.73)

Naszym zadaniem jest taki wybór współczynników {a_k}, żeby zminimalizować sumę błędów kwadratowych pomiędzy zadaną odpowiedzią impulsową h_d(n), a jej wartością przybliżoną h_d^'(n) dla n > M. Stąd mamy:

(3.74)

Minimalizacja błędu ε względem współczynników {a_k} prowadzi do ułożenia układu równań liniowych:

(3.75)

gdzie r_hh(k, l) jest zdefiniowane jako:

(3.76)

Rozwiązanie tego układu równań liniowych prowadzi do otrzymania współczynników mianownika transmitancji. Współczynniki licznika {b_k} otrzymujemy poprzez podstawienie do równania (3.72) w miejsce odpowiedzi impulsowej h(n) zadanej odpowiedzi impulsowej h_d(n) oraz otrzymanych wcześniej współczynników mianownika. Postępujemy więc tak samo jak w metodzie Padé. Metoda minimalizacji błędu średniokwadratowego prowadzi do dobrej estymacji parametrów mianownika {a_k}, lecz nie jest efektywna jeżeli chodzi o dobór współczynników {b_k}, głownie z tego powodu, że ich obliczenie nie bazuje na metodzie minimalizacji błędu średniokwadratowego.

3.5.3. Metoda Shanks'a

W tej metodzie użyjemy minimalizacji błędu średniokwadratowego dla również dla współczynników licznika transmitancji filtru {b_k}. Metoda ta używa metody Prony'ego do ustalania parametrów mianownika transmitancji filtru. Następnie zakłada się, że filtr można rozłożyć na dwa filtry: transmitancja pierwszego zawiera same bieguny obliczone przy pomocy metody Prony'ego, zaś drugi zawiera same zera. Sytuację tę przedstawiono na poniższym rysunku

Rys. 3.6. Rozkład filtru na dwie części w metodzie Shanks'a.

Transmitancja filtru zawierającego same bieguny opisana jest wzorem:

(3.77)

Odpowiedź tego filtru na pobudzenie impulsem jednostkowym wynosi:

(3.78)

Jeśli ciąg v(n) jest podany na wejście filtru, którego transmitancja posiada same zera, to jego odpowiedź jest dana następującym wzorem:

(3.79)

Teraz zdefiniujmy sygnał błędu jako:

(3.80)

I konsekwentnie współczynniki {b_k} mogą zostać wyznaczone w sensie kryterium minimalizacji błędu średniokwadratowego, to jest poprzez minimalizacje wyrażenia:

(3.81)

W wyniku otrzymujemy liniowy układ równań, z którego możemy obliczyć parametry {b_k} w postaci:

(3.82)

gdzie:

(3.83)

3.5.4. Metoda Steiglitz-McBride'a

Metoda opiera się na identyfikacji systemu na podstawie obserwacji odpowiedzi systemu na zadane powodzenie. Zakładamy, że badany układ charakteryzuje się transmitancja w postaci ilorazu dwóch wielomianów:

(3.84)

Kalman zaproponował następujący układ do znajdowania współczynników badanego układu:

Rys 3.7. Algorytm Kalman'a.

Jeżeli x i w są znanymi sekwencjami skończonej długości, to jeżeli następujący sposób minimalizacji błędu zostanie zastosowany:

(3.85)

gdzie sumowanie jest przeprowadzone na całej długości sekwencji wejściowej, zaś konturem całkowania jest okręg jednostkowy. Rozwiązanie tego problemu jest stosunkowo proste. Błąd w chwili j jest dany poprzez:

(3.86)

gdzie wektor współczynników δ i wektor wejścia wyjścia są zdefiniowane następująco:

(3.87)

Następnie sumujemy wszystkie błędy i liczmy gradient względem wektora współczynników. Otrzymany wynik przyrównujemy do zera, aby otrzymać minimum błędu:

(3.88)

Podstawienie (3.86) do (3.88) daje w rezultacie:

(3.89)

Ostatecznie wektor współczynników filtru jest otrzymywany jako:

(3.90)

gdzie Q jest macierzą korelacji o wymiarach 2N na 2N,zaś c wektorem korelacji o długości 2N:

(3.91)

Niestety nie istnieje żadna interpretacja fizyczna układu z rysunku 3.7. Steiglitz i McBride w swojej pracy [] zaproponowali inny układ do identyfikacji nieznanego systemu na podstawie znajomości odpowiedzi układu na zadane pobudzenie .

Rys 3.8. Algorytm zaproponowany przez Steiglitz i McBride.

Interpretacja fizyczna tego algorytmy jest wyraźna. Kiedy błąd e(n) osiągnie minimum filtr N(z)/D(z) będzie najlepszą aproksymacją nieznanego układu H(z). Niestety minimalizacja tego błędu jest tak wysoce nieliniowym problemem, że autorzy nie znaleźli rozwiązania tego problemu wprost. Zamiast tego zaproponowali algorytm iteracyjny, w wyniku otrzymujemy minimalizację błędu dla podanego przez nich sposobu identyfikacji systemu.

W trakcie iteracji powtarzany jest algorytm Kalman'a, przez co jego implementacja jest łatwa. Złożoność obliczeniowa nowej metody jest znacznie większa, lecz dostarcza lepszej aproksymacji. Idea procedury iteracyjnej przedstawiona jest na rysunku 3.9 . Pierwsza iteracja jest przeprowadzana przy pomocy algorytmu Kalman'a przy użyciu oryginalnych ciągów x i w. Rezultat tej operacji jest pierwsze przybliżenie N(z) i D(z). Nazwijmy je N₁(z) i D₁(z). Następnie oryginalne ciągi wejściowy i wyjściowy, x i w, są teraz filtrowane przez filtr cyfrowy o transmitancji 1/N₁(z), dając nowe przefiltrowane ciągi wejściowy i wyjściowy x' i w'. Te sekwencje są następnie użyte, zamiast oryginalnych, do wyznaczenia kolejnych filtrów N₂(z) i D₂(z) przy użyciu algorytmu Kalman'a. Filtr cyfrowy 1/N₂(z) jest użyty do otrzymania nowych ciągów x' i w' i tak dalej. Nie istnieje żaden ostry warunek zbieżności algorytmu. Algorytm uznaje się za zbieżny, jeśli D_i(z)/D_i-1(z) jest bliskie jedności. Dokładność aproksymacji może być polepszona, jeśli pochodne cząstkowe funkcji błędu po współczynnikach filtru będą równe zeru w momencie osiągnięcia zbieżności przez algorytm. Dla przedstawionego powyżej algorytm nie jest to prawdziwe, ponieważ omawiane pochodne są różnymi funkcjami parametrów modelu. W [] jest podana metoda rozwiązanie tego problemu. Podstawowa różnica pomiędzy tymi metodami polega na tym, że w ulepszonej metodzie we wzorze (3.88) wektor wyjścia wejścia q_j zastąpiono wektorem p. danym w następujący sposób:

(3.90)

gdzie x_j^| jest przefiltrowaną przez filtr 1/D(z) sekwencją wejściową x, natomiast v jest przefiltrowaną przez cały filtr N(z)/D(z) zmodyfikowaną sekwencją x_j^|. Tak więc jako we wzorze stosujemy właściwą odpowiedź filtru, a nie jak w metodzie wcześniej sekwencje wyjściową badanego systemu. Różnice w dokładności pomiędzy omówionymi metodami są niewielkie, ale metoda druga daje pewność osiągnięcia globalnego minimum błędu. Wadą metody drugiej jest konieczność zadania warunków początkowych. Algorytm może nie być zbieżny, jeżeli warunki początkowe są dalekie od optimum. Dlatego jako warunki początkowe podaje się współczynniki filtru obliczone na przykład metodą Prony'ego lub przed rozpoczęciem iteracyjnego powtarzania metody 2 wykonać kilka iteracji przy pomocy metody 1.

Rys 3.8. Idea metody iteracyjnej Steiglitz-McBride.

3.6. Porównanie metod projektowania filtrów IIR

W tablicy 3.1. zawarto porównanie metod projektowania filtrów IIR.

Tablica 3.1. Porównanie metod filtrów IIR.

	Prototypu analogowego	Yule-Walkera	p-tej normy	modelowania
Kontrola granic pasm	Dobra, dla filtru eliptycznego bardzo dobra	Słaba	Bardzo dobra	Zależy od danych wejściowych na ogół słaba
Pasmo przejściowe	Zależy od prototypu dla eliptycznego najwęższe	Szerokie	wąskie	Charakterystyka często jest często zniekształcona
Tłumienie	Optymalne dla filtrów Czebyszewa typ 2 i eliptycznych	Zbyt małe	Bardzo dobre	Raczej słabe
Zafalowania	Optymalne dla filtrów Czebyszewa typ 1 i eliptycznych	Zdarzają się duże przerzuty przy granicy pasma	Małe	Mogą być duże
Funkcja wagi	Nie	Nie	Tak	Nie
Estymacja rzędu	Tak	Nie	Nie	Nie
Stabilność	gwarantowana	Gwarantowana	Gwarantowana	Zależy od danych wejściowych
Filtry o wielu pasmach	Nie	Tak	Tak	Tak
Kształt charakterystyki fazowej	Zależny od prototypu	Zniekształcona	Możliwość uzyskania filtru o zadanej charakterystyce fazowej	Zniekształcona
Szacunkowy unormowany czas projektowania	8	1	122	3-14

Największe możliwości w dziedzinie projektowania filtrów ma metoda minimalizacji p. normy. Jej główne zaletą: możliwość uzyskania narzuconego kształtu charakterystyce amplitudowej lub fazowej. Główną wadą jest brak kryterium umożliwiającego oszacowanie rzędu filtru, który spełniłby narzucone przez projektanta wymagania. W tym polu przoduje metoda prototypu analogowego, ponieważ dla każdego rodzaju prototypu istnieje wzór pozwalający oszacować rząd filtru. Oprócz tego filtry eliptyczne cechują się bardzo dobrymi właściwościami, jeżeli chodzi o parametry charakterystyki amplitudowej. Metoda Yule-Walkera generuje filtry o niskiej jakości, lecz jest szybka i wymaga podania prostego zbioru parametrów wejściowych. Metody modelowania mają ograniczone możliwości, jeśli chodzi o projektowanie filtrów cyfrowych. Jest to skutkiem wymaganych przez nie parametrów filtru w dziedzinie czasu, których dobranie, aby otrzymać stabilny filtr o dobrych właściwościach nie jest łatwe.

4. Wybrane problemy implementacji filtrów cyfrowych

Nie można rozpatrywać problemu projektowania filtrów cyfrowych bez rozważenia wpływu implementacji na parametry filtru. Podczas projektowania zakładamy, że zarówno współczynniki filtru, jak i sygnał wejściowy ma nieskończoną precyzję dla każdej wartości pomiędzy -∞, a +∞. W praktyce mogą one przybierać tylko dyskretne własności ponieważ rejestry w których są przechowywane mają skończoną długość. Proces kwantyzacji powoduje nieliniowość równań różnicowych opisujących układ czasu dyskretnego. Te nieliniowe równania są niemal niemożliwe do opisania i analizy. Jednakże, jeśli wielkość błędu kwantyzacji jest mała w stosunku do wartości sygnału i parametrów filtru, prosta teoria aproksymacji bazująca na modelu statystycznym może zostać wprowadzona. Używając modelu statystycznego możliwe jest łatwe uwzględnienie efektu kwantyzacji i otrzymanie wyników, które mogą zostać sprawdzone eksperymentalnie. Głównymi źródłami błędu w filtrach cyfrowych są :

1. kwantyzacja współczynników filtru,

2. szumy konwersji analogowo-cyfrowej,

3. kwantyzacja operacji arytmetycznych,

4. zjawisko pojawiania się cykli granicznych.

Powyżej wymienione błędy zostaną omówione w kolejnych podrozdziałach. Ponadto omówiono metodę pomiaru jakości zaimplementowanego filtru.

4.1. Kwantyzacja współczynników filtru

Kwantyzacja współczynników filtru powoduje zmianę kształtu charakterystyki częstotliwościowej filtru, ponieważ powoduje ona zmianę położenia zer i biegunów filtru. Rozważmy przypadek filtru posiadającego same bieguny. Mianownik jego transmitancji można wyrazić w następującej formie :

(4.1)

gdzie pk są biegunami filtru. Po kwantyzacji mianownik transmitancji ulegnie zmianie do postaci :

(4.2)

gdzie p są biegunami filtru po kwantyzacji. Poprzez Δp_k oznaczmy zmianę położenia bieguna powstała w wyniku dyskretyzacji współczynników. Teraz odnieśmy zniekształcenie Δp_k do błędu kwantyzacji każdego współczynnika transmitancji a_k. Błąd Δp_i może być wyrażony jako :

(4.3)

gdzie pochodne cząstkowe reprezentują zmianę położenia bieguna p_k na wskutek zmiany wartości współczynnika a_k, zaś Δa_k oznacza zmianę konkretnego składnika mianownika transmitancji, powstałą na wskutek kwantyzacji. Pochodne cząstkowe we wzorze (4.3) reprezentują wrażliwość bieguna p_k na zmianę współczynnika a_k. Całkowite przesunięcie bieguna p_i będzie wyznaczone poprzez sumowanie wpływu zmiany wartości wszystkich współczynników. Pochodne cząstkowe można wyliczyć z wyrażenia na pochodną D(z) po współczynniku a_k w punkcie z = p_k :

(4.4)

Interesujące nas pochodne cząstkowe można znaleźć poprzez podzielenie przez siebie odpowiednich pochodnych D(z).Jako wynik otrzymujemy gotowy wzór określający nowe położenie filtru po kwantyzacji :

(4.5)

gdzie symbolem S oznaczono odpowiednie wrażliwości. Identyczne wyrażenie można otrzymać badając wrażliwość zer na zmianę współczynników licznika transmitancji. Wyrażenie (p_i - p_l) w mianowniku wzoru (4.5) reprezentuje wektor w przestrzeni z pomiędzy biegunem p_i, a biegunem p_l . Jeżeli, bieguny leżą blisko siebie, np. jak w przypadku filtrów wąskopasmowych, to odległość pomiędzy nimi, a biegunem p_i jest mała. Te małe odległości będą musiały powodować duży błąd Δp_i. Błąd Δp_i jest minimalizowany, jeżeli odległość | p_i - p_l | jest maksymalizowana. Ten warunek może być spełniony poprzez realizację filtru w postaci sekcji pierwszego lub drugiego rzędu. Jednakże, w ogólności realizacja filtru w postaci sekcji pierwszego rzędu pociąga za sobą konieczność przedstawienia współczynników filtru w postaci zespolonej oraz użycia arytmetyki zespolonej w celu realizacji filtru. Tego problemu można uniknąć grupując sprzężone bieguny zespolone w postaci sekcji drugiego rzędu. A ponieważ sprężone bieguny zespolone zazwyczaj leżą daleko od siebie, błąd Δp_i jest minimalizowany. Należy dodać, że powyżej wyrażone wzory dotyczą tylko realizacji filtru w formie bezpośredniej. Przedstawiona powyżej forma zapisu błędu Δp_i jest niewygodna z kilku powodów. Nie daje ona w sposób czytelny wrażliwości dane bieguna, ponieważ jest ona wyrażona w postaci liczby zespolonej. Ponadto jej obliczenie jest dość skomplikowane. Mitra w swojej pracy [] przestawił wyrażenie wrażliwości biegunów w postaci dwóch składowych : wrażliwości modułu bieguna na zmianę wartości współczynnik i wrażliwości kąta fazowego bieguna. Wszystkie te liczby są rzeczywiste i łatwe w interpretacji. Ponadto podał sposób na ocenę wrażliwości dowolnej struktury filtru. Załóżmy, że mamy daną strukturę filtru cyfrowego zrealizowaną przy pomocy R mnożników o współczynnikach α_k, które są wieloliniową funkcją współczynników transmitancji D(z). Zmiana współczynników α na α_k + Δα_k w związku z kwantyzacją, powoduje zmianę Δa_i współczynników wielomianu D(z) może być opisana jako :

(4.6)

Równanie (4.6) możemy zapisać w postaci mnożenia macierzy :

(4.7)

gdzie C jest macierzą zawierającą wrażliwości współczynników wielomianu D(z) na zmianę wartości współczynników badanej struktury :

(4.8)

Przesunięcie biegunów rozważanej struktury otrzymujemy podstawiając wzór (4.7) do wzoru (4.5) :

(4.9)

gdzie w macierzy S zawarte są wrażliwości położenia biegunów wielomianu D(z) na zmianę współczynników wielomianu.

4.2. Błędy konwersji analogowo - cyfrowej

Przetworniki A/C używane w aplikacjach cyfrowego przetwarzania sygnałów najczęściej używają reprezentacji binarnej wyniku w postaci kodu uzupełnień do dwóch. Dla przetwarzania bipolarnych sygnałów analogowych przetworniki generują bipolarne wyjście w postaci stałoprzecinkowego ułamka. Próbki cyfrowe generowane poprzez konwerter A/C są binarną reprezentacją skwantyzowanej odpowiedzi idealnego układu próbkującego o nieskończonej precyzji. Błąd kwantyzacji pojawia się z powodu skończonej długości słowa wyjściowego. Jeśli słowo wyjściowe ma długość (b + 1) bitów włączając w to bit znaku, całkowita liczba dostępnych poziomów wynosi 2^b+1. Zakres, w którym przetwornik A/C poprawnie przetwarza sygnał analogowy, nazywamy przedziałem pełen skali przetwornika i jest dany następującą zależnością :

(4.10)

gdzie δ jest krokiem kwantyzacji, to znaczy wielkością przedziału pomiędzy dwoma kolejnymi poziomami kwantyzacji. Dla liczby bipolarnej zapisanej na (b+1) bitach w kodzie uzupełnienia do dwóch, przedział prezentowanych liczb wynosi :

(4.11)

Charakterystyka przetwornika analogowo cyfrowego przedstawiona jest na rysunku:

Rys 4.1 Charakterystyka przejściowa 3 bitowego przetwornika A/C []

Jak wynika z charakterystyki przejściowej przetwornika błąd konwersji mieści się w granicach :

(4.12)

Zakładając, że próbka mieszcząca się dokładnie pomiędzy dwoma poziomami kwantyzacji należy do najbliższego wyższego poziomu oraz zakładając, ze wejściowy sygnał analogowy mieści się wewnątrz przedziału pełnej skali przetwornika. Każdy przebieg wyjściowy przekraczający ten przedział, powoduje, że amplituda błędu rośnie liniowo wraz ze wzrostem amplitudy sygnału wejściowego. Powodem tego zjawiska jest nasycenie wyjścia przetwornika do maksymalnej wartości dla 1-2^-b, jeśli sygnał wejściowy jest dodatni, albo do wartości -1 dla ujemnego sygnału wejściowego. W tej sytuacji błąd e[n] nazywany jest szumem nasycenia lub szumem nadmiaru. Obcinanie wyjścia przetwornika A/C powoduje zniekształcenie sygnału, co powoduje wysoko niepożądane efekty, przez co musi być unikane poprzez dokonywanie skalowania sygnału wejściowego, w taki sposób, aby być pewnym, że sygnał wejściowy zmieści się w przedziale pełnej skali przetwornika. Jeżeli założymy, że sygnał wejściowy jest procesem scholastycznym, którego gęstość prawdopodobieństwa jest opisana rozkładem Gaussa o wariacji σ , to założenie, że R_FS = 16σ jest, dla większości aplikacji, więcej niż wystarczające, aby zapewnić brak obcięć przy konwersji. W dalszej części pracy założono, że podczas konwersji nie występuje nasycenie się przetwornika.

Dokładna analiza błędu konwersji jest trudna, ponieważ charakterystyka przejściowa przetwornika jest nieliniowa i zazwyczaj nie znamy kształtu przebiegu wejściowego. Dlatego dużym ułatwieniem w analizie jest przyjęcie, że błąd kwantyzacji e[n] jest sygnałem przypadkowym (szumem), który jest dodawany do sygnału wejściowego. Dla uproszczenia analizy załóżmy, że :

1. ciąg błędu e[n] jest procesem stacjonarnym w szerokim sensie, prawdopodobieństwo pojawienia się konkretnej wartości błędu jest takie samo dla wszystkich możliwych wartości,

2. ciąg błędu jest nieskorelowany z sygnałem wejściowym,

3. ciąg wejściowy x[n] jest stacjonarnym procesem scholastycznym.

Założenie o nieskorelowaniu szumu z sygnałem wejściowym jest prawdziwe dla przetworników używających zaokrąglania w dowolnym formacie zapisu liczb binarnych lub obcinania liczb w formacie uzupełnienia do dwóch. Model statystyczny powoduje, że analiza szumu konwersji A/C staje się bardziej przejrzysta i otrzymane wyniki okazują się użyteczne w większości zastosowań.

Wpływ szumu addytywnego kwantyzacji na sygnał wejściowy jest dany poprzez stosunek sygnału do szumu :

(4.13)

gdzie σ_x² jest wariancją sygnału wejściowego odpowiadającą mocy sygnału, zaś σ_e² jest wariancją szumu kwantyzacji odpowiadającą jego mocy. Dla bipolarnego przetwornika o długości słowa wyjściowego (b + 1) wyrażenie na SNR możemy otrzymać w następującej wygodnej postaci :

(4.14)

To wyrażenie jest używane, aby określić minimalną długość słowa przetwornika, aby osiągnąć określony SNR. Wartość stosunku sygnału do szumu wzrasta o 6 decybeli dla każdego bitu dodanego do długości słowa. Dla stałej długości słowa rzeczywisty SNR zależy od wariacji sygnału wejściowego i przedziału pełnej skali przetwornika R_FS. Załóżmy, że SNR cyfrowego ekwiwalentu analogowego sygnału x[n] opisanego rozkładem Gaussa o wartości średniej równej zeru i wariancji σ_x, jest uzyskany przy pomocy przetwornika A/C o przedziale pełnej skali równym R_FS =Kσ_x. Wtedy :

(4.15)

Stosunek sygnału do szumu możemy poprawić poprzez skalowanie sygnału wejściowego. Rozważmy następujące skalowanie v[n] = A x[n]. Wówczas wariancja przeskalowanego wejścia wynosi σ_v² = A² σ_x². Stąd nowy SNR wynosi :

(4.16)

Dla danej długości słowa przetwornika SNR może być zwiększany poprzez skalowanie w górę sygnału wejściowego poprzez przyjęcie A > 1. Ale zwiększenie sygnału wejściowego zwiększa prawdopodobieństwo, że wartości próbek przekroczą R_FS i w wyniku tego wyrażenie (4.14) nie będzie dłużej prawdziwe, a sygnał cyfrowy będzie silnie zniekształcony. Natomiast skalowanie w dół, dzięki przyjęciu A < 1, spowoduje zmniejszenie wartości SNR. Dlatego jest koniecznością upewnienie się, że przedział wielkości sygnału analogowego jest jak to możliwe najbardziej bliski przedziałowi pełnej skali przetwornika A/C, aby otrzymać możliwie duży stosunek sygnału do szumu.

Aby określić propagację wejściowego szumu kwantyzacji na wyjście filtru cyfrowego, że filtr jest zrealizowany z nieskończoną precyzją. W praktyce, kwantyzacja operacji arytmetycznych powoduje powstawanie błędów wewnątrz struktury filtru, które także propagują na wyjście filtru cyfrowego i objawiają się jako szum. Zakładając, że wewnętrzne źródła szumu są niezależne od wejściowego szumu kwantyzacji, ich efekt będzie analizowany oddzielnie i dodany do wejściowego. Wykorzystując liniowość filtru oraz założenie, że błąd kwantyzacji i sygnał wejściowy jest nieskorelowany, może rozpatrzyć propagację szumu kwantyzacji przez filtr niezależnie do własności sygnału wejściowego. Przyjmijmy, że filtr cyfrowy w odpowiedzi na pobudzenie szumem kwantyzacji wygeneruje przebieg v[n] :

(4.17)

Najważniejszym parametrem sygnału v[n] jest znormalizowana wariancja, określająca stosunek mocy szumu na wyjściu filtru do mocy szumu na jego wejściu :

(4.18)

Alternatywnie wzór (4.18) można zapisać w postaci :

(4.19)

lub

(4.20)

Dokładną wartość znormalizowanej wariancji możemy uzyskać, że wzoru (4.19). Aby uprościć całki zawarte w wyrażeniu, należy dokonać rozkładu transmitancji filtru h(z) na ułamki proste. W wyniku tej operacji wyrażenia pod całkami stają się dość proste i ich funkcje pierwotne są podane w [17]. Wnioskiem ze wzoru (4.19) jest stwierdzenie, że bieguny położone blisko okręgu jednostkowego powodują gwałtowne zwiększenie szumu wyjściowego. Dla rozwiązań wysokiej jakości, długość słowa rejestrów przechowujących próbki sygnału powinna być dłuższa, aby utrzymać szum poniżej określonego poziomu. Prostszą metodą obliczenia wariancji szumu wyjściowego jest policzenie go ze wzoru (4.20). Dla filtru stabilnego i przyczynowego odpowiedź impulsowa gwałtownie maleje do wartości bliskich zeru. Wtedy przybliżona wartość wariancji można uzyskać ze wzoru :

(4.21)

Aby obliczyć przybliżoną wartość wariancji szumu musimy dokonywać sumowania, aż kolejne wartości wariancji będą różniły się od siebie o małą liczbę.

Zmniejszenie szumu kwantyzacji przetwornika A/C możemy osiągnąć poprzez zastosowanie przetwornika sigma-delta. Konwertery te kształtują w ten sposób widmo szumu kwantyzacji, że w użytecznym paśmie jest ono znacznie niższe niż w przypadku klasycznego przetwornika. Efekt ten uzyskuje się dzięki pętli sprzężenia zwrotnego oraz stosowaniu nadpróbkowania

4.3. Kwantyzacja operacji arytmetycznych

Wpływ kwantyzacji operacji arytmetycznych silnie wpływa na jakość realizacji filtrów stałoprzecinkowych. Dwoma najważniejszymi efektami kwantyzacji w tych filtrach są :

1. konieczność odrzucania części wyniku mnożenia,

2. możliwość wystąpienia nadmiaru podczas wielokrotnego dodawania.

Jak wspomnieliśmy wcześniej analityczne rozwiązanie problemu szumu kwantyzacji jest trudne, między innymi z powodu nieliniowości problemu. Dlatego i w tym przypadku używamy modelu statystycznego. Błąd wprowadzany przez obcięcie wyniku mnożenia można uznać za szum dodawany do iloczynu. Model ten jest przedstawiony na rysunku przedstawionym poniżej :

Rys 4.2. Model statystyczny mnożnika.

Wyjście idealnego mnożnika v[n] jest kwantyzowane, aby otrzymać y[n], gdzie:

(4.22)

Dla celów analizy zakładamy, że :

1. ciąg błędu e[n] jest procesem stacjonarnym, którego gęstość prawdopodobieństwa jest jest funkcją stałą,

2. błąd e[n] nie jest skorelowany z sekwencją v[n], ciągiem wejściowym x[n], ani z wszystkimi innymi źródłami szumów.

Przypomnijmy, że założenie, że e[n] jest nieskorelowany z ciągiem v[n] jest prawdziwe tylko w przypadku, kiedy do uzyskania skróconego wyniku mnożenia wykorzystujemy zaokrąglenie lub obcięcie liczby zakodowanej w formacie uzupełnienia do dwóch. Dzięki tym założeniom możemy rozpatrywać każde źródło szumu oddzielnie, a następnie dokonać sumowania otrzymanych wyników. Aby być, dokładnym załóżmy, że liczba mnożników bezpośrednio dołączonych do i-tego węzła sumowania wynosi k_i. Wariancja szumu z każdej operacji sumowania niech wynosi σ₀², dzięki założeniu o stałej gęstości prawdopodobieństwa, widmo mocy dla każdego mnożnika także wynosi N(ω) = σ₀². Superpozycja k_j źródeł szumu wchodzących do i-tego węzła daje całkowite widmo szumu e_i(n) w tym punkcje równe N_j(ω) = k_j σ₀². Wkład każdego źródła szumu w wyjściową gęstość widmową szumu jest dane jako k_iσ₀²|G_i(ω)|², gdzie G_i(ω) jest transmitancja filtru od i-tego węzła do wyjscia filtru. Gęstość widmową całkowitego szumu na wyjściu filtru e(n) jest sumą wszystkich gęstości widmowych źródeł szumów :

(4.23)

Wariancja szumu wyjściowego jest dana wzorem :

(4.24)

Stosując twierdzenie Parsevala oraz zakładając, że odpowiedzią impulsową transmitancji G_i(z) jest g_i(n) otrzymujemy :

(4.25)

Równania (4.23) i (4.25) można rozwiązać w sposób opisany w poprzednim rozdziale.

Rozważmy najprostszy przypadek filtru, którego transmitancja posiada jeden biegun :

Rys 4.3. Model dla filtru o jednym biegunie..

Transmitancja G(z) w tym przypadku jest równa transmitancji całego filtru H(z). Jego odpowiedź impulsowa jest dana następującym wyrażeniem :

(4.26)

Wariancja szumu obliczona ze wzoru (4.25) wynosi:

(4.27)

Moc szumu na wyjściu zależy od czynnika 1/(1-a²),. który silnie zwiększa się, gdy biegun zbliża się do koła jednostkowego. Jeżeli rozważymy połączenie kaskadowe filtrów pierwszego rzędu, to rozmieszczenie kaskad w szeregu ma znaczenie przy redukcji szumu. Podobnie jak w technice analogowej, na początku kaskady umieszczamy sekcję o biegunie położonym najdalej od koła jednostkowego, następnie dołączamy kolejno sekcje o biegunach znajdujących się coraz bliżej okręgu jednostkowego. W ten sposób możemy osiągnąć znaczącą redukcje szumu. Oba powyższe wnioski dotyczą także sekcji wyższego rzędu.

Jeżeli implementujemy filtr używając arytmetyki stałoprzecinkowej musimy brać pod uwagę możliwość przekroczenia, w wyniku wielokrotnego dodawania, zakresu liczb, które możemy przedstawić przy pomocy używanego słowa binarnego. Ten rodzaj zniekształcenia nazywamy błędem nadmiaru (ang. overflow). W przypadku filtru FIR wystąpienie błędu nadmiaru jest niekatastrofalnym, ponieważ wpływa tylko na wynik bieżącej próbki. Natomiast dla filtrów IIR wystąpienie nadmiaru jest katastrofalne, ponieważ wpływa nie tylko na bieżącą próbkę wyjściową, ale także na wszystkie następne. Powodem takiego zachowania się filtrów IIR jest istnienie pętli sprzężenia zwrotnego w strukturze tego filtru. Sposobem na uniknięcie błędu nadmiaru jest takie przeskalowanie współczynników filtru, aby sygnał w każdym węźle filtru nie przekroczył zakresu reprezentowanych liczb. Większość filtrów stałoprzecinkowych posiada współczynniki w postaci ułamków, więc celem skalowania jest osiągnięcie warunku, że w każdym węźle układu wartość bezwzględna sygnału nie przekroczy jedynki. Należy pamiętać, że skalowanie może w znaczący sposób zmniejszyć amplitudę sygnału w niektórych węzłach układu, przez co pogorszyć stosunek sygnału do szumu. Skalowanie sygnału we wszystkich węzłach układu nie konieczne z dwóch powodów :

1. niektóre węzły są oddzielone od siebie tylko poprzez opóźnienia, więc jeśli sygnał w pierwszym węźle nie zawiera błędu nadmiaru, to nie może zawierać go każdy z następnych węzłów,

2. podczas sumowania więcej niż dwóch liczb, jeżeli całkowita amplituda poprawnej sumy jest wystarczająco mała, aby nie wystąpił błąd nadmiaru, wtedy jeśli używamy arytmetyki liczb zapisanych w kodzie uzupełnienia do dwóch, poprawna całkowita suma może być otrzymana bez względu na kolejność wykonywania obliczeń. Właściwość ta jest prawdziwa niezależnie od tego, jak wiele błędów nadmiaru pojawi się w sumach częściowych lub nawet w przypadku, jeżeli którekolwiek z wejść do węzła sumowania będzie obciążone błędem nadmiaru, jako rezultat pomnożenia przez współczynnik większy od jedności

Musimy więc, wybrać pewne węzły w sieci filtru cyfrowego, na które nałożymy ograniczenia na zakres zmian amplitudy sygnału. Tymi węzłami są węzły skoku (ang. branch nodes) posiadające jedno lub więcej wyjść prowadzących do mnożników. Oznaczmy sygnały w tych węzłach jako v_i(n) i transmitancje od wejścia do i-tego węzła poprzez F_i(z), natomiast jako F_i^*(z) transmitancje po przeskalowaniu. Ponadto załóżmy, że dostępny zakres zmian wartości bezwzględnej sygnału wejściowego x(n) i wszystkich węzłów skoku jest ograniczony przez pewną stałą M. Przyjmując zerowe warunki początkowe i pomijając niewielki błąd kwantyzacji otrzymujemy wartość sygnału w i-tym węźle skoku v_i jako splot :

(4.28)

Ponieważ amplituda sygnału węźle v_i ma być ograniczona uzyskujemy :

(4.29)

Jeśli przeskalowane wartości bezwzględnej sygnału w węzłach skoku v_i^*(n) także muszą spełniać ograniczenie na wielkość sygnału, musimy tak przeskalować sieć, aby zapewnić :

(4.30)

Nierówność (4.30) jest warunkiem koniecznym i wystarczającym na uniknięcie występowania błędów nadmiaru w strukturze filtru. Ograniczenie zawarte w (4.30) jest zbyt pesymistyczne dla większości praktycznych sygnałów. Ponadto jego użycie do skalowania powoduje duże obniżenie poziomu sygnału we węzłach filtru, co z kolei powoduje spadek stosunku sygnału do szumu. Kolejną wadą tego sposobu skalowania jest duża trudność wyznaczenia w sposób analityczny odpowiednich współczynników skalujących. Dlatego w praktyce używa się innych warunków skalowania takich jak :

1. norma L1
   

(4.31)

2. norma Czebyszewa
   

(4.32)

3. norma L2

(4.33)

Zauważmy, że we wzorach (4.31) i (4.32) używa się odpowiedzi impulsowej filtru, która jest łatwa do obliczenia. Rozważmy skalowanie filtru IIR w postaci kaskady sekcji drugiego rzędu.


Rys 4.4. Filtr IIR w postaci kaskady sekcji filtrów 2 rzędu.


Piewszym krokiem jest określenie węzłów w każdej sekcji drugiego rzędu, w których mogą pojawić się błędy nadmiaru. Tymi węzłami są węzły d_i(n) oraz y_i(n). W następnym kroku należy wyznaczyć transmitancję F(z) w każdym z tych węzłów. Posiadając te transmitancje możemy wyznaczyć czynnik skalujący G_k dla każdej z sekcji dla każdego z węzłów sieci. Zależnie od typu aplikacji i rodzaju sygnałów wejściowych jedna z trzech metod podanych powyżej może zostać wybrana, jako sposób na wyliczenie G_k. Alternatywną metodą jest wyznaczenie wszystkich trzech wartości G_k dla każdego z węzłów. Wtedy, posiadając obliczone różne wartości G_k dla wszystkich węzłów w sekcji, wybieramy czynnik o największej wartości, nawet jeśli będzie to bardzo ostrożne postępowanie w sensie skalowania. Następnie używamy wyznaczonego współczynnika dla skalowania współczynników tej sekcji.

Proces powinien być powtórzony dla każdej następnej sekcji drugiego rzędu. Jest podstawową sprawą zredukowania ryzyka przeskalowania filtru, poprzez upewnienie się, że wszystkie następne czynniki skalujące biorą pod uwagę skalowanie poprzednich sekcji. Istnieją dwie metody rozwiązania tego problemu :

1. obliczanie współczynników skalujący przy pomocy transmitancji wyznaczonych od wejścia filtru od określonego węzła. W tym przypadku przeskalowane współczynniki poprzednich sekcji muszą być użyte do wyznaczania transmitancji lub odpowiedzi impulsowej,

2. obliczanie czynników skalujący dla każdej sekcji oddzielnie, zmiany czynnika skalującego G_k dla sekcji k poprzez podzielenie go przez współczynnik skalujący poprzedniej sekcji G_k-1. W tej metodzie dokonuje się tylko skalowania współczynników licznika transmitancji sekcji drugiego rzędu. Czynnik skalujący dla bieżącej sekcji jest używany do skalowania współczynników poprzedniej sekcji. Czynnik skalujący dla pierwszej sekcji jest wprowadzany jako czynnik skalujący wejście. Ostatnia sekcja pozostaje nie zmieniona.

Istnieją także inne techniki skalowania. Jedna z nich polega na skalowaniu sygnału wejściowego poprzez stały czynnik skalujący, który obliczany jedna z metod przedstawionych powyżej. Istnieją dwa typowe miejsca umieszczenia tego czynnika. Może zostać użyty do przeskalowania wszystkich współczynników filtru lub użyty do skalowania wszystkich próbek wejściowych. Ta metoda powoduje większe zmniejszenie SNR niż pozostałe. Inną techniką jest skalowanie współczynników filtru podczas jego pracy. Kontrolny bit błędu nadmiaru jest sprowadzany, czy nastąpił nadmiar w którymkolwiek węźle którejkolwiek sekcji podczas wykonywania obliczeń. Jeśli nastąpił, to wyniki obliczeń w węzłach są obniżane, tak aby uniknąć nadmiaru. Ta technika prowadzi do bardzo dobrych wyników, ponieważ nigdy nie dopuści do powstania błędu nadmiaru i ryzyko przeskalowania filtru jest zminimalizowane, to jej koszt obliczeniowy jest bardzo wysoki i zazwyczaj nie stosuje się jej w praktycznych implementacjach na układach stałoprzecinkowych.

4.4. Cykle graniczne


Kwantyzacja powoduje, że filtr będący układem liniowym staję się układem nieliniowym z dwóch powodów:


1. przedział liczb reprezentowanych w systemach cyfrowych jest ograniczona. Jeśli suma dwóch lub więcej liczb daje wynik poza tym ograniczeniem, powstaje błąd nadmiaru, który powoduje zniekształcenie sygnału wyjściowego, zależne od poszczególnej implementacji systemu oraz od wartości sygnału wejściowego. W szczególności może powodować niestabilność systemu, prowadząc do dużych oscylacji na wyjściu, pomimo braku pobudzenia. Tego typu zachowanie się filtru cyfrowego jest nazywane cyklami granicznymi dużej skali lub cyklami granicznymi nadmiaru.

2. Drugi nieliniowy efekt jest spowodowany poprzez wymagane obcięcie długości słowa wyniku mnożenia. Teraz nie jesteśmy zainteresowani modelem rzeczywistego mnożnika, wprowadzonego wcześniej, ale możliwością niestabilności, prowadzącą do pojawienia się oscylacji na wyjściu filtru. Ponieważ powstałe cykle graniczne maja zwykle małe amplitudy, a ich amplitudy są zazwyczaj wielokrotnością kroku kwantyzacji δ, zostały one nazwane cyklami małej skali lub cyklami ziarnistymi (ang. granular).

Aby pokazać charakter cykli granicznych rozważmy filtr posiadający jeden biegun opisany następującym równaniem różnicowym :

(4.34)


gdzie biegun filtru znajduje się w z = a. W wyniku kwantyzacji filtr zostaje opisany następującym nieliniowym równaniem :

(4.35)


gdzie Q oznacza kwantyzację wyniku mnożenia. Załóżmy , że filtr jest zrealizowany w arytmetyce stałoprzecinkowej bazującej na liczbach b bitowych zakodowanych w kodzie uzupełnienia do dwóch. Natomiast kwantyzacja wyniku mnożenia polega na zaokrągleniu otrzymanego wyniku w górę. Odpowiedź impulsowa filtru idealnego zbiega ekwipotencjalnie do zera. Natomiast odpowiedź filtru rzeczywistego w stanie ustalonym jest przebiegiem okresowym o okresie, który zależy od położenia bieguna transmitancji. Jeśli biegun jest dodatni, to na wyjściu filtru uzyskamy przebieg stały. Natomiast w przypadku bieguna ujemnego znak próbki wyjściowej zmienia się co próbkę. Interesującą jest uwaga, że zachowanie się nieliniowego układu w czasie trwania cyklu granicznego odpowiada systemowi liniowemu z biegunem w z =1, gdy biegun jest dodatni, lub z biegunem w z = - 1 dla bieguna ujemnego. Kiedy sygnał wejściowy przyjmuje wartości zerowe, wyjście filtru, po pewnej liczbie iteracji, przechodzi w cykl graniczny. Także przy braku sygnału na wejściu oraz niezerowych warunkach początkowych mogą pojawić się cykle graniczne. Amplituda wyjścia filtru w czasie trwania cykli granicznych mieści się w pewnych granicach zwanych zakresem martwym (ang. dead band). Dla rozważanego filtru pierwszego rzędu zakres martwe dane jest wzorem :

(4.36)


Analiza zachowania się cykli granicznych sekcji drugiego rzędu jest o wiele bardziej złożona i istnieje większa różnorodność możliwych kształtów oscylacji. Rozważmy najczęściej spotykany przypadek, gdy sekcję bikwadratową tworzą bieguny sprzężone z = r e ^±jθ. Wówczas martwy jest dany wzorem :

(4.37)


Widać zatem, że wielkość cykli granicznych zależy od położenia biegunów. Im są one położone bliżej okręgu jednostkowego, tym wielkość cykli granicznych jest większa. Ponadto zwiększenie długości słowa b powoduje zmniejszenie się wielkości zakresu martwego. Zastosowanie przy kwantyzacji wyniku mnożenia, zamiast zaokrąglania, prostego obcinania iloczynu do b bitów pozwala na wyeliminowanie wielu, lecz nie wszystkich, cykli granicznych. Jednak użycie obcięcia powoduje powstanie błędu kwantyzacji o wartości średniej różnej od zera, co jest niepożądane w implementacji filtrów cyfrowych. Charakter cykli granicznych zależy od rodzaju struktury filtru użytej do implementacji filtru cyfrowego. W wypadku struktury równoległej cykl graniczny na wyjściu filtru jest sumą cykli granicznych powstałych w poszczególnych sekcjach drugiego rzędu. Nie istnieje interakcja pomiędzy cyklami powstałymi w poszczególnych sekcjach. W przypadku realizacji w postaci kaskady analiza cykli granicznych jest znacznie trudniejsza. Przy braku pobudzenia pierwsza sekcja generuje cykl graniczny, który jest z jest z kolei filtrowany przez następne sekcje kaskady. Jeśli częstotliwość cyklu granicznego leży w pobliżu częstotliwości rezonansowych pozostałych sekcji, to amplituda cyklu granicznego jest znacznie zwiększona. Dlatego powinniśmy unikać tego typu sytuacji.

W przypadku cykli granicznych powstających na skutek nadmiaru ich unikanie polega zasadniczo na zapobieganiu powstawania tego typu błędu, czy na dokonaniu skalowania filtru. Innym ważnym elementem przeciw działania tego typu błędom, jest zerowanie rejestrów przechowujących wyniki pośrednie w strukturze filtru, po włączeniu (resecie) urządzenia. Niestety te czynności nie gwarantują nie powstawania cykli granicznych nadmiaru. Efektywnym środkiem na uniknięcie osylacji powstałych w wyniku nadmiaru jest zmiana charakterystyki sumatora, w taki sposób, aby przy nastąpieniu nadmiaru wyjście sumatora nasyciło się w pobliżu wartości ±1.

4.5. Metoda pomiaru jakości zaimplentowanego filtru.

Filtr po implementacji jest tylko w przybliżeniu układem liniowym w związku z efektami skończonego słowa cyfrowego używanego wewnątrz struktury filtru. Taki układ może zostać zamodelowany poprzez równoległe połączenie układy liniowego, którego wyjście jest y_L(n) i drugiego, który generuje sygnał szumu kwantyzacji e(n). Taka separacja może zostać dokonana, ponieważ y_L(n) i e(n) są ortogonalne wobec siebie. Podczas, gdy system liniowy jest opisany charakterystyką częstotliwościową H(e^jω), natomiast drugi jest charakteryzowany widmową gęstością mocy Φ_ee(e^jω). Dzięki tak zwanej metodzie szumowej (ang. noise-loading metod) otrzymujemy próbki estymaty tych dwóch wartości

W tej metodzie dokonujemy sekwencji pomiarów, używając jako pobudzenia filtru sygnału v_l(n) będącego złożeniem sygnałów, które są okresowe dla n ≥ 0. W najprostszym przypadku mogą on zostać wygenerowane w następujący sposób :

(4.38)

gdzie wartości widma V_l są dane wzorem:

(4.39)

przy czy posiadaja one stała amplitudę dla wszystkich k i l, natomiast ϕ_l(k) jest zmienną przypadkową o gęstości prawdopodobieństwa równo dystrybuowanej na przedziale [-π,π] i statystycznie niezależnej od k i l. Warunek , że ϕ_l(k) = - ϕ_l(N - k) musi zostać wprowadzony, aby uzyskać rzeczywisty sygnał v_l(n). W wyniku otrzymujemy, że sekwencję v_l(n) jest w przybliżeniu normalnie dystrybuowana.

Testowany system będzie pobudzony przez wyliczone v_l(n) dla l =1 : L, gdzie L jest liczbą powtórzeń, która musi być tak dobrana, aby uzyskać odpowiednią dokładność wyniku. Po pewnym czasie otrzymany wyjściowy sygnał periodyczny opisany przez :

(4.40)

Minimalizując wartość wariancji sygnału błędu e_l(n) ze względu na H(e^jω) uzyskujemy przybliżony wzór na kształt charakterystyki częstotliwościowej :

(4.43)

gdzie Y_l(k) są współczynnikami widma przebiegu y_l(n). Gęstość widmowa szumu kwantyzacji jest dana poprzez wyrażenie :

(4.44)

,której wartość można oszacować wzorem

(4.45)

Dokładność estymacij wzrasta wraz ze wzrostem liczby iteracji L. Dzięki tej metodzie możemy ocenić jakość filtru pracującego w warunkach zbliżonych do rzeczywistych. Otrzymana charakterystyka częstotliwościowa dobrze oddaje rzeczywistą, ponieważ przy pomocy tej metody badamy rzeczywisty filtr bez żadnych założeń upraszczających analizę.







5.0 Projektowanie i implementacja filtrów cyfrowych przy pomocy programu Matlab i CCS firmy Texas Instuments.

Celem pracy było stworzenie narzędzia do projektowania i implementacji filtrów cyfrowych. To zadanie zostało wykonane przy pomocy pakietu MATLAB. Ponadto powinno chciano sprawdzić otrzymane wyniki na gruncie praktyki. Dlatego, przy pomocy narzędzia CCS firmy Texas Instruments, napisano używając języka C, kilkanaście filtrów cyfrowych używających arytmetyki stało- lub zmiennoprzecinkowej. Komunikacja pomiędzy częściami projektu odbywa się przy pomocy plików nagłówkowych.

5.1 Opis aplikacji wykonanych w MATLAB-ie

Pakiet MATLAB jest produktem firmy The Mathworks Inc. działającą na rynku od roku 1983. Jest on interaktywnym środowiskiem do wykonywania obliczeń naukowych i inżynierskich. Umożliwia on testowanie algorytmów, modelowanie i symulację, analizę i wizualizację danych, sygnałów oraz wyników obliczeń. Siła MATLAB'a opiera się na zoptymalizowanych działaniach na wektorach i macierzach rzeczywistych lub zespolonych. Dzięki istnieniu wielu operatorów i funkcji wykonywanie operacji na macierzach jest bardzo efektywne i szybkie. Dlatego podczas pisania programów w MATLAB-ie zaleca się jak najczęstsze ich używanie. MATLAB jest także językiem programowania dostarczającym typowych konstrukcji jak if lub for. Składnia języka jest zapożyczona z języka C. Posiada on kilka zalet pozwalających ułatwiających pracę programiście. Nie trzeba deklarować zmiennych. Zmienne wybranego typu określa użytkownik poprzez podanie jej nazwy i przypisanie do pewnego wyrażenia. MATLAB automatycznie rozpoznaje typ zmiennej, nadając jej typ wyniku wyrażenia, do którego zmienna została przypisana. Wielkość wektorów może być dynamicznie zwiększana podczas wykonywania programu bez żadnej dodatkowej ingerencji programisty, choć nie jest to zalecane, gdy znamy rozmiar wektora lub macierzy, ponieważ powoduje spowolnienie pracy programu. Przy używaniu wielu operacji wektorowych i macierzowych jest dokonywana niejawna indeksacja, co zmniejsza nakład pracy przy pisaniu aplikacji. Ponadto pisząc program w środowisku MATLAB-a można korzystać z prawie pełnych możliwości programowania zorientowanego obiektowo.

Kolejną zaletą MATLAB-a jest możliwość zdefiniowania własnych poleceń i algorytmów obliczeniowych poprzez pisanie M-plików. M-plik jest odpowiednikiem funkcji w języku C. Komunikują się one z przestrzenią roboczą MATLAB-a poprzez listę argumentów wejściowych i wyjściowych. Używanie M-pliku jest bardzo elastyczne, ponieważ wywołując funkcję nie trzeba zdefiniować wszystkich jej parametrów, w drastycznym przypadku nie trzeba określać ich wcale. Do programisty należy obsługa takich przypadków. W jednym M-pliku można zdefiniować wiele funkcji, ale tylko pierwsza z nich może być wywoływana z poziomu przestrzeni roboczej lub przez M-pliki zapisane w innych plikach. Zadeklarowane wewnątrz M-pliku zmienne mają charakter lokalny. Pakiet MATLAB zawiera debugger i profiler, z których pierwszy umożliwiają lokalizację błędów w funkcji, zaś drugi pozwala na optymalizację napisanego kodu. Korzystanie z M-plików zwiększa niezawodność napisanego kodu, szybkość jego wykonywania oraz umożliwia tworzenie bibliotek użytecznych programów. Przez twórców pakietu nazywane są toolbox-ami. Przy projektowaniu filtrów cyfrowych szczególnie użyteczne są dwa toolbox-y:

• Signal Processing Toolbox - zawierających wiele funkcji projektujących filtry cyfrowe IIR i FIR, jak też wiele funkcji obliczających podstawowe parametry filtru jak : charakterystykę częstotliwościową czy odpowiedź impulsową

• Filter Design Toolbox - zawierającym kilka zaawansowanych algorytmów projektowania filtrów, ale przede wszystkim zawierającym model filtru po kwantyzancji, który umożliwia badanie wpływu wielu czynników (np. struktury filtru, formatu zapisu liczb binarnych, czy długości słowa) na parametry zaimplementowanego filtru

Sławne są możliwości MATLAB-a w dziedzinie wizualizacji danych i wyników obliczeń. Pakiet zawiera wiele funkcji kreślących wykresy dwu- i trójwymiarowe. Całość grafiki w pakiecie jest zrealizowana poprzez obiektowo zorientowaną technikę Handle Graphics. Każdy obiekt graficzny posiada własną bazę danych, której pola zawierają informację o niemal wszystkich parametrach tego obieku. Uzytkownik poprzez ich zmianę ma wielką swobodę w kształtowaniu efektu wizualnego. Elementem składowym grafiki jest też GUI - interfejs graficzny użytkownika, który daje możliwość pracy interakcyjnej z wykorzystaniem okienek edycyjnych, przycisków lub menu. Duży ułatwieniem przy tworzeniu interfesu GUI jest program GUIDE znacznie zmniejszający nakład pracy przy tworzeniu interfejsu użytkownika GUI.

Pracując z pakietem MATLAB należy podkreślić wkład pracy, jaki włożyli twórcy w system pomocy. Jest on stworzony w sposób jasny i przejrzysty. Angielski używany do opisu jest zrozumiały, nawet dla użytkownika znającego ten język w stopniu niezaawansowanym. Komentarz do każdego polecenia nie ogranicza się tylko do jego opisu, prawie zawsze jest poparty przykładami. Często podawana jest literatura, z której korzystano przy pisaniu programu.

Przy pomocy pakietu MATLAB wykonano szereg interfejsów GUI, przy pomocy których można zaprojektować filtry cyfrowe za pomocą następujących metod :

1. filtry o skończonej odpowiedzi impulsowej :

1. metodą okien,

2. metodą próbkowania charakterystyki,

3. metodą minimalizacji błędu średniokwadratowego,

4. metodą minimalizacji błędu maksymalnego (Remeza),

5. metodą wymuszonej minimalizacji błędu średniokwadratowego,

6. metodą Butterworth'a,

2. filtry o nieskończonej odpowiedzi impulsowej :

1. metodą prototypu analogowego,

2. metodą Yule-Walkera,

3. metodą p-tej normy,

4. metodą modelowania.

Ponadto stworzono dwa inne interfejsy służące do :

1. sprawdzania , czy zaprojektowane filtry mieszczą się w granicach tolerancji nałożonych na przebieg charakterystyki amplitudowej

2. implementacji filtru

Współczynniki filtrów pomiędzy tymi aplikacjami są przenoszone poprzez zapisywanie ich jako zmienne w przestrzeni roboczej MATLAB'a. Z każdego programu projektującego filtr można wywołać interfejs GUI dokonujący implementacji. Ponieważ tak duża ilość aplikacji jest dość kłopotliwa dla użytkownika stworzono okienko pozwalające uruchamiać wszystkie stworzone programy z jednego miejsca. Obsługujący projekt musi zapamiętać tylko jedną nazwę interfejsu : gmenu. Wygląd tego okna jest przedstawiony na rysunku (5.1) :

Rys 5.1. Wygląd głównego menu projektu.

Przy pomocy menu Typ filtru wybieramy rodzaj cyfrowego, który chcemy zaprojektować: o skończonej odpowiedzi impulsowej (FIR) lub o nieskończonej odpowiedzi impulsowej (IIR). Metodę projektowania filtru wybieramy przy pomocy dwóch menu objętych wspólnym tytułem Rodzaj metody. W zależności od wybranego typu filtru jedno z tych menu jest niedostępne. Po naciśnięciu klawisza Ok jest otwierany interfejs GUI obsługujący wybraną przez użytkownika metodę projektowania.

Klawisze po prawej stronie interfejsu służą do :

1. przycisk Porownanie służy do uruchomienia okna przy pomocy, którego możemy porównać mieszczenie się charakterystyk filtrów, zaprojektowanych przy pomocy różnych metod, w ramach projektowych

2. klawisz Implementacja powodujący pojawienie się okienka w którym jest przeprowadzana implementacja filtru

3. przycisk Koniec zamykający główne menu oraz wszystkie okna otworzone przy jego pomocy.

Teraz zostaną omówione interfejsy GUI obsługujące kolejne metody projektowania filtrów cyfrowych.. Jako pierwsze zostaną omówione interfejsy projektujące filtry FIR. Ich opis rozpocznę od okna projektującego filtry FIR metodą okien czasowych.

Rys 5.2. Interfejs GUI projektujący filtry FIR przy pomocy okien czasowych.

Opis elementów składowych okna :

• Okno edycyjne Dlugosc filtru, służące do podania długości projektowanego filtru. To okno jest nieaktywne, jeżeli używamy estymacji rzędu filtru. Po dokonaniu estymacji w nieaktywnym oknie pojawia się rząd filtru, który powinien spełnić nasze wymagania podane podczas estymacji. Przy projektowaniu filtrów pasmowozaporowych lub górnoprzepustowych rząd filtru musi być parzysty, ponieważ funkcja MATLAB-a użyta w aplikacji tworzy filtry jedynie o symetrycznej odpowiedzi impulsowej. Jeżeli użytkownik będzie próbował zaprojektować wyżej wymienione typy filtrów z nieparzystym rzędem, to zostanie on automatycznie zwiększony o jeden.

• Przy pomocy przycisku Estymuj uruchamiamy dodatkowe okienko, które zostanie omówione później, służące do uzyskania przybliżonego rzędu filtru przy projektowaniu przy pomocy okna Kaisera. Otrzymujemy także wartość współczynnika okna oraz częstotliwości odcięcia. Wszystkie te wartości wpisywane są w odpowiednie pola edycyjne. Po naciśnięciu przycisku Zaprojektuj możemy obejrzeć wyniki estymacji.

• Na element nazwany Czestotliwosci odciecia składają się dwa okienka. Wpisujemy w nie częstotliwości graniczne filtru. W przypadku filtru dolno- i górnoprzepustowego aktywne jest jedno okienko, zaś w przypadku filtru pasmowoprzepustowego i pasmowozaporowego - dwa okienka.

• Za pomocą menu Typ filtr wybieramy typ filtru, który chcemy zaprojektować. Mamy do wyboru cztery typy filtrów: dolno- i górnoprzepustowy , pasmowoprzepustowy i pasmowozaporowy.

• Za pomocą menu Rodzaj okna wybieramy typ okna użytego do projektowania filtru. Od wyboru okna zależy tłumienie filtru w paśmie zaporowym oraz szerokość pasma przejściowego. Przy oknach Kaisera i Czebyszewa należy podać współczynnik okna.

• W okienko edycyjne Wspolczynnik okna wpisujemy współczynnik , potrzebny przy obliczaniu okien Kaisera i Czebyszewa. W przypadku okna Czebyszewa współczynnik ten oznacza wprost tłumienie w paśmie zaporowym, natomiast w przypadku okna Kaisera takiej zależności nie można oddać wprost. Im większy współczynnik okna Kaisera tym większe tłumienie , ale szerokość pasma przejściowego wzrasta

• Naciśnięcie przycisku Zaprojektuj powoduje zaprojektowanie filtru zgodnie z podanymi w innych elementach okna wartościami

• Czestotliwosc probkowania W tym okienku podajemy częstotliwość próbkowania.

• Menu Typy analizy służy do wyboru typu wyświetlanej charakterystyki. Możemy obejrzeć charakterystykę amplitudowa w trzech skalach osi Y:

• logarytmicznej

• podniesionej do kwadratu

• liniowej

Możemy także obejrzeć charakterystykę fazowa filtru, rozkład zer transmitancji filtru, charakterystykę opóźnienia grupowego filtru oraz odpowiedzi filtru na : impuls i skok jednostkowy.

• Jeśli znacznik "Pamietaj" jest zaznaczony, to zostają zapamiętane bieżące współczynniki filtru. Wszystkie charakterystyki będą wykreślane podwójnie, dopóki ten element nie zostanie wyłączony.

• Przycisk "Koniec" kończy pracę programu.

Ponadto w górnym menu znajdują się polecenia pozwalające na eksport współczynników filtru do przestrzeni roboczej MATLAB-a. Można eksportować zarówno filtr bieżący jak i zapamiętany. Istnieje także możliwość przejścia do narzędzia dokonującego implementacji. Znajdują się tam również polecenia pomocy : opis interfejsu oraz przykłady użycia narzędzia, które zostały tak dobrane, aby pokazać możliwości aplikacji.

Jeżeli w wyżej opisanym narzędzi użytkownik naciśnie klawisz Estymuj (estymacja rzędu jest możliwa tylko dla okna Kaisera) spowoduje to pojawienie się następującego okna :

Rys 5.3. Okno estymacji rzędu okna Kaisera.

W oknie znajduje się pięć elementów :

1. W okienko edycyjne Czestotliwosci graniczne pasm wpisujemy, oddzielone spacjami, częstotliwości graniczne pasm naszego filtru w kolejności od najniższej do najwyższej. W przypadku filtrów dolno- i górnoprzepustowego należy wpisać dwie częstotliwości, natomiast w przypadku filtrów : pasmowoprzepustowego i pasmowozaporowego trzy częstotliwości.

2. W menu Typ filtru wybieramy typ filtru spośród czterech standardowych możliwości.

3. Pole edycyjne Zafalowania służy do podania dopuszczalnych odchyłek pomiędzy idealną charakterystyką, a otrzymaną. Możemy więc powiedzieć, że dla pasma przepustowego wpisujemy dopuszczalne zafalowania, zaś dla pasm zaporowych podajemy oczekiwane tłumienie. Obie te wartości należy wpisywać w skali liniowej (nie przy użyciu skali decybelowej). Dla filtrów dolno- i górnoprzepustowego należy wpisać dwie wartości, natomiast dla pozostałych dwóch typów filtru trzy wartości.

4. Zawartość okna edycyjnego zatytułowanego Czestotliwosc probkowania określa częstotliwość próbkowania używaną w systemie, w którym zostanie użyty filtr.

5. Po naciśnięciu przycisku Estymuj obliczone następujące parametry : rząd filtru, jego częstotliwości odcięcia oraz wartość parametru β okna Kaisera. Dla tych wartości zaprojektowany filtr powinien spełnić narzucone na niego wymagania.

Policzone przez ten interfejs GUI parametry są automatycznie umieszczane we właściwych polach aplikacji projektujące filtr FIR przy pomocy okien czasowych.

Kolejnym narzędziem jest program projektujący filtry FIR metodą próbkowania w dziedzinie częstotliwości. Jego interfejs przedstawiono na rysunku (5.4) :

Rys 5.4. Okno programu do projektowania filtrów FIR metodą próbkowania charakterystyki.

Przeznaczenie poszczególnych elementów okna jest następujące :

• Okienko edycyjne Rzad filtru długość projektowanego filtru. Powinna być to całkowita liczba dodatnia. Podobnie, jak w poprzednio omawianym interfejsie filtry, których charakterystyka przy częstotliwości Nyquista równa jest zeru, muszą posiadać parzysty rząd. Jeśli tak nie jest, wartość rzędu jest automatycznie zwiększana o jeden.

• W pole edycyjne Czestotliwosci wpisujemy częstotliwości, w których podamy próbki charakterystyki amplitudowej. Częstotliwości mogą one przyjmować wartości pomiędzy zerem, a połową częstotliwości próbkowania. Pierwsza wartość musi równać się zeru, zaś ostatnia musi być równa połowie częstotliwości próbkowania. Kolejnym wymaganiem jest warunek, ze kolejne częstotliwości muszą być podawane w porządku niemalejącym. Oznacza to, można podać dwa takie same punkty częstotliwości. Takie zestawienie oznacza punkt zagięcia charakterystyki.

• Okno Probki amplitudy służy do pobrania, związanych w omawianymi wyżej częstotliwościami, próbek charakterystyki amplitudowej. Ich ilość musi być równa liczbie częstotliwości podanymi w polu Czestotliwosci.

• Pole edycyjne Liczba probek definiuje liczbę równo rozłożonych próbek interpolujących zadaną charakterystykę filtru pomiędzy podanymi próbkami charakterystyki częstotliwościowej. Ponieważ na tych próbkach dokonuje się w trakcie wykonywania algorytmu odwrotnej szybkiej transformaty Fouriera, liczba podana w tym polu powinna być potęgą liczby dwa. Im większa liczba próbek, tym lepszy otrzymamy wynik.

• Jeżeli zdefiniowaliśmy powtarzające się punkty częstotliwości, to algorytm umieszcza wokół tego punktu obszar, aby pasmo przejściowe było gładkie i w miarę strome. Pole edycyjne Zachodzennie służy do podania wielkości tego obszaru.

• Menu Typ okna pozwala na wybór typu okna, przez które zostanie pomnożony wynik odwrotnej transformaty Fouriera. Pozwala to na polepszenie jakości otrzymanego filtru.

• W okienko Parametr okna wpisujemy współczynnik, potrzebny przy obliczaniu okien Kaisera i Czebyszewa.

• Naciśniecie klawisza Projektuj powoduje zaprojektowanie filtru zgodnie z podanymi w innych elementach okna wartościami.

• Zawartość okna edycyjnego zatytułowanego Czestotliwosc probkowania określa częstotliwość próbkowania używaną w systemie, w którym zostanie użyty filtr.

• Jeżeli zaznaczymy element Pamiętaj zostaną zapamiętane bieżące współczynniki filtru. Wszystkie charakterystyki będą wykreślane podwójnie, dopóki element nie zostanie wyłączony.

• Menu Typy analizy umożliwia użytkownikowi obejrzenie podstawowych charakterystyk zaprojektowanego filtru.

• Przycisk Koniec kończy pracę aplikacji.

Podobnie, jak w poprzednio opisanym programie menu górne zawiera funkcje przesyłające współczynniki filtru do przestrzeni roboczej MATLAB-a lub do programu dokonującego implementacji. Zawiera ono także pomoc do tej aplikacji.

Metody projektowania filtru przy użyciu kryterium błędu średniokwadratowego i kryterium minimalizacji błędu maksymalnego zostały połączone w jednej aplikacji, ze względu na podobną składnie funkcji MATLAB-a je realizujących.

Rys 5.5. Wygląd interfejsu do projektowania filtrów FIR metodą least quadred i Remeza.

Elementy umieszczone w oknie maja następujące przeznaczenie :

• W okienko edycyjne zatytułowane Dlugosc filtru wpisujemy rząd projektowanego filtru. Kiedy używamy metody Remeza istnieje możliwość oszacowania rzędu filtru, dla którego otrzymany filtr spełni nałożone na niego wymagania, co do kształtu charakterystyki amplitudowej.

• Pole edycyjne Czestotliwosci służy do podania wektora par częstotliwości (ponieważ definiują one pasma częstotliwości, na których określimy kształt charakterystyki amplitudowej ), które muszą mieścić się w przedziale od 0 do połowy częstotliwości próbkowania. Częstotliwości muszą zostać wpisane w porządku rosnącym.

• To co wpiszemy w okno Amplitudy zmienia swoje znaczenie w zależności od stanu znacznika Minimalny rzad. W przypadku wyłączonego znacznika, w tym miejscu wpisujemy wektor pożądanych amplitud charakterystyki amplitudowej w częstotliwościach podanych w okienku Czestotliwosci. Pożądana amplituda dla częstotliwości pomiędzy częstotliwościami (f(k),f(k+1)) dla k nieparzystego jest linią ciągłą łączącą punkty (f(k),a(k)) i (f(k+1),a(k+1)). Inaczej mówiąc metoda dla powyższej sytuacji stara się, jak najwierniej oddać zakładany kształt charakterystyki. Natomiast pożądany kształt charakterystyki pomiędzy częstotliwościami (f(k),f(k+1)) dla k parzystego jest nieokreślony. Te przedziały nazywamy pasmami przejściowymi lub pasmami "don't care". Dlatego liczba próbek amplitudy i częstotliwości musi być zgodna. Ponadto obie powinny być parzyste. W przypadku, gdy znacznik jest włączony częstotliwości podane w okienku Czestotliwosci oznaczają granice pasm, zaś wartości podane w okienku Amplitudy pożądane wartości amplitudy w zdefiniowanych pasmach. Wartości częstotliwości powinno być o dwie mniej niż wynosi podwojona liczba wartości amplitudy. Program projektuje filtry o parzystej długości z pasmem przepustowym umieszczonym połowy częstotliwości próbkowania. Jest bowiem niemożliwej zaprojektowanie filtru o liniowej fazie o nieparzystym rzędzie, kiedy charakterystyka częstotliwościowa dla połowy częstotliwości próbkowania równa jest zeru.

• Pole Wagi zmienia swoje znaczenie w zależności od ustawienia znacznika Minimalny rzad. Jeśli jest on wyłączony, to wagi określają jaka uwagę przywiązuje algorytm, aby osiągnąć zadaną amplitudę w określonym paśmie. Im wyższa waga tym, otrzymana charakterystyka bliżej aproksymuje idealną. Ponieważ pasm w których metoda stara się wyznaczyć charakterystykę jest o połowę mniej niż punktów częstotliwości, liczba wag też jest o połowę mniejsza. Natomiast, jeśli znacznik Minimalny rzad jest włączony, to w to pole wpisujemy maksymalną dopuszczalną odchyłkę od idealnej charakterystyki projektowanego filtru. Jeśli dotyczy ona pasma przepustowego to oznacza ona zafalowanie w tym paśmie, jeśli zaś odnosi się do pasma zaporowego to reprezentuje tłumienie. Inaczej niż w poprzednim wypadku liczba wag musi być równa liczbie amplitud, ponieważ teraz one wyznaczają pasma.

• Za pomocą menu Metoda wybieramy metodę projektowania. Do wyboru mamy metodę Remeza lub minimalizacji błędu średniokwadratowego.

• Jeśli znacznik Minimalny rzad jest zaznaczony, nie wpisujemy rzędu filtru, ale jest on automatycznie dobierany. Po naciśnięciu przycisku Projektuj w nieaktywnym okienku Rzad pojawi się wartość rzędu zaprojektowanego filtru. Opcja ta działa tylko dla metody Remeza, ale nie wtedy, gdy projektujemy filtr Hilberta lub dyferencjator.

• Menu Typ filtru pozwala na otrzymanie, oprócz typowego filtru, dwóch specjalnych filtrów: filtru Hilberta i dyferencjatora. Pierwszy z nich, w idealnym przypadku ma stałą charakterystykę amplitudową równa jedności, a jego faza zmienia się liniowo. Czyli taki filtr przepuszcza wszystkie składowe częstotliwościowe sygnału, zmieniając tylko ich fazę. Taki filtr jest wykorzystywany do zamiany sygnału rzeczywistego na sygnał zespolony. Natomiast dyferencjator jest odpowiednikiem układu różniczkującego. W dziedzinie czasu takiemu działaniu odpowiada pomnożenie widma sygnału przez urojoną funkcje narastającą liniowo. Dlatego cyfrowy układ różniczkujący ma liniowo narastającą charakterystykę amplitudową w swoim paśmie przepustowym, zachowując liniową fazę.

Pozostałe elementy okna maja takie samo znaczenie, jak dla poprzednich aplikacji. Współczynniki filtru można przesłać do przestrzeni roboczej lub do narzędzia dokonującego implementacji filtru. Użytkownik ma dostęp do pomocy, na którą składa się opis interfejsu i przykłady użycia.

Kolejnym interfejsem użytkownika jest implementacja wymuszonej metody minimalizacji błędu średniokwadratowego bez użycia pasm przejściowych. Przy użyciu tej metody stworzono dwie aplikacje, tu omówimy jedną z nich.

Rys 5.6. Wygląd okna metody wymuszonej minimalizacji błędu średniokwadratowego.

Opis elementów interfejsu :

• Wpisujemy w okno Dlugosc filtru rząd projektowanego filtru. Jeżeli projektujemy filtr posiadający pasmo przepustowe dla częstotliwości równej połowie częstotliwości próbkowania i wpiszemy w omawiane pole liczbę nieparzystą, to program po naciśnięciu klawisza Projektuj automatycznie zwiększy rząd filtru o jeden.

• Pole Czestotliwości służy do wprowadzenia wektora częstotliwości oddzielających kolejne pasma projektowanego filtru. Muszą one zawierać się pomiędzy zerem, a połową częstotliwości próbkowania. Częstotliwości muszą być tak podane, aby pokrywały całe pasmo częstotliwości. W szczególności w wektorze częstotliwości muszą pojawić się częstotliwości : zerowa i równa połowie częstotliwości próbkowania. Ponadto częstotliwości muszą być podane w porządku rosnącym.

• Element programu, nazwany Amplitudy, służy do wprowadzenia wektora pożądanych amplitud dla poszczególnych pasm. Ich liczba musi być o jeden mniejsza niż liczba częstotliwości podanych w okienku Czestotliwosci.

• Liczby wpisywane w pole Gorne granice określają górne granice, jakie może osiągnąć charakterystyka amplitudowa w konkretnym paśmie. Jeśli amplituda w paśmie wynosi 1, to podanie wartości granicy 1.05 oznacza, że charakterystyka w tym paśmie nie powinna przekroczyć wartości 1.05. Liczba podanych granic musi więc równać się liczbie podanych amplitud. Ponadto wartość granicy musi być większa od wartości pożądanej amplitudy podanej dla danego pasma.

• Okno edycyjne Dolne granice określa dolne granice, jakie może osiągnąć charakterystyka amplitudowa w konkretnym paśmie. Jeśli amplituda w paśmie wynosi 1, to podanie wartości granicy 0.95 oznacza, ze charakterystyka w tym paśmie nie przekroczy granicy 0.95. Liczba podanych dolnych granic musi więc równać się liczbie podanych amplitud. Ponadto wartość granicy musi być mniejsza od wartości amplitudy. Dla pasm zaporowych należy więc podawać liczby ujemne.

Pozostałe elementy okna są takie same, jak we wcześniej omówionych interfejsach GUI.

Uogólniona metoda Butterwortha pozwala na projektowanie zarówno filtrów IIR, jak i filtrów FIR. Niestety możemy projektować tylko filtry dolnoprzepustowe. Interfejs GUI dla tej metody jest przedstawiony na rysunku 5.7.

Rys. 5.7. Interfejs GUI uogólnionej metody Butterwotha.

Elementy okna posiadają następujące znaczenie:

• Liczba wpisana w okno edycyjne Rzad licznika określa rząd licznika transmitancji projektowanego filtru. Powinna być to dodatnia liczba całkowita. Przy projektowaniu filtrów typu FIR rząd licznika musi być parzysty. Jeśli tak nie jest program automatycznie zwiększa rząd licznika

• Wpisujemy w okienko Rzad mianownika rząd mianownika transmitancji projektowanego filtru. Powinna być to dodatnia liczba całkowita.

• W pole Czest. Odcięcia wpisujemy częstotliwość odcięcia projektowanego filtru dolnoprzepustowego. Musi ona zawierać się pomiędzy zerem, a połową częstotliwości próbkowania.

• Dzięki menu Typ filtru możemy projektować filtry IIR lub FIR. Należy pamiętać, że filtry typu FIR muszą mieć parzysty rząd.

• Jeśli znacznik Pokaz przebieg jest włączony, to po naciśnięciu przycisku Projektuj pojawia się dodatkowe okno. W tym oknie wykreślone zostają : amplitudowa charakterystyka filtru, jego opóźnienie grupowe oraz położenie zer i biegunów. Po naciśnięciu dowolnego klawisza okno zostaje zamknięte.

• Przycisk Estymuj pozwala na określenie rzędu filtru, którego granice pasm znajdują się w określonych częstotliwościach. Po naciśnięciu tego klawisza pojawia się okno, w którym musimy podać częstotliwość, przy której charakterystyka amplitudowa ma osiągać 0.95, oraz drugą częstotliwość, dla której charakterystyka opada poniżej wartości 0.05. Czyli pierwsza częstotliwość wyznacza granicę pasma przepustowego, a druga zaporowego.

Pozostałe elementy okna są podobne do omawianych wcześniej.

W pakiecie MATLAB została stworzona funkcja realizująca metodę Remeza, ale o znacznie rozszerzonych możliwościach.

Rys 5.8. Okno interfejsu zaawansowanej metody Remeza.

Oto opis elementów okna interfejsu GUI implementującego zaawansowana metodę Remeza :

• W okienko edycyjne Rzad filtru wpisujemy długość projektowanego filtru. Powinna być to dodatnia liczba całkowita. Należy zwrócić uwagę na parzystość rzędu filtru przy zmianie rodzaju projektowanego filtru!

• Okno Czestotliwosci służy do podania algorytmowi wektora par częstotliwości, które muszą się mieścić w przedziale od 0 do połowy częstotliwości próbkowania. Częstotliwości muszą zostać wpisane w porządku rosnącym. Podane pary częstotliwości określają pasma projektowanego filtru.

• Pole Amplitudy definiuje próbki charakterystyki amplitudowej określone w częstotliwościach podanych w polu opisanym powyżej . Liczba podanych amplitud musi być równa liczbie podanych częstotliwości. Pożądana kształt amplitudy dla częstotliwości pomiędzy częstotliwościami (f(k),f(k+1)) dla k nieparzystego jest linią ciągłą łączącą punkty (f(k),a(k)) i (f(k+1), a(k+1)). Inaczej mówiąc metoda dla powyższej sytuacji stara się jak najwierniej oddać zakładany kształt charakterystyki. Natomiast pożądany kształt charakterystyki pomiędzy częstotliwościami (f(k),f(k+1)) dla k parzystego jest nieokreślony. Te przedziały nazywamy pasmami przejściowymi lub pasmami "don't care". Istnieje jeden wyjątek od tej reguły. Jeżeli w menu Opcje algorytmu wybierzemy opcje Punkty czestotliwosci, to za pomocą znacznika 's' możemy tworzyć pasma składające się z jednego punktu. Z wyjątkiem tego przypadku liczba amplitud i częstotliwości powinna być parzysta.

• Znaczenie liczb wpisywanych w pole Wagi zależy od ustawienia menu Opcje pasm. Jeżeli wybraliśmy opcje Standardowa, liczby wpisane w to pole oznaczają wagi, przypisane do danego pasma. Im większa waga, tym różnica pomiędzy kształtem charakterystyki podanym w polu Ampiltudy, a otrzymaną charakterystyką jest mniejsza. Natomiast, jeśli wybraliśmy opcje Narzucone parametry liczba wpisana w to pole zależy od zawartości okienka Opcje pasm. Kiedy dane pasmo ma przypisaną do siebie literkę 'w', to liczba wpisane w pole Wagi oznacza wagę pasma. Z kolei, jeżeli znacznikiem jest litera 'c', to liczba wpisana w omawiane pole oznacza maksymalną odchyłkę od wartości wpisanej w pole Amplitudy. Czyli. jeśli jest to pasmo przepustowe liczba wpisana w pole Wagi oznacza maksymalne zafalowania zafalowania w tym paśmie, natomiast jeśli mamy do czynienia z pasmem zaporowym liczba określa minimalne tłumienie. Natomiast, jeśli wybrano opcje punkty czestotliwosci z menu Opcje algorytmu to okienko Wagi domyślnie jest nieaktywne. Jeśli zaznaczymy znacznik Niezalezne bledy to nasze okienko jest uaktywniane i spełnia taką samą funkcje jak dla opcji Standardowe menu Opcje algorytmu, czyli służy dla podania wag poszczególnych pasm.

• Menu Rodzaj symetrii zwiąże się z wyborem typu symetrii filtru. Z każdym typem symetrii wiąże się odpowiednia parzystość rzędu filtru oraz specyficzne wymagania na kształt charakterystyki częstotliwościowej. Wszystkie te własności są omówione we wstępie do rozdziału dotyczącego projektowania filtrów FIR. W tym programie można zaprojektować wszystkie typy filtrów FIR. Często, jeśli charakterystyka filtru odbiega od zadanej przyczyną jest niewłaściwy typ symetrii.

• Dzięki menu Typ filtru możemy tworzyć filtry minimalno- i maksymalnofazowe. Filtr minimalnofazowy to taki, którego wszystkie zera znajdują się wewnątrz okręgu jednostkowego. Natomiast wszystkie zera filtru maksymalnofazowego znajdują się na zewnątrz okręgu jednostkowego. Niestety te filtry tracą właściwość liniowej fazy.

• Menu Opcje algorytmu służy do zmieniania właściwości algorytmu. Zawiera trzy opcje :

1. Standardowa - algorytm działa jak algorytm Remeza

2. Narzucone parametry - algorytm tak projektuje filtr, aby pewne pasma spełniały narzucone przez projektanta parametry np. tłumienie

3. Punkty częstotliwości - dzięki tej opcji możemy nadawać pewnym częstotliwościom specjalne właściwości.

• Okno edycyjne Opcje pasm służy do wprowadzania znaczników zmieniających własności pasm. Znaczniki muszą być otoczone poprzez apostrofy. Ich znaczenie zależy od ustawienia menu Opcje algorytmu:

• przy opcji Standardowy okienko jest nieaktywne

• wybór opcji Narzucone parametry powoduje, ze znaczniki wpisane w to okno mówią nam, w jaki sposób traktowane są liczby wpisane w okno Wagi. Znacznik 'w' powoduje, że to pasmo jest pasmem, w którym błąd będzie minimalizowany w zależności od jego wagi. Z kolei znacznik 'c' powoduje, że charakterystyka będzie tak projektowana, aby spełniła narzucone na nią wymagania podane w oknie Wagi. Uwaga! Przynajmniej jedno pasmo filtru musi być zaznaczone jako ważone

• opcja punkty czestotliwosci służy do nadania konkretnym częstotliwościom specjalnych właściwości. Właściwości te nadaje się poprzez podanie przydzielenie każdej z częstotliwości podanej w polu Czestotliwosci jednego z czterech znaczników :

1. 'n' - oznacza zwykły punkt częstotliwości

2. 's' - tworzy jednopunktowe pasmo częstotliwości. Dla częstotliwości, przy której postawiono ten znacznik, musi w wektorze amplitud znajdować się odpowiednia wartość. Ten znacznik wprowadza odstępstwo od zasady, że ilość wprowadzonych częstotliwości musi być parzysta. Znacznik ten powinien znajdować się pomiędzy dwoma pasmami.

3. 'f' - ten znacznik wymusza wartość charakterystyki wewnątrz pewnego pasma. Znacznik musi być poprzedzony znacznikiem 'i'.

4. 'i' - oznacza nieokreślony punkt częstotliwości. Służy do dzielenia pasma na dwa przylegające. Nie może być używany w paśmie przejściowym.

• Jeżeli z menu Opcje algorytmu wybrano opcję standardowe lub punkty częstotliwości, to możemy użyć znacznika Niezalezne bledy. Bez tego znacznika minimalizowany jest całkowity błąd charakterystyki częstotliwościowej w całym zakresie częstotliwości. Dzięki temu znacznikowi dla każdego pasma używany jest odrębny błąd, który podlega odrębnej minimalizacji.

Pozostałe elementy okna maja takie same właściwości jak w interfejsach omówionych wcześniej.

Rys 5.9. Interfejs GUI przeznaczony do projektowania filtrów FIR o paśmie przejściowym w kształcie podniesionego kosinusa.

Na rysunku 5.9 przedstawiono wygląd interfejsu GUI projektującego filtru FIR o paśmie przejściowym w kształcie podniesionej kosinusoidy. Przy pomocy tego narzędzia można projektować tylko filtry dolnoprzepustowe. Poniższe wyliczenie opisuje charakterystyczne elementy tego interfejsu GUI :

• Wpisujemy w okienko Rzad filtru długość projektowanego filtru. Po wpisaniu wartości w to okno zostaje automatycznie zmodyfikowana zawartość okna Opoznienie.

• W pole Czest. Odciecia wpisujemy częstotliwość odcięcia ω₀ projektowanego filtru. Musi ona zawierać się pomiędzy zerem, a połową częstotliwości próbkowania.

• Szerokość pasma Δω, podawana w polu Szerokosc pasma, wyznacza jak szerokie pasmo przejściowe będzie posiadał projektowany filtr. Szerokość pasma przejściowego nie powinna być zbyt wielka. Musi spełniać następujące warunki:

• ω₀ - Δω > 0,

• ω₀ + Δω < fs/2.

• Element okna, nazwany Opoznienie, służy do wprowadzenia wartości opóźnienia delay. Parametr ten mówi po ilu okresach częstotliwości próbkowania na wyjściu filtru pojawi się jego odpowiedź. Opóźnienie musi być liczbą całkowitą. Powinna spelniąć warunki :

• delay > = 0,

• delay < = (rzad filtru + 1) .

Domyślnymi wartościami opóźnienia są :

• rząd filtru podzielony przez dwa dla rzędu parzystego,

• (rząd filtru plus jeden) podzielony przez dwa dla rzędu nieparzystego.

• Dzięki menu Ksztalt przejścia może wybrać kształt pasma przejściowego. Pasmo może mieć kształt pierwiastkowanej lub niezniekształconej podniesionej kosinusoidy.

• Za pomocą menu Typ okna wybieramy typ okna użytego do projektowania filtru. Od wyboru okna zależy tłumienie filtru w paśmie zaporowym oraz szerokość pasma przejściowego. Przy oknach Kaisera i Czebyszewa należy podać współczynnik okna.
  Uwaga! Należy zachować ostrożność używając niestandardowego opóźnienia i funkcji okna różnej od prostokątnej.

• W okno edycyjne Parametr okna wpisujemy współczynnik, potrzebny przy obliczaniu okien Kaisera i Czebyszewa.

Następnym interfejsem GUI, którego wygląd jest przedstawiony na rysunku 5.10, jest narzędzie do projektowania filtrów IIR przy pomocy prototypu analogowego.

Rys 5.10. Wygląd okna interfejsu GUI projektującego filtry IIR metodą prototypu analogowego.

Opis najważniejszych elementów składowych tego interfejsu :

• Włączenie znacznika Minimalny rzad powoduje aproksymacje minimalnego rzędu filtru, przy którym dany projekt spełni nasze wymagania na kształt charakterystyki amplitudowej. Zmiana ustawienia tego przełącznika powoduje zmianę działania pól, do których wpisujemy ramy częstotliwościowe filtru.

• W pole edycyjne Rzad filtru wpisujemy liczbę oznaczającą, jakiego rzędu filtr chcecmy otrzymać. Jeśli jest włączana opcja aproksymacji minimalnego rzędu, to okno jest nieaktywne. Wyświetlany jest tylko otrzymany w wyniku aproksymacji rząd filtru.

• Na ten element Czest. pasma przepustowego składają się dwa okna, w które możemy wpisywać liczby, które mają sens częstotliwości. Ich znaczenie zależy od ustawienia przełącznika Minimalny rzad .Jeśli przełącznik jest wyłączony w okienka wpisujemy częstotliwość(i) odcięcia filtru jako całości, bez dokładnej określania pasm, Natomiast, jeśli przełącznik jest włączony, to w te pola wpisujemy częstotliwości graniczne pasm(a) przepustowego. Tak więc, jeśli projektujemy filtry górno- lub dolnoprzepustowy aktywne jest tylko jedno okno, zaś jeśli chcemy otrzymać filtr pasmowoprzepustowy lub pasmowozaporowy to aktywne są oba okna.

• Okienka pasma zaporowego są aktywne tylko, jeśli przełącznik Minimalny rzad jest aktywny. Służą one do wpisania w nie częstotliwości pasm zaporowych filtru.

• W pole edycyjne Zafalowania wpisujemy zadaną wielkość zafalowań (w dB) w paśmie przepustowym. Okno aktywne, jeśli :

• Włączony jest przełącznik Minimalny rzad

• Używamy do projektowania metody Czebyszewa 1 lub metody eliptyczna

• W okno Tlumienie wpisujemy zadaną wielkość tłumienia (w dB) w paśmie zaporowym. Okno aktywne jeśli :

• Włączony jest przełącznik Minimalny rzad,

• Używamy do projektowania prototypu uzyskanego metodą Czebyszewa 2 lub metodą eliptyczna.

• Menu Os czestotliwosci określa, w jakich jednostkach są podawane częstotliwości

• Unormowane względem połowy częstotliwości próbkowania.

• Podawane są w Hz.

• Podawane są w kHz.

Jeśli oś nie jest unormowana, to wymagane jest podanie częstotliwości próbkowania

• Służy do wyboru prototypu analogowego, który zostanie użyty do projektowania filtru. Można użyć czterech prototypów : Butterwortha, Czebyszewa typ pierwszy, Czebyszewa typ drugi oraz filtru eliptycznego.

• Przy pomocy menu Typ filtru wybieramy typ filtru spośród czterech standardowych możliwości: dolnoprzepustowego, górnoprzepustowego, pasmoprzepustowego i pasmowozaporowego.

Kolejnym programem do projektowania filtrów IIR jest aplikacja implementująca metodę Yule-Walkera.

Rys 5.11. Wygląd okna metody projektującej filtry IIR metodą Yule-Walkera.

Elementy składowe okna spełniają następujące funkcje :

• W okno Rzad filtru wpisujemy długość projektowanego filtru. Powinna być to dodatnia liczbą całkowitą.

• Pole edycyjne Czestotliwosci słuzy do podania wektora par częstotliwości, które muszą mieścić się w przedziale od 0 do połowy częstotliwości próbkowania. W programie posługujemy się częstotliwościami unormowanymi do połowy częstotliwości próbkowania. Częstotliwości muszą zostać wpisane w porządku niemalejącym. Częstotliwości mogą się powtarzać, jeśli w tym miejscu charakterystyki następuje styk dwóch pasm.

• Dzięki polu Amplitudy możemy podać próbki pożądanej charakterystyki amplitudowej w częstotliwościach podanych w polu opisanym powyżej. Ich liczba musi być równa liczbie częstotliwości podanym w okienku Czestotliwosci.

Pozostałe elementy okna pełnią takie same funkcje, jak we wszystkich opisanych wcześniej aplikacjach.

Następną aplikacją jest interfejs implementujący metodę projektowania filtrów IIR przy pomocy kryterium minimalizacji p.-tej normy.

Rys 5.12. Okno metody kryterium minimalizacji p.-tej normy.

Elementy składowe interfejsu GUI dla omawianej metody mają następujące przeznaczenie :

• W okno edycyjne Rzad licznika wpisujemy zadany rząd licznika transmitancji projektowanego filtru. Powinna być to dodatnia liczba całkowita. Jeżeli projektujemy filtr o zadanym kształcie opóźnienia grupowego lub filtr do kompensacji charakterystyki fazowej innego filtru, to omawiane okno jest nieaktywne, ponieważ taki filtr posiada wyłącznie bieguny.

• Pole edycyjne Rzad mianownika zadaje algorytmowi żądany rząd mianownika transmitancji projektowanego filtru. Powinna być to dodatnia liczba całkowita.

• Okno Wektor czestotliwosci służy do podania wektora częstotliwości, które muszą mieścić się w przedziale od 0 do połowy częstotliwości próbkowania. Częstotliwości muszą zostać wpisane w porządku rosnącym.

• W pole zatytułowane Amplitudy wpisujemy wektor amplitud, jakie chcemy, żeby osiągnęła charakterystyka filtru w częstotliwościach podanych w okienku Czestotliwosci. Dlatego też, liczba amplitud musi być równa liczbie częstotliwości. Jeżeli z menu Metoda wybrano Opoznienie podane w tym okienku wartości mają sens opóźnienia grupowego podanego w próbkach.

• Element okna nazwany Wektor czest. granicznych służy do podania algorytmowi informacji o pasmach, jakie ma posiadać projektowany filtr. Częstotliwości podane w tym polu muszą znajdować także w oknie Czestotliwosci. Ilość częstotliwości podanych w omawianym oknie może być mniejsza niż podana w oknie Czestotliwosci. Dzięki takiemu podwojeniu możliwe jest uzyskanie takiego kształtu charakterystyki w danym paśmie, jakie chcemy otrzymać. Kiedy z menu Metoda wybierzemy opcje kompensacja znaczenie liczb wpisywanych w to pole zmienia się. W tym przypadku w to pole należy wpisać częstotliwości graniczne pasm przepustowych filtru, w których chcemy uzyskać płaską charakterystykę opóźnienia grupowego. Jeśli filtr, którego charakterystykę chcemy zmienić, jest filtrem dolnoprzepustowym o granicy pasma przepustowego w częstotliwości 0.4, to w omawiane okienko należy wpisać wektor 0 0.4.

• Liczby wpisane w okno Wagi oznaczają, jaką wagę przywiązuje algorytm do osiągnięcia pożądanej wartości charakterystyki amplitudowej lub opóźnienia grupowego w danym punkcie częstotliwości. Im większa waga, tym różnica pomiędzy zadaną przez projektanta wartością amplitudy (opóźnienia), a otrzymaną w wyniku projektowania będzie mniejsza. Należy pamiętać, że waga jest przypisana do danego punktu częstotliwości, co z kolei powoduje, ze podana liczba wag musi być równa ilości częstotliwości wpisanych w pole Czestotliwosci

• W pole nazwane Zakres norm wpisujemy dwie liczby. Oznaczają one minimalny i maksymalny rząd normy użyty podczas projektowania filtru. Podane liczby muszą być liczbami parzystymi.

• Zawartość pola Gestosc określa liczbę odpowiedzalną za to, ile próbek częstotliwości zostanie użytych podczas projektowania filtru. Im ich gęstość będzie większa, tym odtrzymana charakterystyka lepiej przybliży zadaną przez projektanta, kosztem zwiększonej liczby obliczeń.

• Dzięki menu Metoda możemy wybrać metodę projektowania filtru. Mamy do wyboru cztery metody :

• metodę optymalizacji kształtu charakterystyki amplitudowej w sensie minimalizacji normy rzędu p.

• metodę optymalizacji kształtu charakterystyki amplitudowej w sensie minimalizacji normy rzędu p z dodatkową możliwością regulacji położenia biegunów.

• metodę umożliwiającą projektowanie filtrów o określonym kształcie charakterystyki opóźnienia grupowego,

• stworzenie wszechprzepustowego filtru kompensującego charakterystykę opóźnienia grupowego poprzedzającego go konwencjonalnego filtru.

• Pole Radiks jest aktywne tylko, jeżeli z menu Wybor metody wybrano algorytmy nazwane L-norma 2 i opoznienie. Dzięki temu okienku możemy regulować maksymalną odległość położenia biegunów od początku układu współrzędnych. Wartość wpisana w to pole powinna zawierać się w zakresie od 0 do 1. Przy pomocy tej opcji możemy osiągnąć lepsze zachowanie parametrów filtru, kiedy dokonamy kwantyzacji współczynników.

W górnym menu znajdują się następujące funkcje:

• Menu Filtr zawiera trzy składniki :

• Funkcja Import służącą do pobrania z przestrzeni roboczej współczynników filtru, które posłużą za wartości początkowe dla algorytmu.

• Funkcja Eksport służy do wysłania do przestrzeni roboczej współczynników zaprojektowanego filtru.

• Funkcja Koniec kończy pracę programu.

• Menu Pomoc także zawiera trzy funkcje :

• Funkcja Opis interfejsu, wyświetlająca opis interfejsu programu,

• Funkcje Przyklady, która zawiera kilka przykładów użycia programu,

• Funkcja O programie wyświetlająca krótką informację o programie.

Dla metod realizujących projektowanie filtrów poprzez modelowanie stworzono interfejs GUI przedstawiony na rysunku 5.13.

Rys. 5.13 Wygląd interfejsu metod projektujących filtr IIR przez modelowanie

Opis interfejsu :

• Przycisk Import służy do importowania współczynników filtru lub jego odpowiedniej charakterystyki z przestrzeni roboczej Matlab'a. Dzięki niemu możemy importować do programu :

• Współczynniki filtru. Wtedy możemy używać wszystkich metod projektowania.

• Zespoloną charakterystykę filtru wraz z odpowiadającym mu wektorem częstotliwości. W tym przypadku dostępna jest tylko metoda Invfreqz.

• Odpowiedź impulsową filtru. Jeśli wybierzemy tą charakterystykę dostępne będą dwie metody : Promy'ego i Stmcb.

• Odpowiedz filtru na dowolne pobudzenie rzeczywiste. Po pobraniu tych dwóch przebiegów dostępna będzie tylko metoda Stmcb.

• Okienko Rzad licznika wymusza rząd licznika transmitancji projektowanego filtru. Powinna być to dodatnia liczba całkowita.

• Okienko Rzad mianownika wymusza rząd mianownika transmitancji projektowanego filtru. Powinna być to dodatnia liczba całkowita.

• Przycisk Invfreqz projektuje filtr przekształcając zespoloną charakterystykę częstotliwościowa filtru na jego transmitancję.

• Przycisk Prony projektuje filtr przekształcając odpowiedź impulsową filtru na jego transmitancje przy pomocy metody Prony'ego.

• Przycisk Stmcb projektuje filtr przekształcając odpowiedź filtru na dane pobudzenie na współczynniki filtru przy pomocy metody Steiglitz-McBride'a

Wygląd okna programu do sprawdzania mieszczenia się charakterystyk różnych filtrów w tych samych ramach projektowych jest przedstawiony na rysunku umieszczonym poniżej.

Rys 5.14. Wygląd narzędzia do nakładania ram projektowych na otrzymane charakterystyki fitru.

Poprzez menu Typ filtru określamy jaki typ filtru chcemy zaimportować do układu. Przycisk Import służy do importu współczynników filtru z przestrzeni roboczej Matlab'a. Po naciśnięciu przycisku pojawi się okienko, w które należy wpisać nazwę (w przypadku filtru FIR) lub nazwy (w przypadku filtru IIR) zmiennych zawierających licznik lub mianownik transmitancji. Po zaimportowaniu współczynników wykreślana jest charakterystyka amplitudowa filtru. Można importować dowolną ilość filtrów, przy czym ich charakterystyki będą wykreślane na jednym wykresie, przez co może ich przebieg może nie być czytelny. W oknie edycyjnym Czest. pasm przepustowych wpisujemy wektor par częstotliwości, które wyznaczają granice pasm przepustowych badanych filtrów. Częstotliwości te muszą mieścić się w przedziale od 0 do połowy częstotliwości próbkowania. Liczba wpisanych częstotliwości musi być parzysta, ponieważ każde pasmo ma swoją górną i dolna granicę. Natomiast pole Czest. pasm zaporowych służy do podania wektora par częstotliwości, które wyznaczają pasma zaporowe badanego filtru. Częstotliwości te muszą mieścić się w przedziale od 0 do połowy częstotliwości próbkowania. Liczba wpisanych częstotliwości musi być parzysta, ponieważ każde pasmo ma swoją górna i dolną granicę. Okno Zafalowania pozwala na podanie maksymalnych dopuszczalnych zafalowań w każdym podanym paśmie przepustowym badanego filtru. Liczba podanych zafalowań musi być o połowę mniejsza od liczby podanych częstotliwości pasm przepustowych, gdyż należy zdefiniować ten parametr w każdym podanym paśmie. Zafalowania należy podawać w dB. Natomiast okno Tlumienie określa minimalne tłumienie, podane w mierze decybelowej, w każdym podanym paśmie zaporowym badanego filtru. Liczba podanych tłumień musi być o połowę mniejsza od liczby podanych częstotliwości pasm zaporowych, gdyż każde tłumienie przypisane jest do jednego pasma. W okienku Czest. probkowania podajemy częstotliwość próbkowania używaną przez badane filtry. Należy pamiętać, że częstotliwości używane przy określaniu granic pasm muszą być mniejsze od połowy częstotliwości próbkowania. Przycisk Pokaz nanosi ramy projektowe na charakterystykę filtru.

W górnym menu znajdują się następujące funkcje:

• Menu Obsluga zawiera dwa składniki :

• Funkcja Import służącą do pobrania z przestrzeni roboczej Matlab'a współczynników filtru.

• Funkcja Koniec kończy pracę programu.

• Menu Zblizenia posiada trzy funkcje :

• Funkcje Pasmo przepustowe, dzięki której możemy obejrzeć w powiększeniu wybrane pasmo przepustowe.

• Funkcje Pasmo zaporowe, dzięki któremu możemy obejrzeć w powiększeniu wybrane pasmo zaporowe.

• Funkcje Calosc wyświetlająca całą charakterystykę z naniesionymi ramami projektowymi.

• Menu Pomoc także zawiera dwie funkcje :

• Funkcja Opis interfejsu, wyświetlająca opis interfejsu programu.

• Funkcja O programie wyświetlająca krótką informację o programie.

Narzędzie od dokonania implementacji i zbadania wpływu kwantyzacji na charakterystyki filtrów zostało zrealizowane w postaci interfejsu GUI przedstawionego na rysunku 5.15

Rys 5.15. Wygląd okna programu implementacji filtrów.

Narzędzie dokonuje implementacji filtrów cyfrowych, zarówno typu IIR, jak i filtrów typu FIR. Wybór, który typ filtru cyfrowego jest aktualnie implementowany jest dokonywany za pomocą menu Typ filtru. Import współczynników filtru do programu jest dokonywany poprzez naciśnięcie klawisza Import. Wówczas pojawia się okno dialogowe, w które należy wpisać nazwę zmiennej lub zmiennych zawierających współczynniki filtru. Po dokonaniu importu współczynników może dokonać implementacji filtru poprzez naciśniecie klawisza Implementuj. Strukturę filtru wybiera się przy pomocy dwóch menu rozwijanych objętych wspólnym tytułem Strukura. W programie są dostępne następujące struktury filtrów cyfrowych :

1. dla filtrów o skończonej odpowiedzi impulsowej:

1. struktura bezpośrednia,

2. struktura transponowana,

3. struktura dla filtrów o symetrycznych współczynnikach,

4. struktura dla filtrów o antysymetrycznych współczynnikach,

5. struktura kratowa,

6. struktura kaskadowa.

2. dla filtrów o nieskończonej odpowiedzi impulsowej:

1. struktura bezpośrednia pierwsza,

2. struktura bezpośrednia druga,

3. struktura transponowana pierwsza,

4. struktura transponowana druga,

5. struktura kratowa,

6. struktura zmiennych stanu,

7. struktura równoległa,

8. struktura kaskadowa.

Po uruchomieniu implementacji są sczytywane ustawienia dotyczące kwantyzacji filtru. Jej efekty na działanie filtru jest przeprowadzony na następujących parametrach filtru i elementach składowych struktury filtru:

1. współczynnikach filtru,

2. danych wejściowych,

3. danych wyjściowych,

4. na danych doprowadzanych na wejście mnożnika,

5. na wyjściu mnożnika,

6. na wyjściu sumatora.

W każdym z podanych wyżej miejscu wstawiany jest kwantyzator. Posiada on szereg następujących parametrów :

1. Sposób przedstawiania liczb binarnych. Menu zmieniające tą wartość znajdują się w kolumnie zatytułowanej Format :

1. format stałoprzecinkowy o długości słowa od 2 do 53 bitów, przy czym część ułamkowa liczby może zmieniać od zera do wartości o jeden mniejszej niż długość słowa (opcja stalo),

2. dowolny format zmiennoprzecinkowy o długości słowa od 2 do 64 bitów, przy czym eksponenta może zmieniać się od 0 do 11 bitów (opcja zmienno),

3. liczb zmiennoprzecinkowych pojedynczej precyzji zgodnych z normą IEEE (opcja single),

4. liczb zmiennoprzecinkowych podwójnej precyzji zgodnych z normą IEEE (opcja double).

2. Sposób zaokrąglania używany dla kwantyzacji wartości liczbowych. Menu obsługujące tę opcje znajduje się w kolumnie z nagłówkiem Zaokrąglanie :

1. zaokrąglanie do najbliższej dostępnej kwantyzowanej wartości. Jeżeli zaokrąglana liczba znajduje się dokładnie po środku dwóch najbliżej dostępnych wartości, to jest zaokrąglana w górę (opcja standardowe),

2. zaokrąglanie do najbliższej większej skwantyzowanej wartości (opcja w gore),

3. zaokrąglanie do najbliższej mniejszej skwantyzowanej wartości (opcja w dol),

4. zaokrąglanie do najbliższej większej skwantyzowanej wartości dla liczb ujemnych, zaś dla liczb dodatnich stosuje się zaokrąglanie do najbliższej mniejszej skawntyzowanej wartości (opcja mieszane),

5. zaokrąglanie do najbliższej dostępnej kwantyzowanej wartości. Jeżeli zaokrąglana liczba znajduje się dokładnie po środku dwóch najbliżej dostępnych wartości, to jest zaokrąglana w górę tylko wtedy, gdy najmniej znaczący bit po zaokrągleniu przyjmie wartość 1 (opcja zbiezne).

3. sposób obsługiwania nadmiaru w arytmetyce stałoprzecinkowej. Tę właściwość obsługują menu umieszczone w kolumnie z tytułem Nadmiar

1. Jeżeli liczba, która ma zostać skwantyzowana leży poza obszarem dostępnych wartości, to jest kwantyzowana odpowiednio do największe lub najmniejszej wartości (opcja nasycenie),

2. Jeżeli liczba, która ma zostać skwantyzowana leży poza obszarem dostępnych wartości, to jest przenoszona z powrotem do zakresu dostępnych wartości zgodnie z arytmetyka liczb w kodzie uzupełnień do dwóch (opcja zawijanie).

4. Długość słowa cyfrowego używanego przez kwantyzator. Dla wszystkich formatów składa się ona z dwóch liczb. Pierwsza liczba niezależnie od formatu oznacza długość słowa. Natomiast druga liczba dla formatu stałoprzecinkowego oznacza długość części ułamkowej, zaś dla formatu zmiennoprzecinkowego oznacza długość eksponenty. Menu obsługujące tę opcję znajdują się w kolumnie zatytułowanej Dlugosc.

Dla formatów zmiennoprzecinkowych zgodnych z normami IEEE ustawienie trzech ostatnich właściwości jest niemożliwe, ponieważ są one ustalone przez właściwą normę.

Ponadto kwantyzatory różnią się do siebie interpretowaniem wartości 1. Jeden z nich przyjmuje, że wartość jeden nie podlega kwantyzacji (opcja unitquanizer), zaś drugi dokonuje jej kwantyzacji (opcja unitquanizer). Te możliwości prezentują menu umieszczone w kolumnie opatrzonej tytułem Kwantyzator.

W górnym menu znajdują się następujące funkcje:

1. W menu Obsługa :

1. Funkcja Import importuje współczynniki filtru z przestrzeni roboczej.

2. Funkcja Eksport eksportuje współczynniki skwantyzowanego filtru do przestrzeni roboczej.

3. Funkcja Plik nagłowkowym generuje plik nagłówkowy dla filtru zaprojektowanego w CCS. Plik jest generowany tylko dla filtrów stałoprzecinkowych i filtrów używających pojedynczej precyzji.

4. Funkcja Implementuj dokonuje implementacji.

5. Funkcja Koniec kończy pracę.

2. W menu Filtr :

1. Funkcja Format zmienia sposób reprezentacji liczb binarnych dla całego filtru.

2. Funkcja Zaokraglanie zmienia sposób zaokrąglania dla całego filtru.

3. Funkcja Nadmiar zmienia sposób obsługi nadmiaru dla całego filtru.

4. Funkcja Slowo zmienia długość słowa dla całego filtru.

5. Funkcja Skalowanie dokonuje skalowania filtru IIR zaimplementowanego w postaci kaskadowej.

3. W menu Analizy :

1. Zbiór funkcji Amplitudowa rysuje charakterystykę amplitudową filtru przed i po kwantyzacji w trzech skalach :

1. Logarytmicznej,

2. Liniowej,

3. Kwadratowej.

2. Funkcja Fazowa rysuje charakterystykę fazową filtru przed i po kwantyzacji.

3. Funkcja Bieguny pokazuje położenie zer i biegunów filtru przed i po kwantyzacji.

4. Funkcja Opoznienie grupowe rysuje charakterystykę opóźnienia grupowego filtru przed i po kwantyzacji.

5. Funkcja Impuls jednostkowy rysuje filtru przed i po kwantyzacji na pobudzenie impulsem jednostkowym.

6. Funkcja Skok jednostkowy rysuje odpowiedź filtru przed i po kwantyzacji na pobudzenie skokiem jednostkowym.

7. Funkcja Cykle graniczne sprawdza, czy filtr generuje cykle graniczne.

8. Funkcja Analiza szumowa dokonuje analizy jakości filtru zgodnie z metodą podana w rozdziale 3.

9. Funkcja Stabilnosc sprawdza stabilność skwantyzowanego filtru.

4. Menu Pomoc zawiera pomoc do tego interfejsu GUI.




5.2 Filtry cyfrowe zaprojektowane przy pomocy Code Composer Studio

Code Composer Studio CCS jest narzędziem firmy Texas Instrments służący do generacji kodu na procesory tej firmy. W tej pracy używany CCS dla rodziny procesorów C6000. Ma ona szczególnie duże możliwości z powodu równoległej architektury procesora, która umożliwia wykonywanie do ośmiu operacji równocześnie. Dlatego napisanie efektywnego kodu jest dla programisty dosyć trudne. Dlatego stworzono narzędzia umożliwiające automatyczną optymalizację kodu napisanego w języku C lub w uproszczonym asemblerze o nazwie Linear Asembler, którym pozwala zautomatyzować zagadnienia dotyczące odpowiedniego przydziału rejestrów oraz podziału rozkazów na tak zwane pakiety wykonawcze, to znaczy na grupy rozkazów, które mają zostać wykonane jednocześnie. Według twórców CCS wykonane przez nich narzędzia optymalizacji kodu pozwalają na uzyskanie kodu, którego szybkość wykonania jest niewiele mniejsza od optymalnego kodu napisanego w asemblerze. W celu maksymalizacji efektywności kodu kompilator składa tak wiele instrukcji w równoległe pakiety, jak to tylko możliwe. Aby tego dokonać musi określić zależność pomiędzy kolejnymi instrukcjami. Termin zależność w praktyce oznacza, że jedna instrukcja musi być wykonana przed innymi, w celu zachowania poprawnego działania programu, na przykład zmienna musi być pobrana z pamięci zanim zostanie wykorzystana. Brak zależności pomiędzy instrukcjami oznacza możliwość wykonania ich jednocześnie. Jeśli kompilator nie jest w stanie stwierdzić, czy określone instrukcje a i b są niezależne, to zakłada, że są one od siebie zależne, co prowadzi do tego, że instrukcja b zostanie wykonana dopiero wtedy, gdy swoje działanie zakończy instrukcja a. Dla kompilatora często jest trudnym zadaniem określić czy instrukcje, które wymagają dostępu do pamięci, są od siebie niezależne. Dlatego użytkownik powinien oznaczyć wskaźniki do pamięci słowem kluczowym restict, co oznacza, że w wywołaniu funkcji tylko ten wskaźnik będzie operował na danym zbiorze danych. Kolejnym sposobem rozwiązania tego problemu jest użycie opcji kompilatora -mt, która stara się rozwiązać problemy związane z istniejącymi w programie zależnościami pomiędzy instrukcjami. Podczas procesu optymalizacji dokonywane są następujące najważniejsze czynności:

1. przeprowadzane jest upraszczanie grafu przepływu programu;

2. optymalizowane jest rozmieszczenie zmiennych w rejestrach, aby sięganie do pamięci było jak najrzadsze;

3. usuwa nadmiarowe fragmenty kodu takie jak: nieużywane funkcje, nieużywane przypisania;

4. optymalizuje wykonanie pętli;

5. rozwija pętle, aby przyśpieszyć działanie programu;

6. uprasza wyrażenia, aby ich wykonanie było jak najszybsze, na przykład wyrażenie c = (b+4) - (c-1) jest upraszczane do postaci c  = b - c  + 3;

7. optymalizuje odwołania do pamięci,

8. jeżeli stwierdzi odwołanie do funkcji, której wielkość jest mała, to zastępuje odwołanie do funkcji poprzez wstawienie jej kodu w odpowiednie miejsce; w ten sposób zmniejszony jest rozmiar programu i zwiększana jego szybkość działania;

W procesie optymalizacji najważniejszy problemem jest taka optymalizacja wykonywania pętli, aby w jak największym stopniu wykorzystać moc procesora. Równoległa struktura procesorów rodziny 6000 pozwala na rozpoczęcie nowej iteracji pętli, zanim poprzednia zakończyła swój przebieg. Celem optymalizacji jest jak najwcześniejsze rozpoczęcie wykonywania nowej iteracji pętli. Możliwość przekształcenia oraz uzyskaną szybkość wykonania pętli zależy od minimalnego przedziału iteracji, który jest zdefiniowany jako minimalna liczba cykli, które procesor musi odczekać pomiędzy rozpoczęciem wykonywania kolejnych iteracji pętli równolegle. Im ten przedział jest mniejszy, tym szybciej zostanie wykonana pętla. Liczba zasobów procesora oraz istnienie zależności pomiędzy kolejnymi iteracjami pętli ogranicza minimalny przedział pomiędzy iteracjami. Głównym ograniczeniem jest liczba dostępnych zasobów procesora: ilość rejestrów, możliwość dokonania odczytu i zapisu do pamięci. Dlatego kompilator stara się wygenerować kod, w którym nie występują konflikty w dostępie do jednostek funkcjonalnych procesora pomiędzy wykonywanymi równolegle iteracjami. Ważnym problemem jest także wybór momentu, w którym program dokonuje obliczania wartości pośrednich. Kod programu musi zostać ułożony w taki sposób, aby iteracje pracujące równolegle nie zakłócały się wzajemnie. Kolejnym zagadnieniem związanym z wykonywanie instrukcji kolejnych iteracji w sposób równoległy jest fakt, że często wynik bieżącej iteracji jest wykorzystywany przez wszystkie przyszłe iteracje pętli. Istnienie takich ścieżek przepływu danych pomiędzy kolejnymi iteracjami może silnie wpływać na minimalny przedział iteracji, czasami może nawet uniemożliwić wykonywanie ich w sposób równoległy. Autorzy kompilatora języka podają następujące zalecenia dotyczące pisania pętli w języku C, aby zwiększyć ich efektywność:

• nie należy używać złożonych instrukcji warunkowych wewnątrz pętli,

• nie powinno dokonywać wywoływania funkcji wewnątrz pętli,

• w pętlach nie powinna pojawiać się instrukcja break, która przerywa wykonywanie pętli,

• pętle powinny mieć zmniejszający się licznik pętli,

• licznik pętli nie powinien być modyfikowany wewnątrz pętli,

• kod wewnątrz pętli nie powinien być zbyt duży, ponieważ może zabraknąć zasobów procesora i pętla nie zostanie wykonywania równolegle.

W pakiecie CCS istnieje wolnostojący symulator load6x, który umożliwia uzyskanie informacji o liczbie cykli trwania wszystkich funkcji zawartych w badanym programie. Oto fragment wydruku uzyskany przy pomocy tego symulatora:

Area Name Count Inclusive Incl-Max Exclusive Excl-Max

CF main() 1 349862 349862 4845 4845

CF fir() 128 343040 2680 343040 2680

Kolumna nazwana Area Name zawiera nazwy wszystkich funkcji zawartych w programie. Skrót CF oznacza funkcję napisaną w języku C. Natomiast kolumna nazwana Count pokazuje ilość wywołań danej funkcji w programie. Z kolei w kolumnie opatrzonej tytułem Inclusive zawarta jest informacja o całkowitej liczbie cykli, jaką zużytkował procesor wykonując daną funkcję. Kolumna Incl-Max pokazuje najdłuższy czas wykonywania danej funkcji. Natomiast kolumny Exclusiv i Excl-Max mają podobne znaczenie jak dwie poprzednie, ale od czas ich wykonania jest usunięty czas poświęcony na wykonywanie innych funkcji.

Aby pokazać wyniki optymalizacji kodu przeprowadzono następujące postępowanie. W dwóch projektach: zmiennoprzecinkowego filtru FIR w postaci bezpośredniej i zmiennoprzecinkowego filtru IIR zrealizowanego w postaci kaskadowej, dokonano kompilacji projektu z wyłączoną optymalizacji, a następnie dokonano ponownej kompilacji, tym razem używając optymalizacji na najwyższym poziomie jakości. Wyniki eksperymentu są przedstawione na poniższych wykresach.

Rys. 5.16. Porównanie czasu pracy filtru FIR przed i po optymalizacji.

Rys. 5.17. Porównanie czasu pracy filtru IIR przed i po optymalizacji.

W przypadku filtru FIR uzyskano przyśpieszenie działania o 96%, zaś w przypadku filtru IIR o 76%. Widać więc, że CCS w sposób bardzo efektywny dokonuje optymalizacji szybkości programu napisanego w języku C.

5.3. Zrealizowane struktury filtrów cyfrowych

5.3.1 Filtr FIR w postaci bezpośredniej

Postać bezpośrednia filtru FIR pochodzi bezpośrednio z równania opisującego filtr FIR i dokonuje wprost obliczania splotu sygnału z odpowiedzią impulsową filtru, która jest zapisana w postaci współczynników filtru. Postać tego filtru jest przedstawiona na rysunku 5.18.

Rys.5.18. Struktura bezpośrednia filtru FIR.

Filtr FIR zrealizowano w postaci stało i zmiennoprzecinkowej:


• postać zmiennoprzecinkowa
  
  sum=x*h[0]; //mnożenie wejścia przez pierwszy współczynnik filtru

for (i = nh-1; i != 0; i--){ //dokonanie splotu sygnału z resztą wsp. filtru

sum += h[i]*delay[i-1]; //mnożenie współczynnika z zapamiętaną próbką wejściową

delay[i-1]=delay[i-2]; //przesunięcie zapamiętanych próbek wejściowych

}

delay[0]=x; //wpisanie do bufora kołowego aktualnej próbki wejściowej

*y=sum; //przepisane wyniku splotu na wyjście


• postać stałoprzecinkowa
  
  sum=x*h[0]; //mnożenie wejścia przez pierwszy współczynnik filtru

for (i = nh-1; i != 0; i--){ //dokonanie splotu sygnału z resztą wsp. filtru

sum += h[i]*delay[i-1]; //mnożenie współczynnika z zapamiętaną próbką wejściową

delay[i-1]=delay[i-2]; //przesunięcie zapamiętanych próbek wejściowych

}

delay[0]=x; //wpisanie wartości wejściowej do bufora kołowego

*y=(sum<<9)>>24; //przepisanie wyniku na wejście wraz ze zmianą pozycji ułamka

Oba filtr bardzo dobrze oddają charakterystykę częstotliwościową filtru zaprojektowanego w MATLAB-ie. Ponadto oba filtry mają liniową fazę w swoich pasmach przepustowych, ponieważ ich odpowiedź impulsowa jest symetryczna. W filtrze stałoprzecinkowym sumowanie przeprowadzona na liczbie typu long, co odpowiada liczbie stałoprzecinkowej o długości 40 bitów. Zastosowanie tej liczby daje 9 bitów zapasu, co zapobiega powstawaniu błędów nadmiaru dla filtrów o niezbyt dużych rzędach. Współczynniki filtru stałoprzecinkowego są zawsze ułamkami w formacie 0.15, co znacznie ułatwia napisanie filtru.

5.3.2 Filtry FIR zrealizowane w postaci transponowanej

Struktura transponowana powstaje poprzez zmianę zwrotu krawędzi w grafie przepływowym opisującym filtr cyfrowy zrealizowany w postaci bezpośredniej. Wtedy po zamianie wejścia z wyjściem otrzymana struktura posiada identyczną transmitancję jak struktura w postaci bezpośredniej. Należy zwrócić uwagę, że węzły sumowania istniejące w strukturze bezpośredniej zamieniają się na węzły doprowadzające sygnał do elementów filtru i na odwrót, węzły doprowadzające sygnał stają się teraz węzłami sumowania. Struktura ma takie same wymagania zarówno pod względem liczby operacji, jak i obszaru pamięci potrzebnego do zrealizowania bufora kołowego. O ile filtr w postaci bezpośredniej pamięta próbki wejściowe filtru, to filtr w postaci transponowanej pamięta wyniki pośrednie obliczania splotu.

Rys.5.19 Filtr FIR w postaci transponowanej

Filtr FIR zrealizowano w postaci stało i zmiennoprzecinkowej:

• filtr zmiennoprzecinkowy


*y=h[0]*x+delay[0]; //policzenie próbki wejściowej
for (i = 0; i < nh-2; i++){ //liczenie kolejnych wyników pośrednich

delay[i]=h[i+1]*x+delay[i+1];} //obliczenie wartości bufora kołowego

delay[i]=h[i+1]*x; //ostatnia komórka bufora kołowego


• filtr stałoprzecinkowy


*y=((h[0]*x)>>15)+delay[0]; //policzenie próbki wejściowej

for (i = 0; i < nh-2; i++){ //liczenie kolejnych wyników pośrednich

delay[i]=((h[i+1]*x)>>15)+delay[i+1];} //obliczenie wartości bufora kołowego

delay[i]=(h[i+1]*x)>>15; //ostatnia komórka bufora kołowego

Charakterystyki częstotliwościowe obu filtrów są dobre. Prawdopodobieństwo powstania nadmiaru w filtrze stałoprzecinkowym jest niskie, ale dla dużego rzędu filtru mogą powstawać zniekształcenia spowodowane faktem, że pośrednie wartości przechowywane w komórkach bufora kołowego mogą znacząco różnić się do wartości idealnych

5.3.3 Struktury dla filtrów FIR o symetrycznej odpowiedzi impulsowej

Struktura ta bazuje na właściwości symetrii odpowiedzi impulsowej filtru FIR. Bazując na fakcie, że odpowiedź filtru spełnia warunek h[n] = h[M - 1 - n ], możemy zredukować liczbę operacji mnożenia w strukturze filtru o blisko połowę. Struktura pokazana jest na rysunku 5.20.

Rys.5.20. Strukura filtru FIR wykorzystująca symetrię odpowiedzi impulsowej.

Struktura filtru różni się dla parzystej i nieparzystej długości filtru. Filtr FIR zrealizowano w postaci stało i zmiennoprzecinkowej dla obu parzystości filtru:

• filtr zmiennoprzecinkowy dla parzystej długości filtru:


sum = (delay[2*nh-2]+we)*h[0]; // wejście + ostatnia zapamiętana * h[0]

temp=we; //zapamiętanie wejścia


for (i = 0; i < nh-1; i++){ //pętla licząca splot

sum += (delay[i]+delay[2*nh-3-i])*h[i+1];

}

for(i=0;i<2*nh-1;i++){ // bufor kołowy

temp1=delay[i];

delay[i]=temp;

temp=temp1;}

*y = sum; //przesłanie wyniku

• filtr zmiennoprzecinkowy dla nieparzystej długości filtru:
  
  sum = (delay[2*nh-2]+x)*h[0]; //pomnożenie wejścia i ostatniej komórki bufora przez h[0]

delay[2*nh-2]=delay[2*nh-3]; //zmiana zawartości wykorzystanej komórki bufora
temp=x; //zapamiętanie wejścia
for (i = 1; i < nh-1; i++){ //realizacja splotu zapamiętanych wartości wejścia i h[n]

sum += (delay[i-1]+delay[2*nh-2-i])*h[i]; //mnożenie sumy komórek bufora i odp. h[n]

delay[2*nh-2-i]=delay[2*nh-3-i]; //zmiana bufora po stronie większej od połowy rzędu filtru

temp1=delay[i-1]; //zmiana komórek bufora po stronie mniejszej od połowy

delay[i-1]=temp; // rzędu filtru

temp=temp1;

}

sum += delay[2*nh-2-i]*h[i]; //policzenie jedynej niesymetrycznej próbki

delay[2*nh-2-i]=temp; //wpisanie próbki wejściowej do bufora kołowego

*y = sum; //przesłanie wyniku


• filtr stałoprzecinkowy dla parzystej długości filtru:
  
  sum = (delay[2*nh-2]+we)*h[0]; // wejście +ostatnia zapamiętana * h[0]

temp=we; //zapamiętanie wejścia

for (i = 0; i < nh-1; i++){ //pętla licząca splot

sum += (delay[i]+delay[2*nh-3-i])*h[i+1];

}

for(i=0;i<2*nh-1;i++){ // bufor kołowy

temp1=delay[i];

delay[i]=temp;

temp=temp1;}

*y = (sum<<9)>>24; //przesłanie wyniku ze zmianą pozycji ułamka


• filtr stałoprzecinkowy dla nieparzystej długości filtru:


sum = (delay[2*nh-2]+x)*h[0]; //pomnożenie wejścia i ostatniej komórki bufora przez h[0]

delay[2*nh-2]=delay[2*nh-3]; //zmiana zawartości wykorzystanej komórki bufora
temp=x; //zapamiętanie wejścia
for (i = 1; i < nh-1; i++){ //realizacja splotu zapamiętanych wartości wejścia i h[n]

sum += (delay[i-1]+delay[2*nh-2-i])*h[i]; //mnożenie sumy komórek bufora i odp. h[n]

delay[2*nh-2-i]=delay[2*nh-3-i]; //zmiana bufora po stronie większej od połowy rzędu filtru

temp1=delay[i-1]; //zmiana komórek bufora po stronie mniejszej od połowy

delay[i-1]=temp; // rzędu filtru

temp=temp1;

}

sum += delay[2*nh-2-i]*h[i]; //policzenie jedynej niesymetrycznej próbki

delay[2*nh-2-i]=temp; //wpisanie próbki wejściowej do bufora kołowego

*y = (sum<<9)>>24; //przesłanie wyniku ze zmianą pozycji ułamka

Filtry zarówno w postaci stałoprzecinkowej, jak i w postaci zmiennoprzecinkowej dobrze odzwierciedlają zadaną charakterystykę amplitudową. Współczynniki filtrów stałoprzecinkowych są podawane w formacje 0.15. Niestety znaczna komplikacja pętli obliczające splot powoduje, że te filtry po optymalizacji pracują wolniej od filtrów w postaci bezpośredniej, czy transponowanej.

5.3.4. Struktury flitów FIR wykorzystujące antysymetrię odpowiedzi impulsowej

Podobnie, jak w strukturze omówionej wcześniej wykorzystujemy fakt, że odpowiedź impulsowa filtrów FIR typu III i IV posiada cechę antysymetrii, to jest h[n] = - h[M - 1 - n ]. Różnica w strukturze filtru polega na zamianie dodawania pamiętanych próbek wejściowych na ich odejmowanie. Strukturę filtru przedstawiona na rysunku 2.21.

Rys. 5.21. Struktura filtru wykorzystująca asymetrię odpowiedzi impulsowej.

Ponieważ dla nieparzystego rzędu filtru, nieparzysty współczynnik jest zawsze zerem struktura filtru nie różni się dla filtru parzystego i nieparzystego. Poniżej podano realizację filtru dla liczb zmiennoprzecinkowych i stałoprzecinkowych:

• filtr zmiennoprzecinkowy

sum = (we-delay[2*nh-2])*h[0]; //wejście +ostatnia zapamiętana * h[0]
temp=we; //zapamiętanie wejścia

for (i = 0; i < nh-1; i++){ //policzenie splotu

sum += (delay[i]-delay[2*nh-3-i])*h[i+1];}
for(i=0;i<2*nh-1;i++){ //realizacja bufora kołowego

temp1=delay[i];

delay[i]=temp;

temp=temp1;}
*y = sum; //przypisanie wyjścia

• filtr stałoprzecinkowy
  sum = (we-delay[2*nh-2])*h[0]; // wejście +ostatnia zapamiętana * h[0]
  temp=we; //zapamiętanie wejścia
  for (i = 0; i < nh-1; i++){ //policzenie splotu
  sum += (delay[i]-delay[2*nh-3-i])*h[i+1];
  }
  for(i=0;i<2*nh-1;i++){ //realizacja bufora kołowego
  temp1=delay[i];
  delay[i]=temp;

temp=temp1;}
*y = (sum<<9)>>24; //przesłanie wyniku ze zmianą pozycji ułamka

Filtry dobrze oddają charakterystykę częstotliwościową potrzymana w MATLAB-ie.

5.3.5 Kaskadowa realizacja filtrów FIR

Filtry FIR można przedstawić w postaci kaskady filtrów FIR drugiego rzędu. Poprzez sparowanie zer sprzężonych współczynniki filtru są rzeczywiste, ale ich wartość może przyjmować wartości większe od jedności. Jest to poważny problem przy realizacji filtru stałoprzecinkowego, ponieważ wymusza zastosowanie liczb o formacie różnym od 0.15 lub zastosowanie skalowania filtru. Skalowanie współczynników filtru wymusza stosowanie dodatkowych mnożników pomiędzy stopniami, aby zachować amplitudę sygnału. Rozwiązanie kaskadowe jest stosowane najczęściej dla filtrów o rzędach powyżej 256, w celu zmniejszenia wrażliwości filtru na skutki kwantyzacji współczynników oraz zmniejszenie prawdopodobieństwa nadmiaru. Filtry zrealizowano w postaci stało- i zmiennoprzecinkowej:

• filtr zmiennoprzecinkowy

tempwe=we; //przepisanie wejścia

for (i = 0; i < nh1; i++) {

tempwe=tempwe*skala[i]; //skalowanie wejścia do kolejnej sekcji

tempwy = tempwe*b[i*3+0]; //pierwszy współczynnik sekcji

tempwy += b[i*3+1]*del1[i*2]; //drugi współczynnik sekcji

tempwy += b[i*3+2]*del1[i*2+1]; //trzeci współczynnik sekcji

del1[i*2+1]=del1[i*2]; //bufor kołowy

del1[i*2]=tempwe;

tempwe=tempwy; //przepisanie wyjścia sekcji do wejścia następnej

}

*wy=tempwe; //wyjście z filtru

• filtr stałoprzecinkowy
  tempwe=(we*skalawe)>>15; //skalowanie wejścia;

for (i = 0; i < nh1; i++) { //kolejne sekcje

tempwy = tempwe*b[i*3+0]; //pierwszy współczynnik sekcji

tempwy += b[i*3+1]*del1[i*2]; //drugi współczynnik sekcji

tempwy += b[i*3+2]*del1[i*2+1]; //trzeci współczynnik sekcji

del1[i*2+1]= del1[i*2]; //bufor kołowy

del1[i*2]=tempwe;

tempwe=(tempwy<<9)>>24; //przesuwanie pozycji kropki

tempwe=tempwe<<przes[i]; //zastosowanie skalowania

}

*wy=(tempwy<<9)>>24; //przesłanie wyniku na wyjście z przesunięciem pozycji //kropki

O ile filtr kaskadowy w postaci zmiennoprzecinkowej sprawuje się dobrze, to filtr stałoprzecinkowy przestaje dobrze działać dla filtru o rzędzie większym od 20. Kiedy usuniemy współczynniki skalujące filtr nie filtruje w ogóle , na wyjściu są same zera, natomiast jeśli je zostawimy to na wyjściu pojawia się przebieg przypadkowy. Na wykresie odpowiedzi impulsowej filtru widać, że filtr pracuje początkowo prawidłowo, ale przy wzroście sygnału na wyjściu pojawiają się duże oscylacje. Z takiego zachowania się filtru można wyciągnąć wniosek, że powodem złej pracy filtru jest pojawianie się nadmiaru.

5.3.6. Filtry FIR zrealizowane w postaci kratowej

Postać kratowa filtru cyfrowego wywodzi się z filtracji adaptacyjnej. Filtr FIR można przedstawić w postaci kaskady filtrów kratowych, których sekcja numer l jest opisana następującymi równaniami:

Fl(n) = x(n) + K1 * x(n-1) Gl(n) = Kl * x(n) + x(n-1)

(4.1)

Współczynniki filtru kratowego uzyskuje w trakcie procedury iteracyjnej. Filtry kratowe maja jedno poważne ograniczenie nie mogą posiadać zer na kręgu jednostkowym. Ponadto dla zer położonych poza okręgiem jednostkowym współczynnik K osiąga duże wartości. Dlatego stuk tura kratowa jest strukturą odpowiednia dla filtrów minimalnofazowych, czyli posiadających wszystkie zera wewnątrz okręgu jednostkowego. W tej pracy napisano funkcje, która w prosty sposób przekształca filtr FIR w postać minimalno fazową kosztem dwóch parametrów : liniowości fazy oraz zniekształcenia charakterystyki w paśmie zaporowym. Struktura potrzebuje dwa razy więcej mnożników niż struktura bezpośrednia. Wymagania dotyczące pamięci są takie same.

Rysunek 5.22 przedstawia filtr FIR zrealizowany w postaci kratowej:

Rys.5.22. Struktura kratowa filtr FIR.

Poniżej podano realizację filtru dla liczb zmiennoprzecinkowych i stałoprzecinkowych:

• filtr zmiennoprzecinkowy

f=x; //przepisanie wejścia do wejść filtru

g=x;

for (i = 0; i < nh; i++){ //kolejne sekcje

temp=g; //zachowanie wartości g

g = h[i]*f+delay[i]; //nowa wartość g

f = f+h[i]*delay[i]; //nowa wartość f

delay[i]=temp; //zapamiętanie stare wartości g

}

*wy=f; //wartość wyjścia

• filtr stałoprzecinkowy

f=x<<15; //przekształcenie wejścia z short na long

g=x<<15;

for (i = 0; i < nh; i++){ //kolejne sekcje

temp=g; //zapamiętaj czasowo g

g=delay[i]<<15; //podstaw za g wartość opóźnienia

g += h[i]*((f<<9)>>24); //nowa wartość g

f += h[i]*delay[i]; //nowa wartość f

delay[i]=(temp<<9)>>24; //zapamiętanie wartości g

}

*wy=(f<<9)>>24; //przypisanie wyjścia

Przy realizacji stałoprzecinkowej należy zwrócić uwagę na możliwość powstania nadmiaru w linii obliczającej wartość f. W tej linii dokonuje się tyle operacji dodawania, ile wynosi rząd filtru. W przedstawionej powyżej wartość jest przechowywana w słowie o długości 40 bitów, co daje 9 bitów zapasu.

5.3.7 Pierwsza struktura bezpośrednia filtry IIR

Pierwszą strukturę bezpośrednia wyprowadza się wprost ze transmitancji filtru IIR. Współczynniki tej struktury są wprost współczynnikami przy kolejnych potęgach z w wielomianach transmitancji. Wygląd struktury jest przestawiona na rysunku 5.23.

Rys. 5.23. Struktura bezpośrednia filtru IIR.

Poniżej przedstawiono filtry IIR w postaci bezpośredniej zmiennoprzecinkowe i stałoprzecinkowe:

• filtr zmiennoprzecinkowy

w = x*b[0]; //w=b0*probka wejściowa

for (i=1;i<Nb;i++) //pętla obliczającą cześć MA

w +=b[i]*db[i-1];

for (i=1;i<Na;i++) //obliczanie części AR

w += a[i]*da[i-1];

for (i=Nb-1;i>0;i--) //bufor kołowy cześć MA

db[i]=db[i-1];

for (i=Na-1;i>0;i--) //bufor kołowy cześć AR

da[i]=da[i-1];

db[0]=x; //wejście na początek bufora części MA

da[0]=w; // wejście na początek bufora części MA

*wy=w; //przepisanie rezultatu

• filtr stałoprzecinkowy

w = (x>>SHIFT1)*b[0]; //w=b0*probka wejściowa

for (i=1;i<Nb;i++) //pętla obliczająca cześć MA

w +=b[i]*db[i-1];

for (i=1;i<Na;i++) //obliczanie części AR

w -= a[i]*da[i-1];

y=(w<<SHIFT2)>>24; //wyjście w formacie używanym w filtrze

for (i=Nb-1;i>0;i--) //bufor kołowy cześć MA

db[i]=db[i-1];

for (i=Na-1;i>0;i--) //bufor kołowy cześć AR

da[i]=da[i-1];

db[0]=x>>SHIFT1; //wejście w odpowiednim formacie wpisywane do bufora

da[0]=y; //wyjście wpisywane do bufora

*t=(w<<SHIFT3)>>24; //wyjście w formacie 0.15

Filtr stałoprzecinkowy używa liczb w formacie X.15-X, dlatego wymaga odpowiedniego przesuwania pozycji ułamka dziesiętnego przy wykonywaniu operacji arytmetycznych. Liczba X rośnie wraz ze wzrostem rzędu filtru, powodując stopniową degradację charakterystyki częstotliwościowej, poprzez zmianę położenia zer i biegunów. Jest to szczególnie widoczne w przypadku zer, które zazwyczaj są ułamkami, przez co wraz ze wzrostem liczby X służącej do reprezentacji liczb większych od jedności, ich wartość małe od zera. Ta sama struktura w postaci zmiennoprzecinkowej zachowuje kształt charakterystyki częstotliwościowej..

5.3.8. Druga struktura bezpośrednia filtru IIR

Drugą strukturę bezpośrednią wyprowadza się poprzez zamianę miejscami w pierwszej strukturze bezpośredniej części zawierającej filtr o samych zerach i filtr o samych biegunach. Można wtedy zauważyć, że strukturę można uprościć poprzez połączenie dwóch buforów kołowych w jeden. W ten sposób osiągamy oszczędność pamięci. Struktura bezpośrednia druga jest przedstawiona na rysunku 5.24.

Rys 5.24. Druga struktura bezpośrednia filtru typu IIR.

Tę strukturę także zrealizowano w postaci stało- i zmiennoprzecinkowej:

• filtr zmiennoprzecinkowy

sumwe = x; //wpisanie wejścia do zmiennej wykonującej pierwsze sumowanie

for(i=0;i<na;i++){ //liczenie części AR

sumwe -=a[i]*da[i];

}

sum=sumwe*b[0]; //wpisanie wyniku pośredniego do zmiennej drugiego sumowania

for(i=0;i<nb;i++){ //liczenie części MA

sum +=b[i]*da[i-1];}

*wy=sum; //przypisanie wyjścia

for (i=Na-1;i>0;i--) //wspólny bufor kołowy

da[i]=da[i-1];

da[0]=sumwe; //wynik pierwszego sumowania wpisany do bufora kołowego

• filtr stałoprzecinkowy

pos = (x<<SHIFT1); //wpisanie wejścia do zmiennej wykonującej pierwsze sumowanie

for(i=1;i<na;i++){ //obliczenie części AR

pos+=a[i]*da[i-1];}

pos1=(pos<<SHIFT2)>>24;// zmienna pośredniczącej pomiędzy obiema częściami filtru

pos=pos1*b[0]; //wpisanie wyniku pośredniego do zmiennej dokonującej drugie //sumowanie

for(i=0;i<nb;i++) //liczenie części MA

pos+=b[i]*da[i-1];

*t=(pos<<SHIFT3)>>24; //policzenie wyjścia

for (i=Na-1;i>0;i--) //wspólny bufor kołowy

da[i]=da[i-1];

da[0]=pos1; //wartość pośrednia na początek bufora kołowego

Podobnie jak dla poprzedniej struktury w filtrze stałoprzecinkowym używamy liczb postaci X.15-X. Charakterystyka tego filtru wraz ze wzrostem rzędu jest coraz bardziej zniekształcona, powodu zmniejszania się liczby bitów użytych do reprezentacji współczynników filtru, co z kolei powodu przesuwanie zer i biegunów filtru. Postać zmiennoprzecinkowa tego filtru dobrze oddaje kształt charakterystyki częstotliwościowej.

5.3.9. Pierwsza struktura bezpośrednia transponowana filtru IIR

Podobnie jak dla filtru FIR strukturę transponowaną dla filtru IIR uzyskuje się poprzez dokonanie inwersji zwrotu krawędzi w grafie przepływu ułożonym dla tej struktury oraz zamianę miejscami wejścia i wyjścia. Otrzymaną strukturę przedstawia rysunek 2.25.

Rys 2.25 Pierwsza struktura bezpośrednia transponowana filtru typu IIR.

Poniżej przedstawiono ciała funkcji realizujących tę strukturę filtru:

• filtr zmiennoprzecinkowy

pos1=x+da[0]; //suma wejścia i ostatniego opóźnienia bufora

*wy=pos1*b[0]+db[0]; //policzone wyjście

for(i=1;i<nb-1;i++) //kolejne wartości opóźnień części MA

db[i-1]=pos1*b[i]+db[i];

db[nb-2]=pos1*b[nb-1]; //pierwsze opóźnienie części MA

for(i=1;i<na-1;i++){ //kolejne wartości opóźnień części AR

da[i-1]=pos1*-a[i]+da[i];

da[na-2]=pos1*-a[na-1]; //pierwsze opóźnienie części AR

• filtr stałoprzecinkowy

tempwe = x<<SHIFT1; //przekształcenie wejścia do formatu używanego w filtrze

pos=tempwe+da[0]; //suma wejścia i ostatniego opóźnienia bufora części AR

pos1=(pos<<SHIFT2)>>24; //zmiana długości słowa by umożliwić mnożenie

pos=pos1*b[0]+db[0]; //policzenie wartości wyjściowej

*t=(pos<<SHIFT2)>>24; //przekształcenie wyniku do formatu używanego na wyjściu //filtru

for(i=1;i<nb-1;i++){ //policzenie wartości opóźnień części MA

pos=pos1*b[i]+db[i]; //węzeł sumowania

db[i-1]=(pos<<8)>>8;} //przekształcenie na int

db[nb-2]=pos1*b[nb-1]; //ostatnia komórka bufora

for(i=1;i<na-1;i++){ //policzenie wartości opóźnień części AR

pos=pos1*a[i]+da[i]; //węzeł sumowania

da[i-1]=(pos<<8)>>8;} //przekształcenie na int

da[na-2]=pos1*a[na-1]; //ostatnia komórka bufora

Filtr stałoprzecinkowy ma bardzo złe właściwości. Można na nim realizować filtry, co najwyżej czwartego rzędu. Powodem tego zachowania się filtru są duże zniekształcenia w obliczaniu wartości pośrednich. Filtr zmiennoprzecinkowy zachowuje się poprawnie.

5.3.10. Druga struktura bezpośrednia transponowana filtru IIR

Także tę strukturę uzyskujemy poprzez dokonanie przekształcenia grafu przepływu sygnału drugiej struktury bezpośredniej filtru IIR. Otrzymana struktura jest przedstawiona na rysunku 5.26.

Rys.5.26. Druga struktura transponowana filtru IIR.

Filtry zrealizowano w postaci stało- i zmiennoprzecinkowej:

• filtr zmiennoprzecinkowy

pos=b[0]*x+da[0]; //policzenie wyjścia

*wy=pos; //przypisanie wartości wyjściowej

for(i=1;i<na-1;i++){ //liczenie kolejnych opóźnień

da[i-1]=b[i]*x-a[i]*pos+da[i];

}

da[i-1]=b[i]*x-a[i]*pos; //wartość ostatniego opóźnienia

• filtr stałoprzecinkowy
  

• pos1 = x>>SHIFT1; //ustawienie odpowiedniego formatu słowa wejściowego

pos=b[0]*pos1+(da[0]<<15); //policzenie wartości wyjściowej

y=(pos<<SHIFT2)>>24; //ustawienie formatu zmiennej podanej na mnożniki pętli //AR

*t=y<<SHIFT1; // podanie wyniku na wyjście

for(i=1;i<na-1;i++){ //liczenie zawartości bufora kołowego

pos=b[i]*pos1+a[i]*y+(da[i]<<15);

da[i-1]=(pos<<SHIFT2)>>24; //zamiana na shorta

}

pos=b[i]*pos1+a[i]*y; //ostatnia komórka bufora opóźnień

da[i-1]=(pos<<SHIFT2)>>24;

Ten filtr w postaci zmiennoprzecinkowej zachowuje się poprawnie, natomiast postać stałoprzecinkowa zachowuje się fatalnie, podobnie jak poprzednia struktura transponowana.

5.3.11. Struktura kaskadowa filtru IIR

Każdą opisaną powyżej strukturę bezpośrednią można rozłożyć na kaskadę złożona z filtrów drugiego rzędu, w taki sposób, aby jego transmitancja H(z) mogła zostać wyrażona następującym wzorem:


(4.2)

Współczynniki sekcji drugiego rzędu są rzeczywiste, ponieważ każdą sekcję tworzy albo para biegunów i zer sprzężonych lub jest to para biegunów i zer rzeczywistych. Struktura kaskadowa ma małą wrażliwość na skutki kwantyzacji współczynników filtru, ponieważ współczynniki {b_k} i {a_k} wprost wyznaczają położenie zer i biegunów transmitancji. Ich zaletą jest możliwość łatwego skalowania struktury, aby uniknąć występowania nadmiaru. Przykład typowej struktury kaskadowej pokazuje rysunek 5.27.

Rys. 5.27. Stuktura kaskadowa filtru typu IIR.

Poniżej przedstawiono realizacje struktury kaskadowej stało- i zmiennoprzecinkowej:

• filtr zmiennoprzecinkowy

tempwe=(*r)*skala[0]; //skalowanie wejścia

for (i = 0; i < n; i++) //kolejne sekcje

{

tempwy=((c[6*i+0]*tempwe)+d[2*i+0])*c[6*i+3]; //obliczenie wyjścia

d[2*i+0]=c[6*i+1]*tempwe+tempwy*c[6*i+4]+d[2*i+1];//obliczenie nowych opóźnień

d[2*i+1]=c[6*i+2]*tempwe+tempwy*c[6*i+5];

tempwe=tempwy*skala[i+1]; //wejście do nowe sekcji

}

*r=tempwy; //przypisanie wyjścia

• filtr stałoprzecinkowy


tempwe=*r; //prowadzenie próbki wejściowej do pętli

for (i = 0; i < n; i++) {

tempsum=((c[5*i+0]*tempwe)+d[2*i+0]);//obliczenie wyjścia

tempsum=tempsum<<1; //pomnożenie go przez dwa

tempwy=(tempsum<<9)>>24; //przesunięcie pozycji ułamka wyniku

d[2*i+0]=c[5*i+1]*tempwe+tempwy*c[5*i+3]+d[2*i+1];//1 opóźnienie

d[2*i+1]=c[5*i+2]*tempwe+tempwy*c[5*i+4]; //drugie opóźnienie

tempwe=tempwy; //połączenie między kaskadami

}

*r=tempwy; //wyjście z filtru

Filtr zmiennoprzecinkowy jest w pełni skalowany, natomiast filtr stałoprzecinkowy jest przeskalowany przez czynnik równy dwu, aby wszystkie współczynniki filtru były liczbami w formacie 0.15. Sekcje filtru zostały zrealizowane w postaci filtrów o drugiej strukturze bezpośredniej transponowanej. Oba filtry dobrze oddają charakterystyki częstotliwościowe filtru.

5.3.12 Struktura równoległa filtru IIR

Struktura równoległa powstaje na skutek rozłożenia transmitancji filtru na ułamki proste. W wyniku otrzymujemy równoległe połączenie sekcji drugiego rzędu oraz stałej. Każda sekcja drugiego rzędu zawiera co najwyżej cztery współczynniki. Struktura równoległa nie nadaje się do filtru, dla których najważniejszym parametrem jest tłumienie. Powodem jest fakt, że dla tej struktury na położenie pojedynczego zera mają wpływ współczynniki {b_k} wszystkich sekcji połączonych równoległe.

Rys.5.28. Struktura równoległa filtru IIR.

Poniżej przedstawiono realizację struktury równoległej w postaci kodu napisanego w języku C:

• filtr zmiennoprzecinkowy

tempwe=(*r); //przypisanie wejścia

suma=tempwe*wsp; //wolny współczynnik

for (i = 0; i < n; i++) // kolejne sekcje

{

tempsum=c[4*i+2]*d[2*i+0]+c[4*i+3]*d[2*i+1]+tempwe;//policzenie wartości pośredniej

pos1=tempsum;

pos1=pos1*2; //przemnożenie przez dwa ponieważ a₀=0.5

suma += c[4*i+0]*pos1 + c[4*i+1]*d[2*i+0];//dodanie wyniku do wyjścia

d[2*i+1]=d[2*i]; //bufor kołowy sekcji

d[2*i]=pos1;

}

*r=suma; //przypisanie wyniku

• filtr stałoprzecinkowy

tempwe=(*r); //pobranie wartości wejścia

suma=tempwe*wsp; //wolny współczynnik

for (i = 0; i < n; i++) { //kolejne sekcje

tempsum=tempwe<<15; //zmiana pozycji przecinka w celu poprawnego dodawania

tempsum+=c[4*i+2]*d[2*i+0]+c[4*i+3]*d[2*i+1];//Sumowanie współczynników a

pos1=(tempsum<<9)>>24; //wynik obcięty do 16 bitów

pos1=pos1<<1; //pomnożenie wyniku przez 2

suma += c[4*i+0]*pos1 + c[4*i+1]*d[2*i+0]; //sumowanie współczynników B

d[2*i+1]=d[2*i]; //bufor kołowy sekcji

d[2*i]=pos1;}

*r=(suma<<9)>>24; //przypisanie wyniku na wyjście ze zmianą pozycji ułamka //dziesiętnego

Struktura równoległa została zrealizowana przy pomocy sekcji w postaci drugiej struktury bezpośredniej drugiego rzędu. Aby uzyskać liczby w formacie 0.15 oba filtry zostały przeskalowane przez współczynnik 0.5. Filtry dość dobrze oddają charakterystykę częstotliwościową filtru.

5.3.13. Struktura kratowa filtru IIR

Strukturę kratową dla filtru IIR tworzy się poprzez odpowiednie przekształcenie drugiej struktury bezpośredniej. Cześć odpowiadająca mianownikowi transmitancji jest przekształcana na filtr kratowy.. Po dokonaniu odpowiedniej zmiany we współczynnikach licznika transmitancji otrzymujemy filtr przedstawiony na rysunku 5.29. Struktura kratowa ma minimalne wymagania dotyczące pamięci, ale nie jest optymalna, jeżeli chodzi o liczbę mnożeń potrzebnych do zaimplementowania struktury. Zaletami struktur kratowych jest ich mała wrażliwość na skutki skończonej długości słowa oraz wbudowana stabilność gwarantowana sposobem obliczania współczynników lustrzanych.

Rys 5.28. Struktura kratowa dla filtrów IIR

Filtr IIR w postaci kratowej także zrealizowano w dwóch formatach reprezentacji liczb binarnych:

• filtr zmiennoprzecinkowy

f=r; //wprowadzenie wartości wejścia

sum=0; //zerowanie sumy wsp. ladder

for (i = 0; i < n; i++){ //kolejne sekcje

f -= k[i] * delay[i+1]; //nowa wartość f

g = k[i]*f + delay[i+1]; //nowa wartość g

delay[i]=g; //zapamiętanie obliczonej wartości

sum +=v[i]*delay[i];} //policzenie sumy

delay[n]=f; //zamknięcie pętli AR

sum +=v[n]*delay[n]; //ostatni współczynnik ladder

return sum; //zwróć wynik

• filtr stałoprzecinkowy

f=r<<15; //wprowadzenie wartości wejścia z zamianą na longa

sum=0; //zerowanie sumy wsp. ladder

for (i = 0; i < n; i++){ //kolejne sekcje

f -= k[i] * delay[i+1]; //nowa wartość f

g = k[i]*((f<<9)>>24) + (delay[i+1]<<15); //nowa wartość g

delay[i]=(g<<9)>>24; // zapamiętanie obliczonej wartości z zamianą na //short-a

sum +=v[i]*delay[i]; //policzenie sumy

}

delay[n]=(f<<9)>>24; //zamknięcie pętli AR

sum +=v[n]*delay[n]; //ostatni współczynnik ladder

return (sum<<9)>>24; //zwróć wynik

Struktura kratowa w postaci zmiennoprzecinkowej zachowuje się dobrze, natomiast w postaci stałoprzecinkowej często ma zniekształcone pasmo zaporowe. Współczynniki struktury kratowej dobrze nadają się do realizacji stałoprzecinkowej, bo są zawsze mniejsze od jedności, co znacznie ułatwia używanie arytmetyki stałoprzecinkowej.

5.3.14 Struktura zmiennych stanu filtru IIR

Dotychczas rozważane struktury opierały się na opisie filtru z punktu widzenia zależności pomiędzy wyjściem, a wejściem. Istnieje możliwość opisu układu ciągle zawierającej zależność pomiędzy wejściem i wyjściem, ale wprowadzających dodatkowy zbiór zmiennych, zwanych zmiennymi stanu. Co więcej, równania opisujące układ są podzielone na dwie części:

1. zbiór równań opisujących zależność zmiennych stanu od aktualne wartości wejściowej,

2. zbiór równań przestawiających zależność wyjścia od zmiennych stanu i wartości wejściowych.

Zmienne stanu dostarczają informacji o sygnałach we wszystkich węzłach wewnętrznych układu. W wyniku struktura zmiennych stanu dostarcza bardziej szczegółowego opisu systemu niż wszystkie pozostałe. Prowadzi to otrzymania filtru o minimalnej liczbie zarówno operacji mnożenia, jak i liczbie komórek pamięci potrzebnej do dokonania implementacji filtru. Ponadto szum kwantyzacji dla tej struktury jest najmniejszy.

Rys.5.29. Struktura zmiennych stanu filtru IIR.

Poniżej są przedstawione programy realizujące omówioną strukturę:

• filtr zmiennoprzecinkowy

temp=0;

temp1=*r; //pamiętaj wejście

for (i = 0; i < n; i++) //policz wyjście

temp += C[i]*v[i];

*r=temp+(*r)*D;

temp=v[0]; //policz nowe zmienne stanu

v[0]=v[0]*A[0]+temp1;

for (i = 1; i < n; i++)

v[0 ]+= v[i]*A[i];

for (i = n-1; i > 0; i--)

v[i]=v[i-1];

v[1]=temp;

• filtr stałoprzecinkowy

tempwe=*r; //zapamiętaj wejście

for (i = 0; i < n; i++) //kolejne sekcje

{

w=(tempwe>>1)*d[i]; //policz wyjście sekcji

w+=c[2*i+0]*v[2*i+0]+c[2*i+1]*v[2*i+1];

tempwy=(w<<11)>>24;

w=tempwe<<13; //policz zmienne stanu sekcji

w+=a[2*i+0]*v[2*i+0]+a[2*i+1]*v[2*i+1];

v[2*i+1]=v[2*i+0];

v[2*i+0]=(w<<10)>>24;

tempwe=tempwy; //połączenie pomiędzy sekcjami

}

*r=tempwy;

Strukturę zmiennoprzecinkową jest tworzona wprost, natomiast struktura stałoprzecinkowa jest przedstawiona w postaci kaskady filtrów drugiego rzędu. Powodem takiego postępowania jest fakt, że zmienne stanu przyjmują duże wartości dla realizacji filtru wprost. Natomiast dla realizacji filtru w postaci kaskady dają się przedstawić w formacie 1.14. Oba filtry dobrze oddają zadaną charakterystykę częstotliwościową.




5.3.15 Porównanie szybkości działania zaimplementowanych filtrów

Rys.5.30. Porównanie szybkości działania filtrów zmiennoprzecinkowych FIR.




Rys.5.31. Porównanie szybkości działania filtrów stałoprzecinkowych FIR.

Filtry stałoprzecinkowe zazwyczaj są szybsze niż ich zmiennoprzecinkowe odpowiedniki. Najszybszym filtrem stałoprzecinkowym jest filtr w postaci bezpośredniej, zaś filtrem zmiennoprzecinkowym filtr zrealizowany w postaci zmiennoprzecinkowej. Całkowicie zawiodły filtry wykorzystujące symetrię filtru. Wynika stąd wniosek, że przy pisaniu kodu w CCS bardzo ważny jest jego styl.

Rys.5.32. Porównanie szybkości działania filtrów IIR zmiennoprzecinkowych.

Rys.5.33.Porównanie szybkości stałoprzecinkowych filtrów IIR.

Najszybszym filtrem zmiennoprzecinkowym jest filtr kaskadowy, zaś w postaci stałoprzecinkowej zbliżone rezultat został uzyskany dla filtrów w postaci kaskadowej i zmiennych stanu.




6. Podsumowanie

W pracy zostały przedstawione najbardziej popularne metody projektowania filtrów cyfrowych. Opisano ich sposób podejścia do problemu oraz podano szkielet algorytmu. Dokonano porównania ich właściwości. Znając zalety i ograniczenia poszczególnych sposobów uzyskiwania współczynników filtru projektantowi na wybranie metody najbardziej odpowiadającej jego potrzebom.

Zrealizowano, przy wykorzystaniu pakietu MATLAB 6.0, kilkanaście interfejsów GUI implementujących metody projektowania filtrów cyfrowych. Pozwalają one na otrzymanie współczynników filtru, spełniających wymagania narzucone przez użytkownika. Napisane zostały programy pozwalające zaprojektować

1. filtry o skończonej odpowiedzi impulsowej :

1. metodą okien,

2. metodą próbkowania charakterystyki,

3. metodą minimalizacji błędu średniokwadratowego,

4. metodą minimalizacji błędu maksymalnego (Remeza),

5. metodą wymuszonej minimalizacji błędu średniokwadratowego,

6. metodą Butterworth'a.

2. filtry o nieskończonej odpowiedzi impulsowej :

1. metodą prototypu analogowego,

2. metodą Yule-Walkera,

3. metodą p-tej normy,

4. metodą modelowania.

Stworzono możliwości oceny spełniania przez filtr wymagań na kształt charakterystyki częstotliwościowej, poprzez nałożenie na otrzymaną charakterystykę ram projektowych. Napisano także program umożliwiający implementację filtru w wybranej przez użytkownika strukturze. Użytkownik może na bieżąco określić wpływ skutków kwantyzacji na charakterystykę filtru oraz dobrać strukturę pozwalającą na zminimalizowanie skutków skończonej długości słowa. Ponadto istnieje możliwość sprawdzenia, czy zaprojektowany filtr jest stabilny, dokonać oceny wpływu szumu kwantyzacji na charakterystyki filtru oraz sprawdzić, czy dany filtr nie generuje cykli granicznych. Dla każdej aplikacji w napisanej w MATLAB-ie napisano plik pomocy, objaśniający przeznaczenie poszczególnych elementów okna implementującego daną metodę projektowania filtru cyfrowego

Ponadto zrealizowano szereg praktycznych realizacji struktur filtrów cyfrowych przy wykorzystaniu programu CCS. Używając języka C zostały napisane programy implementujące następujące struktury:

1. dla filtrów FIR:

1. strukturę bezpośrednią,

2. strukturę transponowaną,

3. strukturę kaskadową

4. strukturę wykorzystującą symetrię współczynników filtru dla rzędu parzystego i nieparzystego,

5. strukturę wykorzystującą antysymetrię współczynników filtru,

6. strukturę kratową,

2. dla filtrów IIR:

1. pierwszą strukturę bezpośrednią,

2. drugą strukturę bezpośrednią,

3. pierwszą strukturę bezpośrednią transponowaną,

4. drugą strukturę bezpośrednią transponowaną,

5. strukturę kaskadową,

6. strukturę równoległą,

7. strukturę kratową,

8. strukturę zmiennych stanu.

Każda ze struktur została zaimplementowana dla formatu zapisu liczb binarnych stałoprzecinkowego i zmiennoprzecinkowego. Dokonano oceny jakości otrzymanych filtrów cyfrowych poprzez sprawdzenie szybkości działania filtru oraz zgodności kształtu otrzymanej charakterystyki częstotliwościowej. Filtry cyfrowe zostały wykonane na procesorze sygnałowym TMS320C6711. Komunikacja pomiędzy aplikacjami stworzonymi w MATLAB-ie, a zaimplementowanymi filtrami odbywa się poprzez pliki nagłówkowe.

Stworzony pakiet programów obejmuje cały ciąg procesu projektowania filtrów cyfrowych: od wyboru metody projektującej idealny filtr cyfrowy do oceny jego praktycznej realizacji. Pozwala on użytkownikowi na zapoznanie się z najważniejszymi problemami projektowania filtrów cyfrowych. Dlatego z pewnością będzie przydatny w czasie pracy w laboratorium z przedmiotu „Cyfrowe przetwarzanie sygnałów”.




Bibliografia

1. Alan V. Oppenheim, Roland W. Schafer. „Cyfrowe przetwarzanie sygnałów”. WKŁ, Warszawa 1979.

2. John G. Prolakis, Dimitris G. Manolakis. „Digital Signal Processing. Principles, Algorithms and Applications”. Prentice Hall, Upper Sadle River, 1996

3. Leland B. Jackson. „Digital Filtrering and Signal Procsing”. Kluwer Academic Publishers, Norwell, 1986.

4. Dietrich Schlichtharle. „Digital Filter Basic and Design”. Springer, Berlin, 2000.

5. Friedlander, B., and B. Porat, “The Modified Yule-Walker Method of ARMA Spectral Estimation,” IEEE Transactions on Aerospace Electronic Systems, AES-20, No. 2 (Marzec 1984), pp. 158-173.

6. Zespół „Computer-Based Exerscises for Signal Processing Using MATLAB”. Printice-Hall, Englewood Cliffs, 1994.

7. Richard G. Lyons. „Wprowadzenie do cyfrowego przetwarzania sygnału”. WKŁ, Warszawa, 2000.

8. A.G. Deczky. „Synthesis of digital recursive filters using the minimum P error criterion”. IEEE Trans. on Audio and Electroacoustics, AU-20(2):257, October 1972

9. Steiglitz, K., and L.E. McBride, “A Technique for the Identification of Linear Systems,” IEEE Trans. Automatic Control, Vol. AC-10 (1965), pp. 461-464.

10. Jerzy Seidler, Anatol Badach, Wojciech Molisz. „Metody rozwiązywania zadań optymalizacji”. WNT,Warszawa,1980.

11. Selesnick, I.W., M. Lang, and C.S. Burrus, “Constrained Least Square Design of FIR Filters without Specified Transition Bands,” Proceedings of the IEEE Int. Conf. Acoust., Speech, Signal Processing, Vol. 2 (May 1995),pp. 1260-1263.

12. Selesnick, I.W., M. Lang, and C.S. Burrus. “Constrained Least Square Design of FIR Filters without Specified Transition Bands.” IEEE Transactions on Signal Processing, Vol. 44, No. 8 (August 1996).

13. „Signal Processing Toolbox”. The Mathworks Inc. 2000

14. „Filter Design Toolbox”. The Mathworks Inc. 2000

15. „TMS320C6000 Programmer's Guide”. Texas Instruments, 2001

16. „TMS320C6000 Optimizing C Compiler User's Guide”. Texas Instruments, 2001

17. strona internetowa S.K. Mitry

97

v[n]

y[n]

e[n]

a

x[n]

e[n]

a

+

v[n]

x[n]

+

+

z^-1




























---