Projektowanie filtrów cyfrowych Projektowanie filtra cyfrowego przebiega w następujących etapach: 1. Określenie właściwości filtra (wymagań), 2. Aproksymacja zadanych właściwości, 3. Realizacja systemu. p = 0.4 s = 0.6 rp = 6 [dB] rs = 20 [dB] Dodatkowym wymaganiem może być np. liniowa faza lub niski rząd filtra. 1 Projektowanie dyskretnych filtrów IIR na podstawie filtrów analogowych Filtry analogowe Filtr Butterwortha | H ( j&!) | Charakterystyka częstotliwościowa filtra Butterwortha jest maksymalnie płaska w paśmie przepustowym. Oznacza to, że dla filtra rzędu N pierwsze (2N-1) pochodnych funkcji | H( j&!) |2 ma &! = 0 wartość zero dla . Charakterystyka amplitudowa jest monotoniczna w paśmie przepustowym i zaporowym. Kwadrat charakterystyki amplitudowej definiuje się następująco: 1 | Hc ( j&!) |2 = 1+ ( j&! / j&!c )2N 2 s = j&! Podstawiając 1 | Hc ( j&!) |2 = Hc (s)Hc (-s) = 1+ (s / j&!c )2N Bieguny transmitancji spełniają 1+ (s / j&!c )2N = 0 , skąd j(+k 2Ą ) / n j j( +k 2Ą ) n n C = M e , k = 0, 1,...n -1 C = Me = Me ( , ) j(Ą / 2N )(2k+N-1) 2N sk = j -1&!c = &!ce , k = 0, 1,...2N -1 (*) Hc (s)Hc (-s) Ponieważ (*) podaje bieguny iloczynu do realizacji filtra wybiera się bieguny z lewej półpłaszczyzny, co zapewnia stabilność filtra analogowego. Przykłady: s1,2 = -0.7071 ą j0.7071 1 Hc (s) = (s + 0.7071 + j0.7071)(s + 0.70711 - j0.7071) 1 Hc (s) = s2 + 1.4142s + 1 3 s1,3 = -0.5 ą j0.866, s2 = 1 1 Hc (s) = s3 + 2s2 + 2s +1 s1,4 = -0.7654 ą j1.8478 s2,3 = -1.8478 ą j0.7654 N w = &!c = 24 = 16 16 Hc (s) = s4 + 5.2263s3 + 13.6569s2 + 20.905s + 16 4 Przykład &!c Obliczyć rząd N i pulsację (pulsacja 3dB) analogowego filtra Butterwortha spełniającego następujące &! = 100 rad / s &!s = 160 rad / s wymagania: krawędz pasma przepustowego p , krawędz pasma zaporowego , rp = 1 dB maksymalne tłumienie w paśmie przepustowym , minimalne tłumienie w paśmie zaporowym rs = 30 dB . 5 1 Dzieląc stronami i logarytmując: Ze wzoru | Hc ( j&!) |2 = j&! / j&!c )2N 1+ ( 2N p otrzymujemy: # &! # ś# 10r /10 -1ś# ź# log10 ś# p ź# = log10 ś# 1 ż# ś# ź# ś# s ź# &!s 10r /10 -1 = -rp [dB] # # # # # 1 + (&! / &!c )2N # p # 1 # otrzymujemy: = -r [dB] s # 1 + (&!s / &!c )2N # p # # ś# 1 10r /10 -1ś# # &! ś# 1 101/10 -1 100 ź# ź# N = log10ś# / log10ś# p ź# = log10ś# / log10# ś# = 8.785 ś# ź# ś# ź# ś# ś# s ź# 2 &!s 2 ż#20log + (&! / &!c )2N -1/ 2 10r /10 -1 1030/10 -1ź# #160 # # # # # (1 ) = -rp # # # 10 p # #20log10 + (&!s / &!c )2N -1/ 2 s (1 ) = -r # 1 p log10(&! / &!c ) = log10(10r /10 -1) p 2N ż#log10 + (&! p / &!c )2N rp /10 (1 )= # # 1 #log10(1 + (&!s / &!c )2N )=r /10 p # s log10(&!p / &!c) = log10(10r /10 -1)2N ż#(&! p / &!c )2N =10rp /10 -1 1 # # p &!p /&!c = (10r /10 -1)2N #(&!s / &!c )2N =10rs/10 -1 # &! 100 p &!c = = = 107.7957 1 1 p (10r /10 -1)2N (101/10 -1)2"9 Ponieważ rząd filtra musi być liczbą całkowitą N zostaje zaokrąglone w górę N=9. 6 Transformacja częstotliwości - filtry FGP, FBP, FBS 7 8 Filtry Czebyszewa Kwadrat charakterystyki amplitudowej filtra Czebyszewa typu I definiuje się następująco: 1 | Hc ( j&!) |2 = 2 2 1+ VN ( j&!/ j&!c ) ż#cos(N cos-1 x), | x |d" 1 # cos-1 x = arccos(x) gdzie VN (x) = # jest wielomianem Czebyszewa rzędu N ( ). # #cosh(N cosh-1 x), | x |> 1 V0 (x) = cos(0arccos( x)) = 1 V1(x) = cos(arccos (x)) = x V2 (x) = cos(2arccos(x)) = 2x2 -1 , , VN +1(x) = 2xVN (x) - VN -1(x), N e" 1 Wielomiany Czebyszewa można wyznaczać rekurencyjnie: 2 0 d" VN (x) d" 1 Dla | x |d" 1 , dla | x |> 1 VN (x) rośnie monotonicznie. VN (1) = 1 VN (0) = ą1 ; dla N parzystego i V (0) = 0 dla N nieparzystego. N 9 Bieguny analogowego filtra Czebyszewa są położone na płaszczyznie zespolonej S na elipsie. a&!c b&!c Długość półosi małej wynosi a półosi wielkiej , przy czym: 1 1 -1/ N -1/ N -1 -2 a = (ą1/ N -ą ) b = (ą1/ N + ą ) ą = + 1 + , , . 2 2 Na okręgu o promieniu równym wielkiej półosi elipsy wyznacza się bieguny filtra Butterwortha o tym samym rzędzie, co projektowany filtr Czebyszewa, a następnie przesuwa się je na obwód elipsy. s1,2 = -0.1490 ą j0.9037 s3 = -0.2980 s1,2 = -0.0850 ą j0.9464 s3,4 = -0.2052 ą j0.3920 , , 10 Kwadrat charakterystyki amplitudowej filtra Czebyszewa typu II definiuje się następująco: 1 | Hc ( j&!) |2 = 1 1+ 2 2 VN ( j&!c / j&!) 2 2 VN ( j&!/ j&!c ) W porównaniu z filtem typu I następuje odwrócenie członu w mianowniku oraz j&!c / j&! odwrócenie argumentu wielomianu Czebyszewa na . Bieguny filtra typu II można wyznaczyć znając bieguny filtra typu I na podstawie zależności: 2 &!c II sk = I sk natomiast zera filtra typu II leżą na osi urojonej w punktach V (&!c / &!) = 0 . N 11 s1,2 = -0.1777 ą j1.0773 z1,2 = ą j1.1547 s3 = -3.3553 , ; s1,2 = -0.2860 ą j0.8213 z1,2 = ą j1.1547 s3 = -1.3221 , ; 12 13 Filtr eliptyczny Charakterystyka amplitudowa filtra eliptycznego jest równomiernie falista w paśmie przepustowym &!c i zaporowym. Filtr opisują cztery parametry: N - rząd, - krawędz pasma przepustowego oraz Rp - tłumienie w paśmie przepustowym i Rs - tłumienie w paśmie zaporowym. 14 15 16 17 Transformacja biliniowa j&! Transformacja biliniowa jest przekształceniem algebraicznym odwzorowującym oś urojoną płaszczyzny S w okrąg jednostkowy na płaszczyznie Z. Ponieważ - " d" &! d" " -Ą d" d" Ą , a więc przekształcenie to jest nieliniowe. Transmitancję dyskretną H(z) uzyskuje się z transmitancji analogowej Hc(s) przez podstawienie: #1- z-1 ś# 2 Ą# #1- ś# ź# 2 z-1 ś#ń# s = ś# ź#Ą# Td ś#1+ ź# z-1 # (*), tzn.: H (z) = Hc ó# d ś#1+ z-1 ź#Ą# ó#T # # # Ł# Ś# s = + j&! Rozwiązując (*) względem z i podstawiając otrzymujemy: 1 + (Td / 2)s 1 + Td / 2 + j&!Td / 2 z = = 1 - (Td / 2)s 1 - Td / 2 - j&!Td / 2 | z |< 1 | z |>1 < 0 > 0 Dla , a dla . Tak, więc jeżeli biegun transmitancji analogowej Hc(s) leży w lewej półpłaszczyznie to moduł transmitancji H(z) leży wewnątrz okręgu jednostkowego. Stabilny układ analogowy jest przekształcany w stabilny układ dyskretny. 1 + j&!Td / 2 z = s = j&! | z |= 1 j&! Podstawiając do wyrażenia na z widzimy, że dla całej osi 1 - j&!Td / 2 1 + j&!Td / 2 j e = . 1 - j&!Td / 2 18
&! (*) (tzn. filtr analogowy należy zaprojektować na , a nie ). 19 Przykład: Zaprojektować dolnoprzepustowy, cyfrowy filtr eliptyczny spełniający następujące wymagania: częstotliwość próbkowania Fp=100 [Hz], krawędz pasma przepustowego f1=20 [Hz], krawędz pasma zaporowego f2=25 [Hz], nieliniowość charakterystyki w paśmie przepustowym Rp=0.1 [dB], nieliniowość charakterystyki w paśmie zaporowym Rs=60 [dB]. Projekt
1. Unormowanie częstotliwości filtra cyfrowego f1 i f2 podanych w [Hz] do pulsacji cyfrowej w [rad]: f1 20 f2 25 = 2Ą = 2Ą = 1.2566 = 2Ą = 2Ą =1.5708 1 , 2 Fp Fp 100 100 &! 2. Obliczenie krawędzi pasma [rad/s] filtra analogowego uwzględniające nieliniowość transformacji biliniowej: 2 2 &!1 = tan( / 2) = tan(1.2566 / 2) = 145.3085 1 &!2 = 200 , Td 1/100 3. Określenie parametrów filtra analogowego na podstawie podanych wymagań odnośnie pasma oraz H (s) wyznaczenie jego transmitancji w postaci zer i biegunów. Dla zadanych wymagań rząd analogowego filtra eliptycznego wynosi N=7, rozkład jego zer i biegunów przedstawia rysunek. 1 + (Td / 2)s z = 4. Transformacja biliniowa zer i biegunów filtra analogowego H(s) zgodnie ze wzorem . 1 - (Td / 2)s Otrzymujemy zera i bieguny filtra cyfrowego H(z). 20 filtr eliptyczny Uwaga! j&! Transformacja biliniowa odwzorowuje całą oś płaszczyzny S w okrąg jednostkowy na płaszczyznie Z, j j&! = " e = -1 dlatego jeżeli filtr analogowy ma zera w to przechodzą one w zera . Rozważany filtr ma jedno zero w -1 , ponieważ stopień licznika jego transmitancji H(s) jest o jeden mniejszy od stopnia mianownika. (Dla porównania dolnoprzepustowy filtr Butterwortha ma wielokrotne zero w -1 ). filtr Butterwortha 21 5. Weryfikacja charakterystyk amplitudowych filtra cyfrowego 22 23 24 25 Implementacja Matlaba ? 26 Metoda okien - filtry FIR - liniowa faza W metodzie okien najpierw zadaje się wymaganą charakterystykę częstotliwościową filtra: " j j n Hd (e ) = d "h [n]e- , n=- " na podstawie której wyznacza się jego nieskończoną odpowiedz impulsową Ą 1 j j n hd [n] = d +"H (e )e d , 2Ą -Ą a następnie mnoży się tę odpowiedz impulsową przez okno o skończonej długości (wycina się fragment np. oknem prostokątnym): 1, 0 d" n d" M ż# w[n] = #0, | n |> M h[n] = hd[n]w[n] , . # j H(e ) Na podstawie twierdzenia o okienkowaniu jest splotem widma filtra idealnego z widmem okna: Ą 1 j jŚ j(-Ś) H(e ) = )dŚ d +"H (e )W (e . 2Ą -Ą 27 Idealny filtr dolnoprzepustowy (nieprzyczynowy) 1, | |< c ż# j Hlp (e ) = # #0, c <| |d" Ą sin(cn) ż# c , n `" 0 # 1 1 1 sin(cn) Ą n # j n j n jc n c hlp[n] = = (e - e- jcn ) = = # -c +"e d = [e ] 2Ą 2Ąjn 2Ąjn Ą n #c -c , n = 0 # Ą # Idealny filtr dolnoprzepustowy (przyczynowy) c ż#e- jM / 2, | |< c 1 sin(c (n - M / 2)) # - jM / 2 j n j hlp[n] = e d = Hlp (e ) = # +"e , 2Ą Ą (n - M / 2) #0, c <| |d" Ą # -c 28 Można również wyznaczyć współczynniki filtra przez obliczenie odwrotnej dyskretnej transformaty Fouriera (IDFT) zadanej charakterystyki częstotliwościowej. %4. Wykresy figure, close all, clear all subplot(2,1,1) stem(real(Hf)), N = 32; %długość filtra FDP xlabel('k'), ylabel('Re[H(k)]') Wn = 0.5; %częstotliwość graniczna, 1 odpowiada fp/2 subplot(2,1,2) stem(imag(Hf)), axis([0 N -1 1]) xlabel('k'), ylabel('Im[H(k)]') %1. zadana ch-ka częstotliwościowa figure, hold on %1.1 od 0 do fp/2 stem(real(ht),'filled','r') Hf=zeros(1,round(N/2)); Hf(1:round(Wn*N/2))=1; stem(B ,'b'), Hf=Hf.*exp(-j*(1:length(Hf))*pi); %liniowa faza xlabel('n'), ylabel('Re[h(n)]') %1.2 od fp/2 do fp - symetrie widma legend('obliczeniowo','analitycznie') Hf=[1 Hf 0 conj(fliplr(Hf))]; figure subplot(2,1,1), hold on %2. Współczynniki filtra plot(f,abs(H), 'r:'), %2.1 Przez DFT plot(f,abs(Hb),'b') ht=ifft(Hf); ht(1)=[]; xlabel('f unormowana'), ylabel('|H[e^j^\omega]|') %2.2 Dla porównania - analitycznie legend('obliczeniowo','analitycznie') n=1:(length(ht)-1)/2; B=sin(Wn*n*pi)./(pi*n); B=[fliplr(B) Wn B]; subplot(2,1,2), hold on %3. Ch-ki częstotliwościowe filtrów plot(f,unwrap(angle(H)),'r:'), plot(f,unwrap(angle(Hb)),'b'), M=1024; %zmniejszenie kroku w osi częstotliwości xlabel('f unormowana'), ylabel('faza') htz=zeros(1,M); htz(1:length(ht))=ht; legend('obliczeniowo','analitycznie') htb=zeros(1,M); htb(1:length(B)) =B; figure , hold on H=fft(htz); plot(f,20*log10(abs(H)), 'r:'), Hb=fft(htb); plot(f,20*log10(abs(Hb)),'b'), f=2*(0:M-1)/(M-1); %unormowana os częstotliwości xlabel('f unormowana'), ylabel('|H[e^j^\omega]| [dB]') legend('obliczeniowo','analitycznie') axis([0 1 -40 5]) 29 Zadana charakterystyka częstotliwościowa filtra Uzyskana charakterystyka częstotliwościowa filtra Współczynniki filtra 30 31 Filtr FIR z oknem prostokątnym - największa stromość w paśmie przejściowym i najmniejsze tłumienie pierwszego listka bocznego. 32 Zwiększenie długości okna M powoduje zmniejszenie szerokości listka głównego nie wpływa jednak na położenie (tłumienie) pierwszego listka bocznego. 33 34 35 Okno Kaisera - okno parametryczne 1/ 2 ń# ż# I0 Ą#(1-[(n-ą ) /ą ]2) ó# Ą# Ł# Ś# # , 0 d" n d" M w[n] = I0 ( ) # , #0 # ą = M / 2 , I0(") - zmodyfikowana funkcja Bessela zerowego rzędu pierwszego rodzaju. Okno Kaisera ma dwa parametry: - długość (liczba niezerowych współczynników),
- parametr kształtu .
Wzrost powoduje obniżanie pierwszego listka bocznego (i niestety zwiększanie szerokości listka głównego - tj. pasma przejściowego). Natomiast wzrost M powoduje zawężanie listka głównego i nie wpływa na położenie pierwszego listka bocznego.
Na podstawie doświadczeń numerycznych Kaiser określił wymagania odnośnie M i pozwalające " spełnić zadane wymagania charakterystyki amplitudowej filtra FIR (tj. , ). 36 " = s - p A = -20 log10 0.1102(A - 8.7), A > 50 ż# # = #0.5842(A - 21)0.4 + 0.07886(A - 21), 21 d" A d" 50 #0, A < 50 # A - 8 M = ą2 z dokładnością do 2.285" 37 38 39 Przykład: Dobrać parametry okna Kaisera dla FDP: krawędz pasma przepustowego p = 0.4Ą , krawędz rp = 0.1 dB pasma zaporowego s = 0.6Ą , maksymalna nieliniowość w paśmie przepustowym , rs = 60 dB minimalne tłumienie w paśmie zaporowym .
1. Obliczenie dla obu pasm: p =(1-10-r / 20)/ 2 = 0.0057 p pasmo przepustowe s = 10-rs / 20 = 0.001 pasmo zaporowe p Ponieważ dla FIR projektowanego metodą okien = s wiec wybiera się wymaganie ostrzejsze: = min{ , s} = 0.001 p = 5.6533 M = 37 2. Korzystając ze wzorów projektowych wyznacza się i , a następnie okno o tych parametrach. p + s 3. Odpowiedz impulsową FIR wylicza się dla = . 2 4. Weryfikacja otrzymanych charakterystyk filtra i ewentualna ich korekcja (przez zmianę parametrów okna). 40 Po korekcji parametrów 41 Filtry FGP, FPP, FPS 42 Optymalny filtr FIR - algorytm Parksa-McClellana Przy projektowaniu FIR metodą okien, okno prostokątne jest najlepszą średniokwadratową aproksymacją zadanej charakterystyki częstotliwościowej dla danego rzędu M, tzn. minimalizuje wyrażenie: Ą 2 1 2 j j = Hd (e ) - H(e ) d +" 2Ą -Ą Niestety powyższe kryterium nie uwzględnia oscylacji w punktach nieciągłości (efekt Gibbsa) i uniemożliwia osobne traktowanie pasma przepustowego i zaporowego. Dlatego stosuje się optymalizację min-max (minimalizację błędów maksymalnych) oraz ważone częstotliwościowo kryteria błędów. W algorytmie Parksa-McClellana rząd filtra
s p krawędzie pasma i oraz stosunek 1 /2 1 2 są stałe natomiast (lub ) jest zmienne. 43 44 45