Transcript pptx

Signalbehandling og matematik 2
(Tidsdiskrete signaler og systemer)
Session 12.
Design of digital IIR filters
Ved Samuel Schmidt
[email protected]
http://www.hst.aau.dk/~sschmidt/Mat1/
1
IIR og FIR filtre
• IIR
– Systemer med uendelige impuls respons har altid
mindst en betydende pol (det vil sige ikke nul poler eller ophævede poler)
H ( z) 
• FIR
1
1  az 1
– Systemer med endelige impuls respons har ingen
betydende poler (det vil sige ikke nul poler eller ophævede poler)
M
Eksempel:
General form:
H ( z )   bk z  k
k 0
M
H ( z)  1  bz1
Invers transformation:
h(n)   bk x[n  k ]
k 0
IIR vs FIR filter
• Hvorfor IIR filtre?
– IIR filter har stejlere ”sidelobes” end et FIR filter med
samme antal koefficienter.
– Dermed hurtigere og mindre hukommelses krævende
• Hvorfor ikke ?
– Et IIR filter har ikke lineær fase
– Et IIR filter kan være ustabilt
– ET IIR filter er mere sensitivt i for hold til
afrundingsfejl
Definitioner på filter
Navne af frekevnes variabler i
kontinuær og diskret domæne
Kontinuer signaler
Diskrete signaler
Frekvens:
F (Hz)
f=F/Fs (Normaliseret frekvens)
Vinkel hastighed :
Ω=2πF (Radianer / sekund)
ω=2πf (Radianer / sample)
Konvertering
Ω=ω/T
Ω=ΩT
Afgrænsning
-∞<Ω<∞
-π/T<ω<π/T
T= samplings perioden
Design af digitale IIR filtre ved hjælp af
IIR analoge filtre
Specifikation af filteret i digitalt
domæne
Konverter specifikationer til
analogt
Design filteret i det analoge domæne
Konverter det analoge filter
til det digitale domæne
Implementer filteret i det digital
domæne
Stabile systemer in z og s domænet
• Z: Poler skal være i
enhedscirklen
s: Poler skal være i
venstre halvdel
Im
1
1
*1/3
*1/2 1
*1/3
Re
7
jΩ
*1/2 1
σ
3 metoder til konvertering af analoge
IIR filtre til digitale IIR filtre
• Approksimation af afledte
• Impuls invarians
• Bilineær Transformation
Bilineær Transformation
• Ved den Bilineær Transformation sættes s lig
med
2 1  z 1
s
T 1  z 1
• Det betyder at:
 2  1  z 1  

H ( z )  H a  
1  
 T 1 z 
Karakteristika ved Bilineær
Transformation
• Poler fra højre side i s-planet ligger udenfor enhedscirkelen.
• Modsat poler fra venstre side i s-planet som ligger indenfor
enhedscirklen
Im
1
1
1
Re
10
jΩ
1 σ
Karakteristika ved Bilineær
Transformation (3)
• Sammenhæng mellem ΩT og ω
2
2r sin

,
2
T 1  r  2r cos( )
hvis r  1

2 sin
2
 
 tan 
T 1  cos( ) T
2

  2 tan1
T
2
Karakteristika ved Bilineær
Transformation (4)
Im
1
1
1
Re
12
jΩ
1 σ
Agend
Design of Digital IIR filters
• Design af lowpass filtre
– Design af digitale Butterworth filtre
– Design af digitale Chebyshev filtre
– Design af digitale Elliptic filtre
Analogt lavpas Butterworth filter
• Er et ”all pole” filter
– Kvadreret frekvens amplitude respons
H () 
2
•
•
•
•
1
2N
1   / c 
2
eller
N:filter orden
Ωc: 3dB knæk frekvens
Ωp: Anden knæk frekvens
ε: relateret til dæmpning
ved knæk frekvens se figur.
• Laplace transformation
H ( s ) H ( s) 

1
1   s 2 / c2
H () 

N
1
2N
1   2  /  p 
Poler fra analogt Butterworth filter

1  s / 
2

2 N
c
H ( s ) H ( s) 
0

1
1   s 2 / c2

N

 s2
j ( 2 k 1) N
1/ N
 (1)  e
2
c
jΩ

s2
j j ( 2 k 1) 2N
 ce e
2
c
Ωc
σ

