Skip to content Skip to navigation

Connexions

You are here: Home » Content » Linear-Phase Fir Filter Design By Least Squares

Navigation

Content Actions
  • Download module PDF
  • Add to ...
    Add the module to:
    • My Favorites
    • A lens
    • An external social bookmarking service
    • My Favorites (What is 'My Favorites'?)
      'My Favorites' is a special kind of lens which you can use to bookmark modules and collections directly in Connexions. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need a Connexions account to use 'My Favorites'.
    • A lens (What is a lens?)

      Definition of a lens

      Lenses

      A lens is a custom view of Connexions content. You can think of it as a fancy kind of list that will let you see Connexions through the eyes of organizations and people you trust.

      What is in a lens?

      Lens makers point to Connexions materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

      Who can create a lens?

      Any individual Connexions member, a community, or a respected organization.

    • External bookmarks
  • E-mail the author

Recently Viewed

Linear-Phase Fir Filter Design By Least Squares

Module by: Ivan Selesnick

Summary: This module describes the design of linear-phase FIR filters based on the square error criterion. It includes derivation and an example.

Linear-Phase FIR Filter Design by Least Squares

This module describes the design of linear-phase FIR filters based on the square error criterion. We will see that FIR filters that minimize the square error can be found by solving a linear system of equations. This technique is straight-forward and is applicable to arbitrary desired frequency responses. In addition, linear constraints on the coefficients are easily included.

Linear-phase filter design by least squares has several advantages.

  • Optimal with respect to square error criterion.
  • Simple, non-iterative method.
  • Analytic solutions sometimes possible, otherwise solution is obtained via solution to linear system of equations.
  • Allows the use of a frequency dependent weighting function.
  • Suitable for arbitrary Dω D ω and Wω W ω .
  • Easy to include arbitrary linear constraints.

Derivation

The weighted integral square error (or " L 2 L 2 error") is defined by

ε 2 =∫0πWωAω-Dω2dω ε 2 ω 0 W ω A ω D ω 2 (1)
where
  • Aω A ω : the actual amplitude response
  • Dω D ω : the ideal amplitude response
  • Wω W ω : nonnegative weighting function

The weighting function can be used to assign more importance to specific parts of the frequency response. For example, it is common to weight the stop-band more heavily than the pass-band. Once the length NN and the Type of the filter (I, II, III, IV) is chosen, the goal is to find the filter coefficients hn h n that minimizes ε 2 ε 2 . To develop this approach, we will consider the design of Type I FIR filters. The design of the other 3 linear-phase FIR filter types can be developed similarly.

Recall that for a Type I FIR filter, the amplitude response if given by Aω=∑n=0Mancosnω A ω n 0 M a n n ω where a0=hM a 0 h M an=2hM-n a n 2 h M n 1≤n≤M 1 n M M=N-12 M N 1 2

To obtain the coefficients an a n to minimize ε 2 ε 2 , we can set the derivatives equal to zero, ddak ε 2 =0 a k ε 2 0 0≤k≤M 0 k M The derivatives of ε 2 ε 2 with respect to ak a k can be found as:

