Digital Signal Processing FIR Filter Design

1. Digital Signal Processing FIR Filter Design Marc Moonen Dept. E.E./ESAT, K.U.Leuven

2. 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,…)

3. Review of discrete-time systems 1/10 Linear time-invariant (LTI) systems • 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] u[k] y[k]

4. u[k] y[k] Review of discrete-time systems 2/10 Linear time-invariant (LTI) systems • 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],... = `convolution’

5. u[0],u[1],u[2],u[3],0,0…. y[0],y[1],... h[0],h[1],h[2],0,0,... Review of discrete-time systems 3/10 • Impulse response/convolution: `Toeplitz’ matrix

6. Review of discrete-time systems 4/10 • Z-Transform:

7. Review of discrete-time systems 5/10 Z-Transform : • input-output relation • `shorthand’ notation (for convolution operation/Toeplitz-vector product) • stability bounded input u[k] -> bounded output y[k] --iff --iff poles of H(z) inside the unit circle (for causal,rational systems)

8. 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 : is `frequency response’ is H(z) evaluated on the unit circle

9. Review of discrete-time systems 7/10 Frequency response : • for a real impulse response h[k] Magnitude response is even function Phase response is odd function • example : Nyquist frequency

10. 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

11. Review 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

12. Review of discrete-time systems 10/10 `FIR filters’ (finite impulse response): • `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 • impulse response

13. 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

14. u[k] u[k-1] u[k-2] u[k-3] u[k-4] bo b1 b2 b3 b4 x x x x x y[k] + + + + FIR Filters 1/14 • `direct-form’ realization

15. bo 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] u[k-4] b1 b2 b3 b4 x x x x y[k] + + +

16. bo x + FIR Filters 3/14 • `retiming’ : results in... u[k] u[k-1] u[k-2] u[k-3] b1 b2 b3 b4 x x x x y[k] + + +

17. u[k] bo b1 b2 b3 b4 x x x x x y[k] + + + + FIR Filters 4/14 • `retiming’ :repeated application results in... `Transposed direct-form’ realization

18. FIR Filters 5/14 • `Lattice form’ : derived from combined realization with reversed coefficient vector results in - same magnitude response - different phase response

19. y[k] FIR Filters 6/14 • `Lattice form’ : derivation (I) u[k] u[k-1] u[k-2] u[k-3] u[k-4] bo b1 b2 b3 b4 x x x x x b4 b3 b2 b1 bo x x x x x y[k] + + + + + + + +

20. x x x + + + y[k] FIR Filters 7/14 • `Lattice form’ : derivation (II), this is equivalent to... (=simple proof) u[k] u[k-1] u[k-2] u[k-3] u[k-4] b’o b’1 b’2 b’3 0 x x x x x y[k] 0 b’3 b’2 b’1 b’o x x x x + + + ko + + + +

21. x x + + y[k] FIR Filters 8/14 • `Lattice form’ : derivation (III), retiming leads to... …repeat procedure for shaded graph u[k] u[k-2] u[k-2] u[k-3] b’o b’1 b’2 b’3 x x x x y[k] b’3 b’2 b’1 b’o x x x x + + + ko + + +

22. x x x x x x x x x + + + + + + + + y[k] FIR Filters 9/14 • `Lattice form’ : derivation (IV), end result is... u[k] k4 y[k] ko k1 k2 k3

23. 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)

24. k L 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 frequency response is - real-valued (=zero-phase) transfer function - causal implementation by introducing (group) delay

25. 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 frequency response is - causal implementation of zero-phase filter, by introducing (group) delay k 0 N

26. FIR Filters 13/14 Linear-phase FIR filters Type-1 Type-2 Type-2 Type-4 N=2L=even N=2L+1=odd N=2L=even N=2L+1=odd symmetric anti-symmetric symmetric anti-symmetric h[k]=h[N-k] h[k]=-h[N-k] h[k]=h[N-k] h[k]=-h[N-k] zero at zero at zero at LP/HP/BP LP/BP HP

27. u[k] + + + + + FIR Filters 14/14 Linear-phase FIR filters • efficient direct-form realization. example: bo b1 b2 b3 b4 x x x x x y[k] + + + +

28. Filter 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

29. Filter Design by Optimization • …this is equivalent to i.e. convex `Quadratic Optimization’ problem • this is often supplemented with additional constraints...

30. Filter 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

31. PS: Filter Specification

32. Filter Design using `Windows’ Example : Low-pass filter design • ideal low-pass filter is • hence ideal time-domain impulse response is • truncate hd[k] to N+1 samples : • add (group) delay to turn into causal filter

33. Filter 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.

34. 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) : • 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)

35. Design 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.

36. Filter Specification - Step 1

37. 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.

38. p 1 ( ) ( ) ò w j n = w w h n H e d d p 2 - p w c 1 ò w j n = × w 1 e d p 2 - w c ( ) w ì 2 f sin n c c ï ¹ for n 0 = w í n c ï = 2 f for n 0 î c Window Method • First stage of this method is to calculate the coefficients of the ideal filter. • This is calculated as follows:

39. 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. Using the Hamming Window:

40. Window Method • The third stage is to calculate the set of truncated or windowed impulse response coefficients, h[n]: for Where: for

41. Window 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+(n==0)*eps; % avoiding division by zero [h] = sin(n*2*pi*fc)./(n*pi); % generate sequence of ideal coefficients [w] = 0.54 + 0.46*cos(2*pi*n/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 2*10^4 -70 10]); figure(2); stem(d); % plot coefficient values xlabel('Coefficient number'); ylabel ('Value'); title('Truncated Impulse Response'); figure(3) freqz(d,1,512,44100); % use freqz to plot magnitude and phase response axis([0 2*10^4 -70 10]);

42. Window Method

43. Equiripple Design • Starting point is minimax criterion, e.g. • 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… …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)

44. 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