j 2
s pk   c e e
j ( 2 k 1) 2N
Polerne vil ligge spejlet omkring både den imaginære
akse og den reelle akse.
Design af digitalt lavpass Butterworth
filter
• Steps:
1. Konverter specifikationer fra digitale til analoge
specifikationer (ω til Ω) Formel 7.26
2. Beregn nødvendig filter orden. Formel 7.293. Dan overførselsfunktionen H(s) ud fra poler i venstre
halvplan
4. Bestem Gain ved Ω=0
5. Transformer til z domænet med bilineær
transformation. Formel 7.18
6. Omskriv til simple rationel form (bk og ak
koefficienter)
Konverter specifikationer fra digitale til
analoge specifikationer (Step 1)
• Hvis specifikationen er
opgivet i Hz find den
normaliserede vinkel
hastighed
• ωc=2πFc/Fs
• Omdan knækfrekvensen ωc
til Ωc

2
 
tan 
T
2
• Samme procedure for andre
vinkel hastigheder. F.eks.
– Stopbånds frekvenser (ωs )
Bestem filter orden (Step 2)
Ωc kendt
Ω =3 dB cut off frequency
c
• Hvis filter orden (N) er forud defineret forsæt
til næste slide.
• Hvis en given dæmpning (δ2) er påkrævet ved
Ωs og Ωc findes N ved:
H () 
2
1
2N
1   / c 
1
  22
2N
1   s /  c 

1

2
2
δ2: dæmpning i stopbånd et ved Ωc
 1   s /  c 
2N

 1

log 2  1
 2
 N
2 log s /  c 
Obs: Hvis N ikke er et hel tal rundes op
Bestem filter orden (Step 2)
Ωc ukendt
•
•
•
•
Ωs: Stop bånd knæk frekvens
Ωp: Pass bånd knæk frekvens
δ1: dæmpning i stopbånd et ved Ωp
δ2: dæmpning i stopbånd et ved Ωs
H ( p ) 
1
2




p
2N
1   p /  c 
H ( s ) 
1
2




c
2N
1   s /  c 
2
2
•
H () 
2
1
2N
1   / c 
Løs for Ωc
log
1
𝛿𝑝
2
−1 /
•
𝑁=
•
Ωc kan findes fra
1 2
−1
𝛿𝑐
2 log Ω𝑠 Ω𝑝
Obs: Hvis N ikke er et hel tal rundes op
2
H ( p ) 
1
2
  p 
2N
1   p /  c 
Bestem overførselsfunktionen fra
poler i venstre side af s-planet (Step 3)
• Beregn poler:
j 2
s pk   c e e
j ( 2 k 1) 2N
k  0,1,2....N  1
• Opstil system funktion fra poler i venstre halv plan
1
H ( s ) H ( s) 
altså dem fra H(s)
1   s /  
2
H ( s) 
1

N
 (s  s
k 1
pk
)
1
( s  s p1 )(s  s p 2 )(s  s p 3 )(s  s p 4 )  ( s  s pN )
2 N
c
Bestem gain (Step 4)
• Normalt ønskes i gain på 1 ved DC.
• Find G så H(0)=1;
H (0) 
G
1
(0  s p1 )(0  s p 2 )(0  s p 3 )(0  s p 4 )  (0  s pN )

G  ( s p1 )( s p 2 )( s p 3 )( s p 4 )  ( s pN )
Fra s domæne til z domæme (Step 5)
• Brug bilinear transformation
H (s) 
1

N
 (s  s
pk
)
1
s  s p1 s  s p 2 s  s p3 s  s p 4 s  s pN 
k 1
2 1  z 1
Bilinear transformation s 
T 1  z 1
H ( z) 
1
1

N
2 1  z 1
 2 1  z 1
 2 1  z 1
 2 1  z 1
 2 1  z 1
  2 1  z 1

(
 sk ) 










s

s

s

s


s

1
1
2
3
4
N
1
 T 1  z 1
 T 1  z 1
 T 1  z 1
  T 1  z 1
T
1

z
T
1

z
k 1




 

Fra kompleks z transformation til
tidsdomæne filter (Step 6)
H ( z) 
1
1

N
2 1  z 1
 2 1  z 1
 2 1  z 1
 2 1  z 1
 2 1  z 1
  2 1  z 1

(
 sk ) 
 s1 
 s2 
 s3 
 s4 
 sN 

1
1
1
1
1
1
k 1 T 1  z
 T 1 z
 T 1  z
 T 1  z
 T 1  z
  T 1 z

Simplificer
N
H ( z) 
b z
k
a z
k
k 0
N
k 0
k
b0  b1 z 1  b2 z  2  b3 z 3 bN z  N

