Transcript Moonen - F A R A D A Y

Digital Signal Processing
FIR Filter Design
Marc Moonen
Dept. E.E./ESAT, K.U.Leuven
DSP-II
p. 1
FIR Filter Design
• Review of discrete-time systems
LTI systems, impulse response, transfer function, ...
• FIR filters
Direct form, Lattice, Linear-phase filters
• FIR design by optimization
Weighted least-squares design
Minimax design
• FIR design using `windows’
• Equiripple design
• Software (Matlab,…)
Review of discrete-time systems 1/10
Linear time-invariant (LTI) systems
u[k]
y[k]
• Linear :
input u1[k] -> output y1[k]
input u2[k] -> output y2[k]
hence input a.u1[k]+b.u2[k]-> a.y1[k]+b.y2[k]
• Time-invariant (shift-invariant)
input u[k] -> output y[k]
hence input u[k-T] -> output y[k-T]
Review of discrete-time systems 2/10
Linear time-invariant (LTI) systems
u[k]
• Causal systems:
for all input u[k]=0, k<0 -> output y[k]=0, k<0
• Impulse response :
input 1,0,0,0,... -> output h[0],h[1],h[2],h[3],...
input u[0],u[1],u[2],u[3] -> output y[0],y[1],y[2],y[3],...
y[k ]   u[i].h[k  i]  u[k ] * h[k ]
i
= `convolution’
y[k]
Review of discrete-time systems 3/10
• Impulse response/convolution:
u[0],u[1],u[2],u[3],0,0….
y[0],y[1],...
y[k ]   u[i].h[k  i]  u[k ] * h[k ]
i
0
0 
 y[0]  h[0] 0
 y[1]   h[1] h[0] 0
 u[0]
0


 

 y[2] h[2] h[1] h[0] 0   u[1] 



.
h[2] h[1] h[0] u[2]
 y[3]  0


 y[4]  0
0
h[2] h[1]  u[3]

 

0
0
h[2]
 y[5]  0
h[0],h[1],h[2],0,0,...
`Toeplitz’ matrix
Review of discrete-time systems 4/10
• Z-Transform:
H ( z)   h[i].z i U ( z )   u[i].z i Y ( z )   y[i].z i
i
i
i
0
0
0 
 y[0]
 h[0]
 y[1] 
 h[1] h[0]
 u[0]
0
0






 y[ 2]
h[2] h[1] h[0]
0   u[1] 
1
2
3
4
5
1
2
3
4
5
1 z
z
z
z
z .
z
z
z
z .
1 z
.
y
[
3
]
0
h
[
2
]
h
[
1
]
h
[
0
]



 u[ 2]


 y[ 4]
 0
0
h[2] h[1]  u[3]




0
0
h[2]
 y[5]
 0
 


Y ( z)


H ( z).1 z1





Y ( z )  H ( z ).U ( z )
z2
z3





Review of discrete-time systems 5/10
Z-Transform :
• input-output relation
• `shorthand’ notation
Y ( z )  H ( z ).U ( z )
(for convolution operation/Toeplitz-vector product)
• stability
bounded input u[k] -> bounded output y[k]
--iff
h[k ]  

k
--iff poles of H(z) inside the unit circle
(for causal,rational systems)
Review of discrete-time systems 6/10
Frequency response :
• given a system with impulse response h[k]
• given an input signal = complex sinusoid
• output signal :
u[k ]  e
jk
,  k  
y[k ]   h[i].x[k  i]   h[i].e j ( k i ) e jk  h[i].e ji x[k ].H ( z) z e j
i
i
j
H (e )
i
is `frequency response’
is H(z) evaluated on the unit circle
Review of discrete-time systems 7/10
Frequency response :
• for a real impulse response h[k]
Magnitude response H (e j )
Phase response
H (e j )
• example :
1
Nyquist frequency
0.5
0
-4
-2
0
2
-2
0
2
5

e jk  ...,1,1,1,1,1,...
4
0
-5
-4
is even function
is odd function
4
Review of discrete-time systems 8/10
`Popular’ frequency responses for filter design :
low-pass (LP)
high-pass (HP)
band-pass (BP)