ddak ε 2 =∫0πddakWωAω-Dω2dω=2∫0πWωAω-DωddakAωdω=2∫0πWωAω-Dωcoskωdω a k ε 2 ω 0 a k W ω A ω D ω 2 2 ω 0 W ω A ω D ω a k A ω 2 ω 0 W ω A ω D ω k ω (2)
Therefore, ddak ε 2 =0 a k ε 2 0 becomes ∫0πWωAωcoskωdω=∫0πWωDωcoskωdω ω 0 W ω A ω k ω ω 0 W ω D ω k ω or ∫0πWω∑n=0Mancosnωcoskωdω=∫0πWωDωcoskωdω ω 0 W ω n 0 M a n n ω k ω ω 0 W ω D ω k ω or ∑n=0Man∫0πWωcosnωcoskωdω=∫0πWωDωcoskωdω n 0 M a n ω 0 W ω n ω k ω ω 0 W ω D ω k ω If we define
Qkn=1π∫0πWωcosnωcoskωdω Q k n 1 ω 0 W ω n ω k ω (3)
and
bk=1π∫0πWωDωcoskωdω b k 1 ω 0 W ω D ω k ω (4)
then the derivative conditions can be written as
∑n=0MQknan=bk n 0 M Q k n a n b k (5)
where 0≤k≤M 0 k M This is a linear system of equations, and can be written in matrix form as Q00Q01…Q0MQ10Q11…Q1M⋮⋮⋮⋮QM0QM1…QMMa0a1⋮aM=b0b1⋮bM Q 0 0 Q 0 1 … Q 0 M Q 1 0 Q 1 1 … Q 1 M ⋮ ⋮ ⋮ ⋮ Q M 0 Q M 1 … Q M M a 0 a 1 ⋮ a M b 0 b 1 ⋮ b M or Qa=b Q a b Therefore, the Type 1 FIR filter that minimizes the square error can be obtained by solving this linear system of equations. a=Q-1b a Q b

Structure of the Matrix Q

It turns out that the matrix QQ has a special structure. Note that

cosnωcoskω=12cosk-nω+12cosk+nω n ω k ω 1 2 k n ω 1 2 k n ω (6)
so that
Qkn=1π∫0πWωcosnωcoskωdω=12π∫0πWωcosk-nωdω+12π∫0πWωcosk+nωdω Q k n 1 ω 0 W ω n ω k ω 1 2 ω 0 W ω k n ω 1 2 ω 0 W ω k n ω (7)
or
Qkn=12 Q 1 kn+12 Q 2 kn Q k n 1 2 Q 1 k n 1 2 Q 2 k n (8)
where Q 1 Q 1 and Q 2 Q 2 are defined as
Q 1 kn=1π∫0πWωcosk-nωdω Q 1 k n 1 ω 0 W ω k n ω (9)
and
Q 2 kn=1π∫0πWωcosk+nωdω Q 2 k n 1 ω 0 W ω k n ω (10)
Accordingly, we can write
Q 1 kn=qk-n Q 1 k n q k n (11)
Q 2 kn=qk+n Q 2 k n q k n (12)
where
qn=1π∫0πWωcosnωdω q n 1 ω 0 W ω n ω (13)
With this notation, the matrices Q 1 Q 1 and Q 2 Q 2 are written as
Q 1 =q0q1…qMq1q0…qM-1⋮⋮⋮⋮qMqM-1…q0 Q 1 q 0 q 1 … q M q 1 q 0 … q M 1 ⋮ ⋮ ⋮ ⋮ q M q M 1 … q 0 (14)
and
Q 2 =q0q1…qMq1q2…qM+1⋮⋮⋮⋮qMqM+1…q2M Q 2 q 0 q 1 … q M q 1 q 2 … q M 1 ⋮ ⋮ ⋮ ⋮ q M q M 1 … q 2 M (15)
Note that we have used q-n=qn q n q n here. The matrix Q 1 Q 1 is a symmetric Toeplitz matrix (constant along its diagonals), and the matrix Q 2 Q 2 is a Hankel matrix (constant along its anti-diagonals). Consequently,
  1. the matrices can be stored with less memory than arbitrary matrices ( 2M+1 2 M 1 numbers instead of M+12 M 1 2 numbers),
  2. there are fast algorithms to compute the solution to 'Toeplitz plus Hankel' systems with computational complexity OM2 O M 2 instead of OM3 O M 3 . (In fact, the complexity can be reduced further, but with higher overhead.)

Relation to the DTFT

To express qk q k and bk b k using the inverse Fourier transform, extend Dω D ω and Wω W ω symmetrically, so that D-ω=Dω D ω D ω and W-ω=Wω W ω W ω . Then we can write