1  a1 z 1  a2 z  2  a3 z 3  a N z  N
k
Invers z-transformation
y[n]  a1 y[n 1]  a2 y[n  2]  a3 y[n  3]  aN y[n  N ]  b0 x[n]  b1x[n 1]  b2 x[n  2]  b3 x[n  3]  bN x[n  N ]
Eksempel
• Konstruer et digitalt butterworth lavpas filter
(fc=40 Hz, Fs=200 sps. og δ2=20dB dæmpning
ved fs=60Hz)
Eksempel step 1
• Normaliserede vinkel hastighed:
– Knækfrekevns
ωc=2π 40/200=0.4 π
ωc=2πFc/Fs
– Stopbåndets hjørne frekvens
ωs=2π 60/200=0.6 π
• Fra digital til analog vinkel hastighed (T=1)
2
 0.4
tan
T
 2

  1.4531rad / s

2
 0.6 
 s  tan
  2.7528rad / s
T
 2 
c 

2
 
tan 
T
2
Eksempel: Bestem filter orden (Step 2)
• Ønsket dæmpning ved Ωs 20dB
c  1.4531rad / s
 s  2.7528rad / s
 2  20dB  10
20
20
 0.1
 1

log
 1
 0.01 
N
 3.6  4
2 log2.753/1.4531
Filter orden N=4
 1

log 2  1
 2
 N
2 log s / c 
Bestem overførselsfunktionen fra
poler i venstre side af s-planet (Step 3)
j 2
• Beregn poler: s pk   c e e
j 2
s pk  1.4531 e e
j ( 2 k 1) 8
j ( 2 k 1) 2N
k  0,1,2....N  1
k  0,1,2,3
S-Plane

2
 (2k  1) 8
0.5
j
• Absolut værdi:
• Vinkeler:
1
1.4531

5

(
2
*
0

1
)

2
8
8
0


2
 (2 *1  1) 8 
7
8
-0.5
-1

-1.5
-1
-0.5
0

0.5
1
1.5
Bestem overførselsfunktionen fra
poler i venstre side af s-planet
(Step 3 forsat)
Opstil system funktion fra poler i venstre halvplan
S-Plane
4 poler i venstre halvplan:
4 poler i venstre halvplan:
0.5
j
(0.556072 j1.342475)
(1.322475 j 0.556072)
1
0
-0.5
-1
-1.5
-1
-0.5
0
0.5
1
1.5

