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łaszczyz
nie 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łaszczyz
nie 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łaszczyz
nie 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
�
&!
Zależność pomiędzy a wyznacza się z definicji przekształcenia:
j� / 2
�#
e - e- j� / 2 ś#
ś# ź#
j� / 2 2 2 j ź#
ś#
�#1- �#1 �# ś#
2 e- j� ś# 2 - e- j� / 2e- j� / 2 ś# 2 e- j� / 2(e - e- j� / 2 )
= j
ś# ź# ś# ź# ś# ź#
j&! = = =
Td ś# j� / 2 j� / 2 ź#
e + e-
Td ś#1+ j� ź# Td ś# Td ś# j� / 2 j� / 2 j� / 2 ź#
e- # 1 + e- j� / 2e- j� / 2 ź# e- (e + e- ) ś# ź#
�# �# # �# #
ś# ź#
2
�# #
�# ś#
j2 sin(� / 2) 2
= ś# ź# = j tan(� / 2)
ź#
Td ś# cos(� / 2) Td
�# #
2
&! = tan(� / 2)
(*), � = 2arctan(&!Td / 2)
Td
Ze względu na nieliniowość transformacji biliniowej projekt filtra analogowo powinien uwzględnić
�
&!
(*) (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łaszczyz
nie 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 j�c n
c
hlp[n] = = (e - e- j�cn ) = =
�#
-�c
+"e d� = [e ]�
2Ą 2Ąjn 2Ąjn Ą n
�#�c
-�c
, n = 0
�#
Ą
�#
Idealny filtr dolnoprzepustowy (przyczynowy)
�c
ż#e- j�M / 2, | � |< �c
1 sin(�c (n - M / 2))
�# - j�M / 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