qn=12π∫-ππWωcosnωdω q n 1 2 ω W ω n ω (16)
As sine is an anti-symmetric function
∫-ππWωsinnωdω=0 ω W ω n ω 0 (17)
so we can write
qn=12π∫-ππWωcosnω+ⅈsinnωdω q n 1 2 ω W ω n ω n ω (18)
or
qn=12π∫-ππWωⅇⅈnωdω q n 1 2 ω W ω n ω (19)
which we recognize as the inverse discrete-time Fourier transform
qn=DTFT-1 Wω q n DTFT W ω (20)
Similarly,
qn=DTFT-1 WωDω q n DTFT W ω D ω (21)

Low-Pass: Weighted Square Error

The weighting function Wω W ω can be used to improve the FIR low-pass filter because

  1. it allows you to eliminate Gibbs phenomenon by deleting a neighborhood around the band edge, and
  2. it allows you to assign different weights to the pass-band and the stop-band.
For example, if the ideal low-pass amplitude response is as shown in Figure 1.

Figure 1
Figure 1 (lowpass1.png)

and if the weighting function is as shown, Figure 2.

Figure 2
Figure 2 (lowpass2.png)

where ω p < ω o < ω s ω p ω o ω s , then the square error criterion becomes

ε 2 =∫0 ω p Aω-12dω+K∫ ω s πA2ωdω ε 2 ω 0 ω p A ω 1 2 K ω ω s A ω 2 (22)
To find the matrix QQ and the vector bb to this weighting function it is useful to recall
1π∫ ω 1 ω 2 cosnωdω= ω 2 πsinc ω 2 πn- ω 1 πsinc ω 1 πn 1 ω ω 1 ω 2 n ω ω 2 sinc ω 2 n ω 1 sinc ω 1 n (23)
Then
qk=1π∫0πWωcoskωdω=1π∫0 ω p coskωdω+Kπ∫ ω s πcoskωdω= ω p π+K1- ω s πifk=0 ω p πsinc ω p πk-K ω s πsinc ω s πkifk≠0 q k 1 ω 0 W ω k ω 1 ω 0 ω p k ω K ω ω s k ω ω p K 1 ω s k 0 ω p sinc ω p k K ω s sinc ω s k k 0 (24)
Similarly,
bk=1π∫0πWωDωcoskωdω=1π∫0 ω p coskωdω= ω p πsinc ω p πk b k 1 ω 0 W ω D ω k ω 1 ω 0 ω p k ω ω p sinc ω p k (25)

Weighted Low-Pass Example

In the following example, we design a Type 1 FIR low-pass filter of length 31, with band-edges ω p =0.26π ω p 0.26 , ω s =0.34π ω s 0.34 and a stop-band weight of K=10 K 10 . (See Figure 3)

Figure 3: The second figure shows the pass-band edge detail.
Figure 3 (wated.png)

Compared to the Impulse Response Truncation (IRT) method the peak error is reduced. The filter was designed using the following Matlab code. 


	  
	  % WEIGHTED LEAST SQUARE LOWPASS FILTER
	  
	  % filter length
	  N = 31;  
	  M = (N-1)/2;
	  
	  % set band-edges and stop-band weighting
	  wp = 0.26*pi;
	  ws = 0.34*pi;
	  K = 10;

	  % normalize band-edges for convenience
	  fp = wp/pi;
	  fs = ws/pi;

	  % construct q(k)
	  q  = [fp+K*(1-fs), fp*sinc(fp*[1:2*M])-K*fs*sinc(fs*[1:2*M])];

	  
	  % construct Q1, Q2, Q
	  Q1 = toeplitz(q([0:M]+1));
	  Q2 = hankel(q([0:M]+1),q([M:2*M]+1));
	  Q  = (Q1 + Q2)/2;

	  % construct b
	  b  = fp*sinc(fp*[0:M]');

	  % solve linear system to get a(n)
	  a  = Q\b;

	  % form impulse response h(n)
	  h  = [a(M+1:-1:2)/2; a(1); a(2:M+1)/2];

Comments, questions, feedback, criticisms?

Send feedback