G
H ( s) 
( s  (0.556072 j1.342475))(s  (0.556072 j1.342475))(s  (1.322475 j 0.556072))(s  (1.322475 j 0.556072))
Tip: multiplikation af kompleks konjugerede
s  (a  jb)s  (a  jb  s2 - 2 a s + b2 + a 2
H ( s) 
G
(s 2 + 1.11214s + 2.11145)(s 2 + 2.68495s+ 2.11145)
Bestem gain (Step 4)
• Find G så H(0)=1;
H (0) 
G
1
( 2.11145)( 2.11145)

G  4.4582
H (s) 
4.4582
(s 2 + 1.11214 s + 2.11145)(s 2 + 2.68495s + 2.11145)
Fra s domæne til z domæme (Step 5)
• Bilinear transformation
H ( s) 
4.4582
(s 2 + 1.11214s + 2.11145)(s 2 + 2.68495s+ 2.11145)

2 1  z 1 2 1  z 1
s

1
T 1 z
1  z 1
Bilinear transformation
H ( z) 

  2 1  z 1

  1  z 1

 
2

 2 1  z 1
 + 1.11214 1  z 1



4.4582
  2 1  z 1

 + 2.11145 
  1  z 1



 
2

 2 1  z 1
 + 2.68495 1  z 1



 + 2.11145




Fra kompleks z transformation til
tidsdomæne filter (Step 6)
H ( z) 
H ( z) 
H ( z) 

  2 1  z 1

  1  z 1

21  z
4  8z
1
-1

2
 

2
 2 1  z 1

+
1.11214
1


 1 z



4.4582
  2 1  z 1

 + 2.11145 
  1  z 1






 4z-2 + 1.11214 2  2 z 2


2
 2 1  z 1

+
2.68495
1


 1 z
 
 21  z 
4.45821  z 1
+ 1.11214 2 1  z 1 1  z 1 + 2.111451  z 1
 
2
2






+ 2.6849512 1  z 1 1  z 1 + 2.111451  z 1


2
4
4.45821  z 1
+ 2.111451  2 z 1  z 2 4  8z-1  4z-2 + 2.6849512  2 z 2 + 2.111451  2 z 1  z 2








4
4.45821  z 1
H ( z) 
8.3357420 3.777088z-1  3.887169z-2 11.48136 3.77709z 1  0.741554z 2



4
1


 + 2.11145



4
4.45821  z 1
H ( z) 
95.706- 74.851z-1  65.078z-2 - 17.483z-3 + 2.8825z-4


4
0.0104* 4.45821  z 1
H ( z) 
-4
1 - 0.7820923z -1  0.679976z-2 - 0.182675z-3 + 0.0301187z
0.04636 0.18546z-1 + 0.27819z-2 + 0.18546z-3 + 0.04636z-4
H ( z) 
1 - 0.78209z -1  0.67997z-2 - 0.18267z-3 + 0.03011z-4
Mål:
b0  b1 z 1  b2 z 2  b3 z 3 bN z  N
H ( z) 
1  a1 z 1  a2 z 2  a3 z 3 aN z  N
Test af eksempel
• Test i matlab
b=[0.0464
0.1855
a=[1.0000
-0.7821
freqz(b,a,1000,200)
0.2782
0.6800
0.1855
-0.1827
0.0464];
0.0301];
• Sammenlign med ”butter” i matlab
[b a]=butter(4,[40/100])
b =[0.0466
0.1863
0.2795
a =[1.0000
-0.7821
0.6800
0.1863
-0.1827
0.0466]
0.0301]
Chebyshev filter type I
• Overførselsfunktion
H  
2
1
1   2TN2  /  p 
• Hvor ε er relateret til
ripples i pasbåndet
• Hvor TN er et N ordens
polynomium
 cos(N cos1 x) x  1
TN ( x)  
1
cosh(N cosh x) x  1
Chebyshev filter type I
Poler
• Polerne ligger på en
ellipse
 2 1
r1   p
2
 2 1
r2   p
2
Hvor β er relateret til ε
1/ N
 1   2  1
 




• Polernes location:
xk  r2 cos(k )
yk  r1 sin(k )
Hvor vinklen φk er:
k  2  (2k  1) N
Chebyshev filter type II
• Overførselsfunktion
H 
•
2
1

1   2 TN2  s /  p / TN2  s / 

Bemærk indeholder også
nulpunkter
• Hvor ε er relateret til
ripples i stopbåndet
• Hvor TN er samme
Chebyshev
polynomium
 cos(N cos1 x) x  1
TN ( x)  
1
cosh(N cosh x) x  1

Bestem filter orden (N) på Chebyshev
filter
N
• Hvor
• N: filter orden
• ε: Ripple i pass bånd
• δ2 : Dæmpning i stopbåndet
• Ωs : Knæk frekvens
• Ωp: Start på stopbånd
• Hvor TN er samme Chebyshev
polynomium


 
  1
log 1   22  1   22 1   2 /  2
log  s /  p  


s / p
2
Transformation
af lavpas filtre til andre filter typer.
• Transformer et eksisterende filter til ønskede
egenskaber.
1
0.5
4
0
Simple transformation
Spejl poler omkring IMG.
aksen
-0.5
-1
-1
0
Real Part
Imaginary Part
Imaginary Part
1
-0.5
-1
-1
-0.5
0
0.5
Real Part
1
50
Magnitude (dB)
Magnitude (dB)
4
0
1
0
-50
-100
-150
-200
-250
0.5
0
0.2
0.4
0.6
0.8
1
Normalized Frequency ( rad/sample)
0
-50
-100
0
0.2
0.4
0.6
0.8
1
Normalized Frequency ( rad/sample)
Transformation i analogt domæne
Filter type
Lavpas>Lavpas
Transformation
p
s
s
' p
Lavpas>Højpas
 p ' p
Lavpas>Båndpas
Lavpas>Stopbånd
p
Oprindelig knæk frekvens
' p Ny knæk frekvens
s
s
s 2  l u
s  p
s ( u   l )
s ( u   l )
s  p 2
s  l u
Ny knæk frekvens
' p
' p
l ,  u
l ,  u
l :
Laveste knækfrekvens
u :
Højeste knækfrekvens
Eksempel
• Konverter vores lavpas filter til et højpas filter
med knæk frekvens ved ωc=0.5 π
H lp ( s) 
s
4.4582
(s 2 + 1.11214s + 2.11145)(s 2 + 2.68495s+ 2.11145)
 p ' p
 p  2.7528rad / s
2
 0.5 
' p  tan
  2 rad / s
T
 2 
s
5.5056
s
s
H hp ( s) 
4.4582
  5.5056
  5.5056 2

 5.5056
 5.5056

 
+
1.11214
+
2.11145
+ 2.68495
+ 2.11145





 s 
  s 

 s 
 s 



2

H hp ( s) 
4.4582s 4
4.4582s 4

2.11145s2  6.1244s 30.3116 2.11145s2 + 14.7550s 30.3116 4.45822s4 + 44.0858s3 + 218.362s2 + 632.855s+ 918.696



Test af eksempel i Matlab
4.4582s 4
H hp ( s) 
4.45822s4 + 44.0858s3 + 218.362s2 + 632.855s+ 918.696
0
44.0858
0
218.3620
0
632.8550
0
10
-5
10
Magnitude
b=[4.4582
a=[4.4582
freqs(b,a)
-10
10
-15
10
-2
10
-1
10
0
10
Frequency (rad/s)
1
10
2
10
0];
918.6960];
Transformation i det digitale domæne
Eksempel digitalt
• Konverter vores lavpas filter til et højpas filter med knæk
frekvens ved ωc=0.5 π
-4
0.046365 0.18546z-1 + 0.27819z-2 + 0.18546z-3 + 0.046365z
H ( z) 
-3
-4
1 - 0.7820923z -1  0.679976z-2 - 0.182675z
+ 0.0301187z
• Oprindelig knækfrekvens ωp=0.4π
• Ny knækfrekvens ω’p=0.5π
z 1  a
1
• Transformation: z  
, hvor
1
1  az
a
cos p  ' p / 2
cos p  ' p / 2
2
 0.1584
3
 z 1  0.1584 
 z 1  0.1584 
 z 1  0.1584 
 z 1  0.1584 







0.046365 0.18546 
+ 0.27819 
+ 0.18546 
+ 0.046365 
1 
1 
1 
1 
1

0
.
1584
z
1

0
.
1584
z
1

0
.
1584
z
1

0
.
1584
z








H ( z) 
2
3
4
1
1
1
1
 z  0.1584 
 z  0.1584 
 z  0.1584 
 z  0.1584 
  0.679976 
 - 0.182675 
 + 0.0301187 

1 - 0.7820923 
1 
1 
1 
1 
 1  0.1584z 
 1  0.1584z 
 1  0.1584z 
 1  0.1584z 
H ( z) 



3


2

2




 0.1584 1  0.1584z  + 0.0301187 z
3

0.046365 0.18546 z 1  0.1584 1  0.1584z 1 + 0.27819 z 1  0.1584 1  0.1584z 1 + 0.18546 z 1  0.1584 1  0.1584z 1 + 0.046365 z 1  0.1584


1 - 0.7820923 z 1  0.1584 1  0.1584z
  0.679976 z
1 3
H ( z) 
1

2
 0.1584 1  0.1584z

1 2

- 0.182675 z 1
0.0754 - 0.0801 z -1  0.1131z -2 - 0.1045z -3 - 0.0197z -4
1.0000  0.9342z -1 1.0437z -2 0.4071z -3 0.0419z -4
3
1
1

 0.1584
4
4
4
Test af digitalt eksempel i Matlab
0.0754 - 0.0801z -1  0.1131z-2 - 0.1045z-3 - 0.0197z-4
H ( z) 
1.0000  0.9342z-1 1.0437z-1 0.4071z-3 0.0419z-4
Magnitude (dB)
b=[0.0754
-0.0801
a=[1.0000
0.9342
freqz(b,a,1000,200)
0.1131
1.0437
-0.1045
0.4071
-0.0197];
0.0419];
0
Ups: lille fejl
-20
-40
-60
0
10
20
30
40
50
60
Frequency (Hz)
70
80
90
100
0
10
20
30
40
50
60
Frequency (Hz)
70
80
90
100
Phase (degrees)
200
0
-200
-400
Sammenligning mellem filtre
• 4. ordens filter (fc=40 Hz, fs =80 Hz ,Fs=200 sps)
Magnitude Response (dB)
0
Butter
-10
Optimal filter
-20
Chev 1
Elliptic
Chev 2
Magnitude (dB)
-30
-40
-50
-60
-70
-80
-90
-100
0
10
20
30
40
50
Frequency (Hz)
60
70
80
90
Update til næste år
• Mere fokus på transformation mellem filter
typer
• Mere fokus på Butterworth poler
• Ny opgaver mindre bi linear