band-stop

multi-band

Review of discrete-time systems 9/10
Rational transfer functions (`IIR filters’):
•
1
N
B( z ) b0  b1 z  ...  bN z
H ( z) 

1
N
A( z ) 1  a1 z  ...  aN z
• N poles (zeros of A(z)) , N zeros (zeros of B(z))
• corresponds to difference equation
y[k ]  a1. y[k 1]  ... aN . y[k  N ]  b0 .u[k ]  b1.u[k 1]  ... bN .u[k  N ]
y[k ]  b0 .u[k ]  b1.u[k 1]  ... bN .u[k  N ]  a1. y[k 1]  ... aN . y[k  N ]
Review of discrete-time systems 10/10
`FIR filters’ (finite impulse response):
•
1
H ( z)  B( z)  b0  b1z  ... bN z  N
•
•
•
•
`Moving average filters’ (MA)
N poles at the origin z=0 (hence guaranteed stability)
N zeros (zeros of B(z)), `all zero’ filters
corresponds to difference equation
y[k ]  b0 .u[k ]  b1.u[k 1]  ...  bN .u[k  N ]
• impulse response
h[0]  b0 , h[1]  b1 ,...,h[ N ]  bN , h[ N  1]  0,...
Review of discrete-time systems
`FIR filter’ (finite impulse response) design
• this lecture
+++ : phase control (linear phase)
guaranteed stability
design flexibility
minor coefficient sensitivity/quantization/round-off
problems,….
- - - : long filters, hence expensive
FIR Filters 1/14
1
H ( z)  B( z)  b0  b1z  ... bN z
N
y[k ]  b0 .u[k ]  b1.u[k 1]  ...  bN .u[k  N ]
u[k]
• `direct-form’
realization
u[k-1] u[k-2] u[k-3]

bo

b1

b2
u[k-4]

b3
b4
x
x
x
x
+
+
+
+
y[k]
x
FIR Filters 2/14
• `retiming’ = select subgraph (shaded)
remove delay element on all inbound arrows
add delay element on all outbound arrows
u[k]
u[k-1] u[k-2] u[k-3]

bo

b1

b2
u[k-4]

b3
b4
x
x
x
x
+
+
+
+
y[k]
x
FIR Filters 3/14
• `retiming’ : results in...
u[k]
u[k-1]
u[k-2]

bo

b1
b2
x
x
+
+
y[k]

u[k-3]

b3
b4
x
x
+
+
x
FIR Filters 4/14
• `retiming’ : repeated application results in...
`Transposed direct-form’ realization
u[k]
bo
b1
x
y[k]
+
b2
x

+
b3
x

+
b4
x

+
x

FIR Filters 5/14
y[k ]  b0 .u[k ]  b1.u[k 1]  ...  bN .u[k  N ]
• `Lattice form’ : derived from combined realization with
y[k ]  bN .u[k ]  bN 1.u[k 1]  ... b0 .u[k  N ]
reversed coefficient vector results in
- same magnitude response
- different phase response
FIR Filters 6/14
• `Lattice form’ : derivation (I)
u[k]
bo x
y[k]
y[k]

u[k-1]

u[k-2]

u[k-3]
u[k-4]

b1 x
b2 x
b3 x
b4 x
x b3
x b2
x b4
x b1
x bo
+
+
+
+
+
+
+
+
FIR Filters 7/14
• `Lattice form’ : derivation (II), this is equivalent to...
u[k]
b’o x
y[k]
+
x
ko
x
y[k]
+

u[k-1]

u[k-2]

u[k-3]
u[k-4]

b’1 x
b’2 x
b’3 x
0 x
x b’3
x b’2
x 0
x b’1
x b’o
+
+
+
k0 
+
+
+
+
b4 '
, b0  b0 , b1'  ...,b2'  ...,b3'  ...,
b0
+
(=simple proof)
FIR Filters 8/14
• `Lattice form’ : derivation (III), retiming leads to...
u[k]
b’o x
y[k]
+
x
ko
x
y[k]
+


u[k-2]

u[k-2]

u[k-3]
b’1 x
b’2 x
b’3 x
x b’3
x b’2
x b’1
x b’o
+
+
+
+
+
+
…repeat procedure for shaded graph
FIR Filters 9/14
• `Lattice form’ : derivation (IV), end result is...
u[k]
k4
x
y[k]
+
+
x
+
x
ko
+
x
k1
x
y[k]
+
x
k2
x

+
k3
x

+
x

+

FIR Filters 10/14
`Lattice form’ :
• ki’s are `reflection coefficients’
• procedure for computing ki’s from bi’s corresponds to
`Schur-Cohn’ stability test (control theory):
all zeros of B(z) are stable (i.e. lie inside unit circle)
iff all reflection coefficients statisfy |ki|<1 (i=1,…,N-1)
(ps: procedure breaks down if |ki|=1 is encountered)
FIR Filters 11/14
Linear-phase FIR filters
• Non-causal zero-phase filters :
example: symmetric impulse response
h[-L],….h[-1],h[0],h[1],...,h[L]
h[k]=h[-k], k=1..L
L
frequency response is
j
H (e ) 
L
 h[k ].e
k  L
 j .k
k
L
 ...   a[k ].cos(.k )
k 0
- real-valued (=zero-phase) transfer function
- causal implementation by introducing (group) delay
FIR Filters 12/14
Linear-phase FIR filters
• Causal linear-phase filters :
example: symmetric impulse response & N even
h[0],h[1],….,h[N]
N=2L (even)
h[k]=h[N-k], k=0..L
0
N
frequency response is
N
H (e )   h[k ].e
j
k 0
 j .k
 ...  e
 jL
L
. a[k ].cos(.k )
k 0
- causal implementation of zero-phase filter, by
L
 jL

e
introducing (group) delay z
j
z e
k
FIR Filters 13/14
Linear-phase FIR filters
Type-1
N=2L=even
symmetric
h[k]=h[N-k]
e
 jN / 2
Type-2
N=2L+1=odd
anti-symmetric
h[k]=-h[N-k]
L
. a[k ].cos(.k )
k 0
LP/HP/BP
e
 jN / 2

L
cos( ). a[k ].cos( .k ) e
2 k 0
zero at   
LP/BP
Type-2
N=2L=even
symmetric
h[k]=h[N-k]
 jN / 2
L 1
Type-4
N=2L+1=odd
anti-symmetric
h[k]=-h[N-k]
sin( ). a[k ].cos(.k ) j.e
k 0
zero at   0, 
 jN / 2

L
sin( ). a[k ].cos(.k )
2 k 0
zero at   0
HP
FIR Filters 14/14
Linear-phase FIR filters
u[k]
• efficient direct-form realization.
example:
+
bo
y[k]








+
+
+
+
b1 b2 b3 b4
x
x
x
x
x
+
+
+
+
Filter Design by Optimization
(I) Weighted Least Squares Design :
• select one of the basic forms that yield linear phase
L
e.g. Type-1
j
 jN / 2
 jN / 2
H (e )  e
. a[k ].cos(.k )  e
. A( )
k 0
• specify desired frequency response (LP,HP,BP,…)
H d ( )  e jN / 2 .Ad ( )
e.g.: LP
Ad ( )  1,    P
( P  passband edge)
Ad ( )  0,  S    
( S  stopband edge)
• optimization criterion is


F (a0 ,...,aL )   W ( ) H (e j )  H d ( ) d   W ( ) A( )  Ad ( ) d

2
2

where W ( )  0 is a weighting function
Filter Design by Optimization
• …this is equivalent to F (a0 ,...,aL )  xT .Q.x  2 xT .b  
xT  a0
a1 ... aL 

Q   W ( ).c( ).cT ( )d
0

b   W ( ).Ad ( ).c( )d
0
cT ( )  1 cos( ) ... cos(L )
  constant
i.e. convex `Quadratic Optimization’ problem
• this is often supplemented with additional constraints...
Filter Design by Optimization
Example: Low-pass (LP) design
Ad ( )  1,    P
(passband edge)
Ad ( )  0,  S    
(stopband edge)
optimization criterion is
F (a0 ,...,aL ) 
P

0

A( )  1 d   .  A2 ( )d  xT .Q.x  2 xT .b  
2
S
an additional constraint may be imposed to control the
pass-band ripple …
A() 1   P ,   P ( P is passband ripple)
… as well as the stop-band ripple
A( )   S ,S     ( S is stopband ripple)
PS: Filter Specification
1.2
1  P
1
Passband Ripple
0.8
 P S
0.6
Passband Cutoff -> <- Stopband Cutoff
0.4
Stopband Ripple
S
0.2
0
0
0.5
1
1.5
2
2.5
3
Filter Design using `Windows’
Example : Low-pass filter design
• ideal low-pass filter is
1
  C
H d ( )  
0  C    
• hence ideal time-domain impulse response is
1
hd [k ] 
2

sin(c k )
 H d (e ).e d  ...   . c k
j
jk
• truncate hd[k] to N+1 samples :
hd [k ]
h[k ]  
 0
N /2k  N /2
ot herwise
• add (group) delay to turn into causal filter

Filter Design using `Windows’
Example : Low-pass filter design (continued)
• it can be shown that this corresponds to the solution of weighted leastsquares optimization with the given Hd, and weighting function W ( )for
1
all freqs.
• truncation corresponds to applying a `rectangular window’ :
h[k ]  hd [k ].w[k ]
1
w[k ]  
0
N /2 k  N /2
ot herwise
• +++: simple procedure (also for HP,BP,…)
• - - - : truncation in the time-domain results in `Gibbs effect’ in the
frequency domain, i.e. large ripple in pass-band and stop-band, which
cannot be reduced by increasing the filter order N.
Filter Design using `Windows’
Remedy : apply windows other than rectangular window:
• time-domain multiplication with a window function w[k] corresponds to
frequency domain convolution with W(z) :
h[k ]  hd [k ].w[k ]
H ( z)  H d ( z) *W ( z)
• candidate windows : Han, Hamming, Blackman, Kaiser,…. (see
textbooks, see DSP-I)
• window choice/design = trade-off between side-lobe levels (define peak
pass-/stop-band ripple) and width main-lobe (defines transition
bandwidth)
Design Procedure
•
To fully design and implement a filter five
steps are required:
(1)
(2)
(3)
(4)
(5)
Filter specification.
Coefficient calculation.
Structure selection.
Simulation (optional).
Implementation.
Filter Specification - Step 1
| H(f)|
pass-band
stop-band
1
fc : cut-off frequency
f(norm)
fs/2
(a)
| H(f)|
(dB)
pass-band
transition band
| H(f)|
(linear)
stop-band
p
1  
p
1
0
1  
p
pass-band
ripple
-3
stop-band
ripple
s
fsb : stop-band frequency
fc : cut-off frequency
fpb : pass-band frequency
(b)

fs/2
s
f(norm)
Coefficient Calculation - Step 2
•
There are several different methods
available, the most popular are:
– Window method.
– Frequency sampling.
– Parks-McClellan.
•
We will just consider the window method.
Window Method
•
•
First stage of this method is to calculate
the coefficients of the ideal filter.
This is calculated as follows:
1
(
)

hd n
2
1

2

H ( )e j n d 



c
j n

1
e
 d

c
2 fc sin (n c )


n c

2 fc
for n  0
for n  0
Window Method
•
Second stage of this method is to select a window function
based on the passband or attenuation specifications, then
determine the filter length based on the required width of the
transition band.
Window Type
Rectangular
Hanning
Normalised Transition
Width (f(Hz))
Passband Ripple(dB)
Stopband Attenuation
(dB)
0.7416
21
0.0546
44
0.9
N
3 .1
N
Hamming
3.3
N
0.0194
53
Blackman
5.5
N
0.0017
74
0.0274
50
0.000275
90
Kaiser
2.93
   4.54
N
5.71
   8.96
N
Using the Hamming
Window:
N
3. 3
3 .3

 8kHz  132
(1.2  1.4)kHz
f
Window Method
•
The third stage is to calculate the set of
truncated or windowed impulse response
coefficients, h[n]:
h(n )  hd (n )  W (n )
Where:
for

N 1
N 1
n
for N  odd
2
2
N
N
for N  even
 n
2
2
 2n 
W (n )  0.54  0.46 cos

 N 
 2n 
 0.54  0.46 cos

 133 
for
66  n  66
Window Method
•
Matlab code for calculating coefficients:
close all;
clear all;
fc = 8000/44100;
N = 133;
n = -((N-1)/2):((N-1)/2);
n = n+(n==0)*eps;
% cut-off frequency
% number of taps
[h] = sin(n*2*pi*fc)./(n*pi);
[w] = 0.54 + 0.46*cos(2*pi*n/N);
d = h.*w;
% generate sequence of ideal coefficients
% generate window function
% window the ideal coefficients
[g,f] = freqz(d,1,512,44100);
% transform into frequency domain for plotting
figure(1)
plot(f,20*log10(abs(g)));
axis([0 2*10^4 -70 10]);
% plot transfer function
figure(2);
stem(d);
xlabel('Coefficient number');
ylabel ('Value');
title('Truncated Impulse Response');
figure(3)
freqz(d,1,512,44100);
axis([0 2*10^4 -70 10]);
% avoiding division by zero
% plot coefficient values
% use freqz to plot magnitude and phase response
Truncated Impulse Response
Window Method
0.4
0.3
Value
0.2
0.1
0
20
40
60
80
Coefficient number
100
120
140
Magnitude (dB)
0
0
-20
-40
-60
0
0.2
0.4
0.6
0.8
1
1.2
Frequency (Hz)
1.4
1.6
1.8
2
4
x 10
0
Phase (degrees)
-0.1
-2000
-4000
-6000
0
0.5
1
Frequency (Hz)
1.5
2
4
x 10
Equiripple Design
• Starting point is minimax criterion, e.g.
F (a0 ,...,aL )  mina0 ,...,aL max0  W (). A()  Ad ()  mina0 ,...,aL max0  E()
• Based on theory of Chebyshev approximation and the `alternation
theorem’, which (roughly) states that the optimal ai’s are such that the
`max’ (maximum weighted approximation error) is obtained at L+2
extremal frequencies…
max0  E()  E(i ) for i  1,..,L  2
…that hence will exhibit the same maximum ripple (`equiripple’)
• Iterative procedure for computing extremal frequencies, etc. (Remez
exchange algorithm, Parks-McClellan algorithm)
• Very flexible, etc., available in many software packages
• Details omitted here (see textbooks)
Software
• FIR Filter design abundantly available in
commercial software
• Matlab:
b=fir1(n,Wn,type,window), windowed linear-phase FIR
design, n is filter order, Wn defines band-edges, type is
`high’,`stop’,…
b=fir2(n,f,m,window), windowed FIR design based on
inverse fourier transform with frequency points f and
corresponding magnitude response m
b=remez(n,f,m), equiripple linear-phase FIR design with
Parks-McClellan (Remez exchange) algorithm