Title: Digital Signal Processing FIR Filter Design
1Digital Signal Processing FIR Filter Design
- Marc Moonen
- Dept. E.E./ESAT, K.U.Leuven
2FIR 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,)
3Review of discrete-time systems 1/10
- Linear time-invariant (LTI) systems
- Linear
- input u1k -gt output y1k
- input u2k -gt output y2k
- hence input a.u1kb.u2k-gt
a.y1kb.y2k - Time-invariant (shift-invariant)
- input uk -gt output yk
- hence input uk-T -gt output yk-T
uk
yk
4Review of discrete-time systems 2/10
- Linear time-invariant (LTI) systems
- Causal systems
- for all input uk0, klt0 -gt output
yk0, klt0 - Impulse response
- input 1,0,0,0,... -gt output
h0,h1,h2,h3,... - input u0,u1,u2,u3 -gt output
y0,y1,y2,y3,... -
-
- convolution
5Review of discrete-time systems 3/10
- Impulse response/convolution
-
-
-
-
-
Toeplitz matrix
6Review of discrete-time systems 4/10
7Review of discrete-time systems 5/10
- Z-Transform
- input-output relation
- shorthand notation
- (for convolution operation/Toeplitz-vector
product) - stability
- bounded input uk -gt bounded output yk
- --iff
-
- --iff poles of H(z) inside the unit circle
- (for causal,rational systems)
-
-
8Review of discrete-time systems 6/10
- Frequency response
- given a system with impulse response hk
- given an input signal complex sinusoid
- output signal
-
-
- is frequency
response - is H(z) evaluated on
the unit circle
9Review of discrete-time systems 7/10
- Frequency response
- for a real impulse response hk
- Magnitude response is
even function - Phase response
is odd function - example
Nyquist frequency
10Review 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
11Review of discrete-time systems 9/10
- Rational transfer functions (IIR filters)
-
-
- N poles (zeros of A(z)) , N zeros (zeros of B(z))
- corresponds to difference equation
12Review of discrete-time systems 10/10
- FIR filters (finite impulse response)
-
-
- Moving average filters (MA)
- N poles at the origin z0 (hence guaranteed
stability) - N zeros (zeros of B(z)), all zero filters
- corresponds to difference equation
- impulse response
13Review 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
14FIR Filters 1/14
15FIR Filters 2/14
- retiming select subgraph (shaded)
- remove delay element
on all inbound arrows - add delay element on
all outbound arrows
uk
uk-4
uk-3
uk-2
uk-1
b4
b3
b2
b1
x
x
x
x
yk
16FIR Filters 3/14
uk
uk-1
uk-3
uk-2
b1
b4
b3
b2
x
x
x
x
yk
17FIR Filters 4/14
- retiming repeated application results in...
- Transposed direct-form realization
18FIR Filters 5/14
- Lattice form derived from combined
realization with - reversed coefficient vector results in
- - same magnitude response
- - different phase response
19FIR Filters 6/14
- Lattice form derivation (I)
uk
uk-1
uk-2
uk-3
uk-4
b1
b2
b3
bo
b4
x
x
x
x
x
b3
b2
b4
b1
bo
x
x
x
x
x
yk
20FIR Filters 7/14
- Lattice form derivation (II), this is
equivalent to... -
(simple proof)
uk
uk-1
uk-2
uk-3
uk-4
b1
b2
b3
bo
0
x
x
x
x
x
b3
b2
yk
0
b1
bo
x
x
x
x
ko
21FIR Filters 8/14
- Lattice form derivation (III), retiming leads
to... - repeat
procedure for shaded graph
uk
uk-3
uk-2
uk-2
bo
b1
b2
b3
x
x
x
x
b3
b2
yk
b1
bo
x
x
x
x
ko
22FIR Filters 9/14
- Lattice form derivation (IV), end result
is... -
uk
k4
yk
ko
k1
k2
k3
23FIR Filters 10/14
- Lattice form
- kis are reflection coefficients
- procedure for computing kis from bis
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
kilt1 (i1,,N-1) - (ps procedure breaks down if ki1 is
encountered) -
24FIR Filters 11/14
- Linear-phase FIR filters
- Non-causal zero-phase filters
- example symmetric impulse response
- h-L,.h-1,h0,h1,...,hL
- hkh-k, k1..L
-
- frequency response is
-
- - real-valued (zero-phase) transfer
function - - causal implementation by introducing
(group) delay
25FIR Filters 12/14
- Linear-phase FIR filters
- Causal linear-phase filters
- example symmetric impulse response N even
- h0,h1,.,hN
- N2L (even)
- hkhN-k, k0..L
- frequency response is
-
- - causal implementation of zero-phase
filter, by - introducing (group) delay
k
N
0
26FIR Filters 13/14
- Linear-phase FIR filters
- Type-1 Type-2 Type-2
Type-4 - N2Leven N2L1odd N2Leven
N2L1odd - symmetric anti-symmetric symmetric
anti-symmetric - hkhN-k hk-hN-k
hkhN-k hk-hN-k -
-
- zero at
zero at zero at - LP/HP/BP LP/BP
HP
27FIR Filters 14/14
- Linear-phase FIR filters
- efficient direct-form realization.
- example
bo
b4
b3
b2
b1
x
x
x
x
x
yk
28Filter Design by Optimization
- (I) Weighted Least Squares Design
- select one of the basic forms that yield linear
phase - e.g. Type-1
- specify desired frequency response (LP,HP,BP,)
- e.g. LP
-
- optimization criterion is
- where is a weighting function
29Filter Design by Optimization
- this is equivalent to
- i.e. convex Quadratic Optimization problem
- this is often supplemented with additional
constraints...
30Filter Design by Optimization
- Example Low-pass (LP) design
-
- optimization criterion is
- an additional constraint may be imposed to
control the pass-band ripple - as well as the stop-band ripple
31PS Filter Specification
32Filter Design using Windows
- Example Low-pass filter design
- ideal low-pass filter is
- hence ideal time-domain impulse response is
- truncate hdk to N1 samples
- add (group) delay to turn into causal filter
33Filter Design using Windows
- Example Low-pass filter design (continued)
- it can be shown that this corresponds to the
solution of weighted least-squares optimization
with the given Hd, and weighting function
for all freqs. - truncation corresponds to applying a rectangular
window - 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.
34Filter Design using Windows
- Remedy apply windows other than rectangular
window - time-domain multiplication with a window function
wk corresponds to frequency domain convolution
with 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)
35Design Procedure
- To fully design and implement a filter five steps
are required - (1) Filter specification.
- (2) Coefficient calculation.
- (3) Structure selection.
- (4) Simulation (optional).
- (5) Implementation.
36Filter Specification - Step 1
37Coefficient 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.
38Window Method
- First stage of this method is to calculate the
coefficients of the ideal filter. - This is calculated as follows
39Window 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.
Using the Hamming Window
40Window Method
- The third stage is to calculate the set of
truncated or windowed impulse response
coefficients, hn
for
Where
for
41Window Method
- Matlab code for calculating coefficients
close all clear all fc 8000/44100
cut-off frequency N 133 number of taps n
-((N-1)/2)((N-1)/2) n n(n0)eps
avoiding division by zero h
sin(n2pifc)./(npi) generate sequence of
ideal coefficients w 0.54
0.46cos(2pin/N) generate window function d
h.w 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))) plot transfer
function axis(0 2104 -70 10) figure(2) ste
m(d) plot coefficient values xlabel('Coeffic
ient number') ylabel ('Value') title('Truncated
Impulse Response') figure(3) freqz(d,1,512,4410
0) use freqz to plot magnitude and phase
response axis(0 2104 -70 10)
42Window Method
43Equiripple Design
- Starting point is minimax criterion, e.g.
- Based on theory of Chebyshev approximation and
the alternation theorem, which (roughly) states
that the optimal ais are such that the max
(maximum weighted approximation error) is
obtained at L2 extremal frequencies - 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)
44Software
- FIR Filter design abundantly available in
commercial software - Matlab
- bfir1(n,Wn,type,window), windowed linear-phase
FIR design, n is filter order, Wn defines
band-edges, type is high,stop, - bfir2(n,f,m,window), windowed FIR design
based on inverse fourier transform with frequency
points f and corresponding magnitude response m - bremez(n,f,m), equiripple linear-phase FIR
design with Parks-McClellan (Remez exchange)
algorithm