2.161: Lecture 17 Notes

Massachusetts Institute of Technology Department of Mechanical Engineering 2.161 Signal Processing - Continuous and Discrete Lecture 17 Reading: 1 • Class Handout: Frequency-Sampling Filters. • Proakis and Manolakis: Secs. 10.2.3, 10.2.4 • Oppenheim, Schafer and Buck: 7.4 • Cartinhour: Ch. 9 Frequency-Sampling Filters In the frequency-sampling filters the parameters that characterize the the filter are the values of the desired frequency response H(ejω ) at a discrete set of equally spaced sampling frequencies. In particular, let ωk = 2π k N k = 0, . . . , N − 1 (1) as shown below for the cases of N even, and N odd. Note that when N is odd there is no sample at the Nyquist frequency, ω = π. The frequency-sampling method guarantees that the resulting filter design will meet the given design specification at each of the sample frequencies. Im (z ) Im (z ) z - p la n e z - p la n e M  2 F 1 0 -1 1 M R e (z ) 2 F 1 1 -1 N = 1 0 (e v e n ) 17–1 1 R e (z ) N = 1 1 (o d d ) (b ) (a )  For convenience denote the complete sample set {Hk } as Hk = H(ejωk ) k = 1, . . . , N − 1. For a filter with a real impulse response {hn } we require conjugate symmetry, that is ¯k HN −k = H and further, for a filter with a real, even impulse response we require {Hk } to be real and even, that is HN −k = Hk . Within these constraints, it is sufficient to specify frequency samples for the upper half of the z-plane, that is for  2π k = 0, . . . , N 2−1 N odd k ωk = N N even. k = 0, . . . , 2 N and use the symmetry constraints to determine the other samples. If we assume that H(ejω ) may be recovered from the complete sample set {Hk } by the cardinal sinc interpolation method, that is H(ejω ) = N −1  Hk k=0 sin (ω − 2πk/N ) ω − 2πk/N then H(ejω ) is completely specified by its sample set, and the impulse response, of length N , may be found directly from the inverse DFT, {hn } = IDFT {Hk } where N −1 2πkn 1  hn = Hk ej N N k=0 n = 0, . . . , N − 1 As mentioned above, this method guarantees that the resulting FIR filter, represented by {hn }, will meet the specification H(ejω ) = Hk at ω = ωk = 2kπ/N . Between the given sampling frequencies the response H(ejω ) will be described by the cardinal interpolation. 1.1 Linear-Phase Frequency-Sampling Filter The filter described above is finite, with length N , but is non-causal. To create a causal filter with a linear phase characteristic we require an impulse response that is real and symmetric about its mid-point. This can be done by shifting the computed impulse response to the right by (N − 1)/2 samples to form H  (z) = z −(N −1)/2 H(z) but this involves a non-integer shift for even N . Instead, it is more convenient to add the appropriate phase taper to the frequency domain samples Hk before taking the IDFT. The non-integer delay then poses no problems: 17–2 • Apply a phase shift of πk(N − 1) N to each of the samples in the upper half z-plane  k = 0, . . . , (N − 1)/2  jφk Hk = Hk e k = 0, . . . , N/2 φk = − (2) (for n odd) (for n even) • Force the lower half plane samples to be complex conjugates.  k = 1, . . . , (N − 1)/2 (for n odd)   HN −k = H̄k k = 1, . . . , N/2 − 1 (for n even) • Then the linear-phase impulse response is {hn } = IDFT {Hk } 1.2 A Simple MATLAB Frequency-Sampling Filter The following is a MATLAB script of a tutorial frequency-sampling filter h = firfs(samples) that takes a vector samples of length N of the desired frequency response in the range 0 ≤ ω ≤ (N − 1)π/N , and returns the linear-phase impulse response {hn } of length 2N − 1. ------------------------------------------------------------------------% 2.161 Classroom Example - firfs - A simple Frequency-Sampling Linear-Phase FIR % Filter based on DFT interpolation. % Usage : h = firfs(samples) % where samples - is a row vector of M equi-spaced, real values % of the freq. response magnitude. % The samples are interpreted as being equally spaced around % the top half of the unit circle at normalized (in terms of % the Nyquist frequency f_N) frequencies from % 0 to 2(M-1)/(2M-1) x f_N, % or at frequencies 2k/(2N-1)xf_N for k = 0...M-1 % Note: Because the length is odd, the frequency response % is not specified at f_N. % h - is the output impulse response of length 2M-1 (odd). % The filter h is real, and has linear phase, i.e. has symmetric % coefficients obeying h(k) = h(2M+1-k), k = 1,2,...,M+1. %------------------------------------------------------------------------­ function h = firfs(samples) % % Find the length of the input array... % The complete sample set on the unit circle will be of length (2N-1) 17–3 % N = 2*length(samples) -1; H_d = zeros(1,N); % % We want a causal filter, so the resulting impulse response will be shifted % (N-1)/2 to the right. % Move the samples into the upper and lower halves of H_d and add the % linear phase shift term to each sample. % Phi = pi*(N-1)/N; H_d(1) = samples(1); for j = 2:N/2-1 Phase = exp(-i*(j-1)*Phi); H_d(j) = samples(j)*Phase; H_d(N+2-j) = samples(j)*conj(Phase); end % % Use the inverse DFT to define the impulse response. % h = real(ifft(H_d)); The following MATLAB commands were used to generate a filter with 22 frequency samples, generating a length 43 filter. h=firfs([1 1 1 1 0.4 0 0 0 0 0.8 2 2 2 2 0.8 0 0 0 0 0 0 0 ]); freqz(h,1) The filter has two pass-bands; a low-pass region with a gain of unity, and a band-pass region with a gain of two. Notice that the band-edges have been specified with transition samples, this is discussed further below. The above commands produced the following frequency response for the filter. Magnitude (dB) 0 −20 −40 −60 −80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Normalized Frequency ( π rad/sample) 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Normalized Frequency ( π rad/sample) 0.8 0.9 1 Phase (degrees) 0 −500 −1000 −1500 −2000 17–4 1.3 The Effect of Band-Edge Transition Samples One of the advantages of the frequency-sampling filter is that the band-edges may be more precisely specified than the window method. For example, low-pass filters might be specified by h = firfs([1 1 1 1 1 0.4 0 0 0 0 0 0]); with one transition value of 0.4, or h = firfs([1 1 1 1 0.7 0.2 0 0 0 0 0 0]); with a pair of transition specifications. The frequency-sampling filter characteristic will pass through these points, and they can have a significant effect on the stop-band characteristics of the filter. The figure below shows the effect of varying the value of a single transition point in a filter of length N = 33. 10 t = 0.6 t = 0.4 t = 0.2 0 Magnitude (dB) −10 −20 −30 −40 −50 −60 −70 −80 0 0.1 0.2 0.3 0.4 0.5 0.6 Frequency f/FN 0.7 0.8 0.9 1 The values shown are for t = 0.6, 0.4 and 0.2. There is clearly a significant improvement in the stop-band attenuation for for the case t = 0.4. Similarly the following figure compares the best of these single transition values (t = 0.4) with a the response using two transition points (t1 = 0.59, t2 = 0.11). The filter using two transition points shows a significant improvement in the stop-band over the single point case, at the expense of the transition width. 10 t1 = .59, t2 = 0.11 t = 0.4 0 −10 Magnitude (dB) −20 −30 −40 −50 −60 −70 −80 0 0.1 0.2 0.3 0.4 0.5 0.6 Frequency f/fN 17–5 0.7 0.8 0.9 1 Rabiner et al. (1970) did an extensive linear-programming optimization study to deter­ mine the optimum value of band edge transition values, and tabulated the results for even and odd filters of different lengths. The results show that for one transition point topt ≈ 0.4, and for two points topt ≈ 0.59, and 0.11. 2 FIR Filter Design Using Optimization These methods allow much greater flexibility in the filter specification. In general they seek the filter coefficients that minimize the error (in some sense) between a desired frequency response Hd ( ej ω ) and the achieved frequency response H( ej ω ). The most common optimiza­ tion method is that due to Parks and McClellan (1972) and is widely available in software filter design packages (including MATLAB). The Parks-McClellan method allows • Multiple pass- and stop-bands. • Is an equi-ripple design in the pass- and stop-bands, but allows independent weighting of the ripple in each band. • Allows specification of the band edges. tr a n s itio n e q u i- r ip p le 1 + @ 1 1 1 -@ 1 p a s s -b a n d s to p -b a n d e q u i- r ip p le @ 0 2 M M c F s M For the low-pass filter shown above the specification would be 1 − δ1 < H( ej ω ) < 1 + δ1 −δ2 < H( ej ω ) < δ2 in the pass-band 0 < ω ≤ ωc in the stop-band ωs < ω ≤ π. where the ripple amplitudes δ1 and δ2 need not be equal. Given these specifications we need to determine, the length of the filter M + 1 and the filter coefficients {hn } that meet the specifications in some optimal sense. If M + 1 is odd, and we assume even symmetry hM −k = hk k = 0 . . . M/2 17–6 and the frequency response function can be written jω H( e ) = h0 + 2 M/2  hk cos(ωk) k=1  M/2 = ak cos(ωk) k=0 Let Hd ( ej ω ) be the desired frequency response, and define a weighted error   E( ej ω ) = W ( ej ω ) Hd ( ej ω ) − H( ej ω ) where W ( ej ω ) is a frequency dependent weighting function, but by convention let W ( ej ω ) be constant across each of the critical bands, and zero in all transition bands. In particular for the low-pass design ⎧ ⎪δ2 /δ1 in the pass-band ⎨ jω W(e ) = 1 in the stop-band ⎪ ⎩ 0 in the transition band. This states that the optimization will control the ratio of the pass-band to stop-band ripple, and that the transition will not contribute to the error criterion. Let Ω be a compact subset of the frequency band from 0 to π representing the pass- and stop-bands. The goal is to find the set of filter parameters {ak }, k = 0, . . . , M/2 + 1 that minimize the maximum value of the error E( ej ω ) over the pass- and stop-bands. ⎞ ⎤⎤ ⎡ ⎡ ⎛ M/2   min max E( ej ω ) = min ⎣ max ⎣ W ( ej ω ) ⎝Hd ( ej ω ) − ak cos(ωk)⎠ ⎦⎦ over ak over Ω over ak over Ω k=0 where Ω is the disjoint set of frequency bands that make up the pass- and stop-bands of the filter. The solution is found by an iterative optimization routine. We do not attempt to cover the details of the algorithm here, and merely note: • The method is based on reformulating the problem as one in polynomial approximation, using Chebyshev polynomials, where cos(ωk) = Tk (cos(ω)) where Tk (x) is a polynomial of degree k, (see the Class Handout on Chebyshev filter design). Consequently jω H( e ) = M/2  ak cos(ωk) = k=0 M/2  k=0 17–7 ak (cos(ω))k • The algorithm uses Chebyshev’s alternation theorem to recognize the optimal solution. In general terms the theorem is stated: Define the error E(x) as above, namely   E( ej ω ) = W ( ej ω ) Hd ( ej ω ) − H( ej ω ) and the maximum error as E( ej ω )∞ = argmaxx∈Ω E( ej ω ) A necessary and sufficient condition that H( ej ω ) is the unique Lth-order polynomial minimizing E( ej ω )∞ is that E( ej ω ) exhibit at least L + 2 extremal frequencies, or “alternations”, that is there must exist at least L+2 values of ω, ωk ∈ Ω, k = [0, 1, . . . , L + 1], such that ω0 < ω1 < . . . < ωL+1 , and such that   E( ej ωk ) = −E( ej ωk+1 ) = ± E( ej ω )∞ . Note that the alternation theorem is simply a way of recognizing the optimal equi­ ripple solution. For example, the following figure is from a Parks-McClellan low-pass filter with length M + 1 = 17. | H ( e jM ) | 1 + @ tr a n s itio n 1 1 1 -@ 1 p a s s -b a n d s to p -b a n d a lte r n a tio n fr e q u e n c ie s @ 2 0 M c M F s M From above, H( ej ω ) is written as a polynomial of degree M/2, jω H( e ) = M/2  ak (cos(ω))k k=0 so that L = M/2 and the pass- and stop-bands must exhibit at least M/2 + 2 = 10 points of alternation. These 10 points are shown in the figure. • The Parks-McClellan algoritm uses the Remez exchange optimization method. 17–8 See Proakis and Manolakis Sec. 10.2.4 or Openheim, Schafer and Buck Sec. 7.4 for details. MATLAB Parks-McClellan Function: The Parks-McClellan algorithm is implemented in the MATLAB function b = firpm(M,F,A,W) where b is the array of filter coefficients, M is the filter order (M+1 is the length of the filter), F is a vector of band edge frequencies in ascending order, A is a set of filter gains at the band edges, and W is an optional set of relative weights to be applied to each of the bands. For example, consider a band-pass filter with two pass-bands as shown below: H d (e ) jM p a s s 1 p a s s 0 .7 0 s to p M s to p 1 tr a n s itio n M M 2 3 tr a n s itio n M s to p M 4 5 tr a n s itio n M M 6 7 tr a n s itio n M F 8 M There are five distinct bands in this filter, separated by four transition regions. The filter would require the following specifications: F = [0 ω1 ω2 ω3 ω4 ω5 ω6 ω7 ω8 1] A = [0 0 1 1 0 0 0.7 0.7 0 0] W = [10 1 10 1 10] where the errors in the stop-bands have been weighted 10 times more heavily than in the pass-bands. See the MATLAB help/documentation for more details. Example 1 Design a length 33 Parks-McClellan band-pass filter with the following band specifications: H d (e jM ) p a s s 1 0 0 s to p 0 .2 F tr a n s itio n s to p 0 .4 F 0 .7 F 17–9 tr a n s itio n 0 .8 5 F F M Weight the stop-band ripple ten times more heavily than the pass-band. Solution: h=firpm(32,[0 0.2 0.4 0.7 0.85 1],[0 0 10 10 0 0],[10 1 10]) freqz(h,1) M a g n itu d e ( d B ) 5 0 0 -5 0 0 -1 0 0 P h a s e (d e g re e s ) 1 0 0 0 0 .2 0 .4 N o r m a liz e d F r e q u e n c y (x p 0 .2 0 .4 N o r m a liz e d F r e q u e n c y (x p 0 .6 r a d /s a m p le ) 0 .8 1 0 .8 1 5 0 0 0 -5 0 0 -1 0 0 0 -1 5 0 0 0 17–10 0 .6 r a d /s a m p le )