330 likes | 581 Views
ELEN E4810: Digital Signal Processing Topic 9: Filter Design: FIR. Windowed Impulse Response Window Shapes Design by Iterative Optimization. 1. FIR Filter Design. FIR filters no poles (just zeros) no precedent in analog filter design Approaches windowing ideal impulse response
E N D
ELEN E4810: Digital Signal ProcessingTopic 9: Filter Design: FIR • Windowed Impulse Response • Window Shapes • Design by Iterative Optimization Dan Ellis
1. FIR Filter Design • FIR filters • no poles (just zeros) • no precedent in analog filter design • Approaches • windowing ideal impulse response • iterative (computer-aided) design Dan Ellis
Least Integral-Squared Error • Given desired FR Hd(ejw), what is the best finite ht[n] to approximate it? • Can try to minimize Integral Squared Error (ISE) of frequency responses: best in what sense? = DTFT{ht[n]} Dan Ellis
Least Integral-Squared Error • Ideal IR is hd[n] = IDTFT{Hd(ejw)}, (usually infinite-extent) • By Parseval, ISE • But: ht[n] only exists for n = -M..M , minimized by making ht[n] = hd[n], -M ≤ n ≤ M not altered byht[n] Dan Ellis
Least Integral-Squared Error • Thus, minimum mean-squared error approximation in 2M+1 point FIR is truncated IDFT: • Make causal by delaying by M points h't[n] = 0 forn < 0 Dan Ellis
w -p -wc wc p Approximating Ideal Filters • From topic 6,ideal lowpasshas:and: (doubly infinite) Dan Ellis
Approximating Ideal Filters • Thus, minimum ISEcausal approximation to an ideal lowpass Causal shift 2M+1 points Dan Ellis
Freq. Resp. (FR) Arithmetic • Ideal LPF has pure-real FR i.e.q(w) = 0, H(ejw) = |H(ejw)| Can build piecewise-constant FRs by combining ideal responses, e.g. HPF: d[n] i.e. H(ejw) = 1 1 w – HLP(ejw) = 1 for |w| < wc wouldn’t work if phases were nonzero! hLP[n] 1 w = hHP[n] 1 = d[n] - (sinwcn)/pn w -p -wc wc p Dan Ellis
Gibbs Phenomenon • Truncated ideal filters have Gibbs’ Ears: Increasing filter lengthnarrower ears(reduces ISE)but height the same (11% overshoot) not optimal by minimax criterion Dan Ellis
Where Gibbs comes from • Truncation of hd[n] to 2M+1 points is multiplication by a rectangular window: • Multiplication in time domain is convolution in frequency domain: wR[n] Dan Ellis
Where Gibbs comes from • Thus, FR of truncated response is convolution of ideal FR and FR of rectangual window (pdc.sinc): Hd(ejw) Ht(ejw) DTFT{wR[n]}periodic sinc... Dan Ellis
Where Gibbs comes from • Rectangular window: • Mainlobe width (1/L) determinestransition band • Sidelobe heightdetermines ripples doesn’t vary with length “periodic sinc” Dan Ellis
2. Window Shapes for Filters • Windowing (infinite) ideal response FIR filter: • Rectangular window has best ISE error • Other “tapered windows” vary in: • mainlobe transition band width • sidelobes size of ripples near transition • Variety of ‘classic’ windows... Dan Ellis
Window Shapes for FIR Filters bigsidelobes! • Rectangular: • Hann: • Hamming: • Blackman: double widthmainlobe reduced 1stsidelobe triple width mainlobe Dan Ellis
2p 2M+1 Window Shapes for FIR Filters • Comparison on dB scale: Dan Ellis
Adjustable Windows • Have discrete main-sidelobe tradeoffs... • Kaiser window = parametric, continuous tradeoff: • Empirically, for min. SB atten. of a dB: modified zero-order Bessel function requiredorder transitionwidth Dan Ellis
H(ejw) w p4 kHz 2p 8 kHz 0.15 p600 Hz Windowed Filter Example • Design a 25 point FIR low-pass filter with a cutoff of 600 Hz (SR = 8 kHz) • No specific transition/ripple req’s compromise: use Hamming window • Convert the frequency to radians/sample: Dan Ellis
Windowed Filter Example • Get ideal filter impulse response: • Get window: Hamming @ N = 25 M = 12 (N = 2M+1) • Apply window: M Dan Ellis
3. Iterative FIR Filter Design • Can derive filter coefficients by iterative optimization: • Gradient descent / nonlinear optimiz’n desiredresponse H(ejw) Goodness of fit criterion error e Filter coefs h[n] Update filterto reduce e Estimate derivatives e/h[n] Dan Ellis
errormeasurement region exponent:2 least sq’s minimax actualresponse error weighting desiredresponse Error Criteria = W(w)·[D(ejw) – H(ejw)] Dan Ellis
Minimax FIR Filters • Iterative design of FIR filters with: • equiripple (minimax criterion) • linear-phase symmetric IR h[n] = (–)h[-n] • Recall, symmetric FIR filters have FR with pure-real i.e. combo of cosines of multiples of w (type I) n M a[0] = h[M]a[k] = 2h[M - k] Dan Ellis
Minimax FIR Filters • Now, cos(kw) can be expressed as a polynomial in cos(w)k and lower powers • e.g. cos(2w) = 2(cosw)2 - 1 • Thus, we can find a’s such that • Mth order polynomial in cosw • a[k]s are simply related to a[k]s Mth order polynomial in cosw Dan Ellis
Minimax FIR Filters Mth order polynomial in cosw • An Mth order polynomial has at most M - 1 maxima and minima: has at most M-1 min/max (ripples) 1 5th order polynomial in cosw 3 2 4 Dan Ellis
Alternation Theorum • Key ingredient to Parks-McClellan: is the unique, best, weighted-minimax order 2M approx. to D(ejw) • has at least M+2 “extremal” freqs over w subset R • error magnitude is equal at each extremal: • peak error alternates in sign: Dan Ellis
Alternation Theorum • Hence, for a frequency response: • Ife(w) reaches a peak error magnitude e at some set of extremal frequenceswi • And the sign of the peak error alternates • And we have at least M+2 of them • Thenoptimal minimax (10th order filter, M = 5) Dan Ellis
Alternation Theorum • By Alternation Theorum, M+2extrema of alternating signs optimal minimax filter • But has at most M-1 extrema need at least 3 more from band edges • 2 bands give 4 band edges can afford to “miss” only one • Alternation rules out transition band edges, thus have 1 or 2 outer edges Dan Ellis
Alternation Theorum • For M = 5 (10th order): • 8 extrema (M+3, 4 band edges) - great! • 7 extrema (M+2, 3 band edges) - OK! • 6 extrema (M+1, only 2 transition band edges)NOT OPTIMAL Dan Ellis
Parks-McClellan Algorithm • To recap: • FIR CAD constraintsD(ejw), W(w) e(w) • Zero-phase FIRH(w) = SkakcoskwM-1 min/max • Alternation theorumoptimal ≥M+2 pk errs, alter’ng sign • Hence, can spot ‘best’ filter when we see it – but how to find it? ~ Dan Ellis
Parks-McClellan Algorithm ~ ~ • Alternation [H(w)-D(w)]/W(w) must = ±e at M+2 (unknown) frequencies {wi}... • Iteratively update h[n] with Remez exchange algorithm: • estimate/guess M+2 extremals {wi} • solve for a[n], e ( h[n] ) • find actual min/max in e(w) new {wi} • repeat until |e(wi)|is constant • Converges rapidly! Dan Ellis
Parks-McClellan Algorithm • In Matlab, >> h=remez(10, [0 0.4 0.6 1], [1 1 0 0], [1 2]); filter order (2M) band edges ÷ p desired magnitudeat band edges error weightsper band Dan Ellis