# Connexions

You are here: Home » Content » Linear phase l_p filter design

### Recently Viewed

This feature requires Javascript to be enabled.

# Linear phase l_p filter design

Module by: Ricardo Vargas. E-mail the author

Linear phase FIR filters are important tools in signal processing. As will be shown below, they do not require the user to specify a phase response in their design (since the assumption is that the desired phase response is indeed linear). Besides, they satisfy a number of symmetry properties that allow for the reduction of dimensions in the optimization process, making them easier to design computationally. Finally, there are applications where a linear phase is desired as such behavior is more physically meaningful.

## Four types of linear phase filters

The frequency response of an FIR filter h(n)h(n) is given by

H ( ω ) = n = 0 N - 1 h ( n ) e - j ω n H ( ω ) = n = 0 N - 1 h ( n ) e - j ω n
(1)

In general, H(ω)=R(ω)+jI(ω)H(ω)=R(ω)+jI(ω) is a periodic complex function of ωω (with period 2π2π). Therefore it can be written as follows,

H ( ω ) = R ( ω ) + j I ( ω ) = A ( ω ) e j φ ( ω ) H ( ω ) = R ( ω ) + j I ( ω ) = A ( ω ) e j φ ( ω )
(2)

where the magnitude response is given by

A ( ω ) = | H ( ω ) | = R ( ω ) 2 + I ( ω ) 2 A ( ω ) = | H ( ω ) | = R ( ω ) 2 + I ( ω ) 2
(3)

and the phase response is

φ ( ω ) = sin I ( ω ) R ( ω ) φ ( ω ) = sin I ( ω ) R ( ω )
(4)

However A(ω)A(ω) is not analytic and φ(ω)φ(ω) is not continuous. From a computational point of view Equation 2 would have better properties if both A(ω)A(ω) and φ(ω)φ(ω) were continuous analytic functions of ωω; an important class of filters for which this is true is the class of linear phase filters


[11].

Linear phase filters have a frequency response of the form

H ( ω ) = A ( ω ) e j φ ( ω ) H ( ω ) = A ( ω ) e j φ ( ω )
(5)

where A(ω)A(ω) is the real, continuous amplitude response of H(ω)H(ω) and

φ ( ω ) = K 1 + K 2 ω φ ( ω ) = K 1 + K 2 ω
(6)

is a linear phase function in ωω (hence the name); K1K1 and K2K2 are constants. Figure 1 shows the frequency response for a linear phase FIR filter. The jumps in the phase response correspond to sign reversals in the magnitude resulting as defined in Equation 3.

Consider a length-NN FIR filter (assume for the time being that NN is odd). Its frequency response is given by

H ( ω ) = n = 0 N - 1 h ( n ) e - j ω n = e - j ω M n = 0 2 M h ( n ) e j ω ( M - n ) H ( ω ) = n = 0 N - 1 h ( n ) e - j ω n = e - j ω M n = 0 2 M h ( n ) e j ω ( M - n )
(7)

where M=N-12M=N-12. Equation Equation 7 can be written as follows,

H ( ω ) = e - j ω M [ h ( 0 ) e j ω M + ... + h ( M - 1 ) e j ω + h ( M ) + h ( M + 1 ) e - j ω + ... + h ( 2 M ) e - j ω M ] H ( ω ) = e - j ω M [ h ( 0 ) e j ω M + ... + h ( M - 1 ) e j ω + h ( M ) + h ( M + 1 ) e - j ω + ... + h ( 2 M ) e - j ω M ]
(8)

It is clear that for an odd-length FIR filter to have the linear phase form described in Equation 5, the term inside braces in Equation 8 must be a real function (thus becoming A(ω)A(ω)). By imposing even symmetry on the filter coefficients about the midpoint (n=Mn=M), that is

h ( k ) = h ( 2 M - k ) h ( k ) = h ( 2 M - k )
(9)

Equation 8 becomes

H ( ω ) = e - j ω M h ( M ) + 2 n = 0 M - 1 h ( n ) cos ω ( M - n ) H ( ω ) = e - j ω M h ( M ) + 2 n = 0 M - 1 h ( n ) cos ω ( M - n )
(10)

Similarly, with odd symmetry (i.e. h(k)=h(2M-k)h(k)=h(2M-k)) equation Equation 8 becomes

H ( ω ) = e j ( π 2 - ω M ) 2 n = 0 M - 1 h ( n ) tan -1 ω ( M - n ) H ( ω ) = e j ( π 2 - ω M ) 2 n = 0 M - 1 h ( n ) tan -1 ω ( M - n )
(11)

Note that the term h(M)h(M) disappears as the symmetry condition requires that

h ( M ) = h ( N - M - 1 ) = - h ( M ) = 0 h ( M ) = h ( N - M - 1 ) = - h ( M ) = 0
(12)

Similar expressions can be obtained for an even-length FIR filter,

H ( ω ) = n = 0 N - 1 h ( n ) e - j ω n = e - j ω M n = 0 N 2 - 1 h ( n ) e j ω ( M - n ) H ( ω ) = n = 0 N - 1 h ( n ) e - j ω n = e - j ω M n = 0 N 2 - 1 h ( n ) e j ω ( M - n )
(13)

It is clear that depending on the combinations of NN and the symmetry of h(n)h(n), it is possible to obtain four types of filters [11], [10], [5]. Table 1 shows the four possible linear phase FIR filters described by Equation 5.

 N Odd Even Symmetry A ( ω ) = h ( M ) + 2 ∑ n = 0 M - 1 h ( n ) ⋅ cos ω ( M - n ) φ ( ω ) = - ω M A ( ω ) = h ( M ) + 2 ∑ n = 0 M - 1 h ( n ) ⋅ cos ω ( M - n ) φ ( ω ) = - ω M Odd Symmetry A ( ω ) = 2 ∑ n = 0 M - 1 h ( n ) ⋅ sin ω ( M - n ) φ ( ω ) = π 2 - ω M A ( ω ) = 2 ∑ n = 0 M - 1 h ( n ) ⋅ sin ω ( M - n ) φ ( ω ) = π 2 - ω M N Even Even Symmetry A ( ω ) = h ( M ) + 2 ∑ n = 0 N 2 - 1 h ( n ) ⋅ cos ω ( M - n ) φ ( ω ) = - ω M A ( ω ) = h ( M ) + 2 ∑ n = 0 N 2 - 1 h ( n ) ⋅ cos ω ( M - n ) φ ( ω ) = - ω M Odd Symmetry d A ( ω ) = 2 ∑ n = 0 N 2 - 1 h ( n ) ⋅ sin ω ( M - n ) φ ( ω ) = π 2 - ω M d A ( ω ) = 2 ∑ n = 0 N 2 - 1 h ( n ) ⋅ sin ω ( M - n ) φ ( ω ) = π 2 - ω M

## IRLS-based methods

Section 1 introduced linear phase filters in detail. In this section we cover the use of IRLS methods to design linear phase FIR filters according to the lplp optimality criterion. Recall from Section 1 that for any of the four types of linear phase filters their frequency response can be expressed as

H ( ω ) = A ( ω ) e j ( K 1 + K 2 ω ) H ( ω ) = A ( ω ) e j ( K 1 + K 2 ω )
(14)

Since A(ω)A(ω) is a real continuous function as defined by Table 1, one can write the linear phase lplp design problem as follows

min a D ( ω ) - A ( ω ; a ) p p min a D ( ω ) - A ( ω ; a ) p p
(15)

where aa relates to hh by considering the symmetry properties outlined in Table 1. Note that the two objects from the objective function inside the lplp norm are real. By sampling Equation 15 one can write the design problem as follows

min a k | D ( ω k ) - A ( ω k ; a ) | p min a k | D ( ω k ) - A ( ω k ; a ) | p
(16)

or

min a k | D k - C k a | p min a k | D k - C k a | p
(17)

where DkDk is the kk-th element of the vector DD representing the sampled desired frequency response D(ωk)D(ωk), and CkCk is the kk-th row of the trigonometric kernel matrix as defined by Table 1.

One can apply the basic IRLS approach described in (Reference) to solve Equation 17 by posing this problem as a weighted least squares one:

min a k w k | D k - C k a | 2 min a k w k | D k - C k a | 2
(18)

The main issue becomes iteratively finding suitable weights ww for Equation 18 so that the algorithm converges to the optimal solution aa of the lplp problem Equation 15. Existence of adequate weights is guaranteed by Theorem (Reference) as presented in (Reference); finding these optimal weights is indeed the difficult part. Clearly a reasonable choice for ww is that which turns Equation 18 into Equation 17, namely

w = | D - C a | p - 2 w = | D - C a | p - 2
(19)

Therefore the basic IRLS algorithm for problem Equation 17 would be:

1. Initialize the weights w0w0 (a reasonable choice is to make them all equal to one).
2. At the ii-th iteration the solution is given by
ai+1=[CTWiTWiC]-1CTWiTWiDai+1=[CTWiTWiC]-1CTWiTWiD
(20)
3. Update the weights with
wi+1=|D-Cai+1|p-2wi+1=|D-Cai+1|p-2
(21)
4. Repeat the last steps until convergence is reached.

It is important to note from Appendix (Reference) that Wi=diag(wi)Wi=diag(wi). In practice it has been found that this approach has practical defficiencies, since the inversion required by Equation 20 often leads to an ill-posed problem and, in most cases, convergence is not achieved.

As mentioned before, the basic IRLS method has drawbacks that make it unsuitable for practical implementations. Charles Lawson considered a version of this algorithm applied to the solution of ll problems (for details refer to [8]). His method has linear convergence and is prone to problems with proportionately small residuals that could lead to zero weights and the need for restarting the algorithm. In the context of lplp optimization, Rice and Usow [12] built upon Lawson's method by adapting it to lplp problems. Like Lawson's methods, the algorithm by Rice and Usow updates the weights in a multiplicative manner; their method shares similar drawbacks with Lawson's. Rice and Usow defined

w i + 1 ( ω ) = w i α ( ω ) | ϵ i ( ω ) | β w i + 1 ( ω ) = w i α ( ω ) | ϵ i ( ω ) | β
(22)

where

α = γ ( p - 2 ) γ ( p - 2 ) + 1 α = γ ( p - 2 ) γ ( p - 2 ) + 1
(23)

and

β = α 2 γ = p - 2 2 γ ( p - 2 ) + 2 β = α 2 γ = p - 2 2 γ ( p - 2 ) + 2
(24)

L. A. Karlovitz realized the computational problems associated with the basic IRLS method and improved on it by partially updating the filter coefficient vector. He defines

a ^ i + 1 = [ C T W i T W i C ] - 1 C T W i T W i D a ^ i + 1 = [ C T W i T W i C ] - 1 C T W i T W i D
(25)

and uses a^a^ in

a i + 1 = λ a ^ i + 1 + ( 1 - λ ) a i a i + 1 = λ a ^ i + 1 + ( 1 - λ ) a i
(26)

where λ[0,1]λ[0,1] is a partial step parameter that must be adjusted at each iteration. Karlovitz's method [7] has been shown to converge globally for even values of pp (where 2p<2p<). In practice, convergence problems have been found even under such assumptions. Karlovitz proposed the use of line searches to find the optimal value of λλ at each iteration, which basically creates an independent optimization problem nested inside each iteration of the IRLS algorithm. While computationally this search process for the optimal λλ makes Karlovitz's method impractical, his work indicates the feasibility of IRLS methods and proves that partial updating indeed overcomes some of the problems in the basic IRLS method. Furthermore, Karlovitz's method is the first one to depart from a multiplicative updating of the weights in favor of an additive updating on the filter coefficients. In this way some of the problems in the Lawson-Rice-Usow approach are overcome, especially the need for restarting the algorithm.

S. W. Kahng built upon the findings by Karlovitz by considering the process of finding an adequate λλ for partial updating. He applied Newton-Raphson's method to this problem and proposed a closed form solution for λλ, given by

λ = 1 p - 1 λ = 1 p - 1
(27)

resulting in

a i + 1 = λ a ^ i + 1 + ( 1 - λ ) a i a i + 1 = λ a ^ i + 1 + ( 1 - λ ) a i
(28)

The rest of Kahng's algorithm follows Karlovitz's approach. However, since λλ is fixed, there is no need to perform the linear search at each iteration. Kahng's method has an added benefit: since it uses Newton's method to find λλ, the algorithm tends to converge much faster than previous approaches. It has indeed been shown to converge quadratically. However, Newton-Raphson-based algorithms are not guaranteed to converge globally unless at some point the existing solution lies close enough to the solution, within their radius of convergence [1]. Fletcher, Grant and Hebden[6] derived the same results independently.

Burrus, Barreto and Selesnick [4], [2], [3] modified Kahng's methods in several important ways in order to improve on their initial and final convergence rates and the method's stability (we refer to this method as BBS). The first improvement is analogous to a homotopy [9]. Up to this point all efforts in lplp filter design attempted to solve the actual


lplp problem from the first iteration. In general there is no reason to believe that an initial guess derived from an unweighted l2l2 formulation (that is, the l2l2 design that one would get by setting w0=1^w0=1^) will look in any way similar to the actual lplp solution that one is interested in. However it is known that there exists a continuity of lplp solutions for 1<p<1<p<. In other words, if a2a2 is the optimal l2l2 solution, there exists a pp for which the optimal lplp solution apap is arbitrarily close to a2a2; that is, for a given δ>0δ>0

a 2 - a p δ for some p ( 2 , ) a 2 - a p δ for some p ( 2 , )
(29)

This fact allows anyone to gradually move from an lplp solution to an lqlq solution.

To accelerate initial convergence, the BBS method of Burrus et al. initially solves for l2l2 by setting p0=2p0=2 and then sets pi=σ·pi-1pi=σ·pi-1, where σσ is a convergence parameter defined by 1σ21σ2. Therefore at the ii-th iteration

p i = min ( p d e s , σ p i - 1 ) p i = min ( p d e s , σ p i - 1 )
(30)

where pdespdes corresponds to the desired lplp solution. The implementation of each iteration follows Karlovitz's method with Kahng's choice of λλ, using the particular selection of pp given by Equation 30.

To summarize, define the class of IRLS algorithms as follows: after ii iterations, given a vector aiai the IRLS iteration requires two steps,

1. Find wi=f(ai)wi=f(ai)
2. Find ai+1=g(wi,ai)ai+1=g(wi,ai)

The following is a summary of the IRLS-based algorithms discussed so far and their corresponding updating functions:

1. Basic IRLS algorithm.
• wi=|D-Cai|p-2wi=|D-Cai|p-2
• Wi=diag(wi)Wi=diag(wi)
• ai+1=CTWiTWiC-1CTWiTWiDai+1=CTWiTWiC-1CTWiTWiD
2. Rice-Usow-Lawson (RUL) method
• wi=wi-1α|D-Cai|α2γwi=wi-1α|D-Cai|α2γ
• Wi=diag(wi)Wi=diag(wi)
• ai+1=CTWiTWiC-1CTWiTWiDai+1=CTWiTWiC-1CTWiTWiD
• α=γ(p-2)γ(p-2)+1α=γ(p-2)γ(p-2)+1
• γγ constant
3. Karlovitz' method
• wi=|D-Cai|p-2wi=|D-Cai|p-2
• Wi=diag(wi)Wi=diag(wi)
• ai+1=λCTWiTWiC-1CTWiTWiD+(1-λ)aiai+1=λCTWiTWiC-1CTWiTWiD+(1-λ)ai
• λλ constant
4. Kahng's method
• wi=|D-Cai|p-2wi=|D-Cai|p-2
• Wi=diag(wi)Wi=diag(wi)
• ai+1=1p-1CTWiTWiC-1CTWiTWiD+p-2p-1aiai+1=1p-1CTWiTWiC-1CTWiTWiD+p-2p-1ai
5. BBS method
• pi=min(pdes,σ·pi-1)pi=min(pdes,σ·pi-1)
• wi=|D-Cai|pi-2wi=|D-Cai|pi-2
• Wi=diag(wi)Wi=diag(wi)
• ai+1=1pi-1CTWiTWiC-1CTWiTWiD+pi-2pi-1aiai+1=1pi-1CTWiTWiC-1CTWiTWiD+pi-2pi-1ai
• σσ constant

Much of the performance of a method is based upon whether it can actually converge given a certain error measure. In the case of the methods described above, both convergence rate and stability play an important role in their performance. Both Karlovitz and RUL methods are supposed to converge linearly, while Kahng's and the BBS methods converge quadratically, since they both use a Newton-based additive update of the weights.

Barreto showed in [2] that the modified version of Kahng's method (or BBS) typically converges faster than the RUL algorithm. However, this approach presents some peculiar problems that depend on the transition bandwidth ββ. For some particular values of ββ, the BBS method will result in an ill-posed weight matrix that causes the lplp error to increase dramatically after a few iterations as illustrated in Figure 2 (where f=ω/2πf=ω/2π).

Two facts can be derived from the examples in Figure 2: for this particular bandwidth the error increased slightly after the fifth and eleventh iterations, and increased dramatically after the sixteenth. Also, it is worth to notice that after such increase, the error started to decrease quadratically and that, at a certain point, the error became flat (thus reaching the numerical accuracy limits of the digital system).

The effects of different values of σσ were studied to find out if a relationship between σσ and the error increase could be determined. Figure 3 shows the lplp error for different values of ββ and for σ=1.7σ=1.7. It can be seen that some particular bandwidths cause the algorithm to produce a very large error.

Our studies (as well as previous work from J. A. Barreto [2]) demonstrate that this error explosion occurs only for a small range of bandwidth specifications. Under most circumstances the BBS method exhibits fast convergence properties to the desired solution. However at this point it is not understood what causes the error increase and therefore this event cannot be anticipated. In order to avoid such problem, I propose the use of an adaptive scheme that modifies the BBS step. As pp increases the step from a current lplp guess to the next also increases, as described in Equation 30. In other words, at the ii-th iteration one approximates the l2σil2σi solution (as long as the algorithm has not yet reached the desired pp); the next iteration one approximates l2σi+1l2σi+1. There is always a possibility that these two solutions lie far apart enough that the algorithm takes a descent step so that the l2σi+1l2σi+1guess is too far away from the actual l2σi+1l2σi+1


solution. This is better illustrated in Figure 4.

The conclusions derived above suggest the possibility to use an adaptive algorithm [13] that changes the value of σσ so that the error always decreases. This idea was implemented by calculating temporary new weight and filter coefficient vectors that will not become the updated versions unless their resulting error is smaller than the previous one. If this is not the case, the algorithm "tries" two values of σσ, namely

σ L = σ * ( 1 - δ ) and σ H = σ * ( 1 + δ ) σ L = σ * ( 1 - δ ) and σ H = σ * ( 1 + δ )
(31)

(where δδ is an updating variable). The resulting errors for each attempt are calculated, and σσ is updated according to the value that produced the smallest error. The error of this new σσ is compared to the error of the nonupdated weights and coefficients, and if the new σσ produces a smaller error, then such vectors are updated; otherwise another update of σσ is performed. The modified adaptive IRLS algorithm can be summarized as follows,

1. Find the unweighted approximation a0=CTC-1CTDa0=CTC-1CTD and use p0=2σp0=2σ (with 1σ21σ2)
2. Iteratively solve Equation 25 and Equation 26 using λi=1pi-1λi=1pi-1 and find the resulting error εiεi for the ii-th iteration
3. If εiεi-1εiεi-1,
• Calculate Equation 31
• Select the smallest of εσLεσL and εσHεσH to compare it with εiεi until a value is found that results in a decreasing error
Otherwise iterate as in the BBS algorithm.

The algorithm described above changes the value of σσ that causes the algorithm to produce a large error. The value of σσ is updated as many times as necessary without changing the values of the weights, the filter coefficients, or pp. If an optimal value of σσ exists, the algorithm will find it and continue with this new value until another update in σσ becomes necessary.

The algorithm described above was implemented for several combinations of σσ and ββ; for all cases the new algorithm converged faster than the BBS algorithm (unless the values of σσ and ββ are such that the error never increases). The results are shown in Figure 5.a for the specifications from Figure 2. Whereas using the BBS method for this particular case results in a large error after the sixteenth iteration, the adaptive method converged before ten iterations.

Figure 5.b illustrates the change of σσ per iteration in the adaptive method, using an update factor of δ=0.1δ=0.1. The lplp error stops decreasing after the fifth iteration (where the BBS method introduces the large error); however, the adaptive algorithm adjusts the value of σσ so that the lplp error continues decreasing. The algorithm decreased the initial value of σσ from 1.75 to its final value of 1.4175 (at the expense of only one additional iteration with σ=1.575σ=1.575).

One result worth noting is the relationship between l2l2 and ll solutions and how they compare to lplp designs. Figure 6 shows a comparison of designs for a length-21 Type-I linear phase low pass FIR filter with transition band defined by f={0.2,0.24}f={0.2,0.24}. The curve shows the l2l2 versus ll errors (namely ε2ε2 and εε); the values of pp used to make this curve were p={2,2.2,2.5,3,4,5,7,10,15,20,30,50,60,100,150,200,400,}p={2,2.2,2.5,3,4,5,7,10,15,20,30,50,60,100,150,200,400,} (Matlab's firls and firpm functions were used to design the l2l2 and ll filters respectively). Note the very small decrease in εε after pp reaches 100. The curve suggests that a better compromise between ε2ε2 and εε can be reached by choosing 2<p<2<p<. Furthermore, to get better results one can concentrate on values between p=5p=5 and p=20p=20; fortunately, for values of pp so low no numerical complications arise and convergence is reached in a few iterations.

## References

1. Aoki, Masanao. (1971). Introduction to Optimization Techniques. The Macmillan Company.
2. Barreto, Jose A. (1992). Approximation by the Iterative Reweighted Least Squares Method and the Design of Digital FIR Filters in One Dimension. Masters thesis. Rice University.
3. Burrus, C. S. and Barreto, J. A. (1992, May). Least -power Error Design of FIR Filters. In Proc. IEEE Int. Symp. Circuits, Syst. ISCAS-92. (Vol. 2, pp. 545-548). San Diego, CA
4. Burrus, C. S. and Barreto, J. A. and Selesnick, I. W. (1994, November). Iterative Reweighted Least-Squares Design of FIR Filters. IEEE Transactions on Signal Processing, 42(11), 2926-2936.
5. Burrus, Charles S. and McClellan, James H. (1994). Computer-based Exercises for Signal Processing. Prentice Hall.
6. Fletcher, R. and Grant, J. A. and Hebden, M. D. (1972, Apr). The Calculation of Linear Best Approximations. The Computer Journal, 14(118), 276-279.
7. Karlovitz, L. A. (1970). Construction of Nearest Points in the , even and norms, I. Journal of Approximation Theory, 3, 123-127.
8. Lawson, C. L. (1961). Contributions to the Theory of Linear Least Maximum Approximations. Ph. D. Thesis. UCLA.
9. Nocedal, Jorge and Wright, Stephen J. (1999). Springer series in operations research: Numerical Optimization. New York, NY: Springer-Verlag.
10. Oppenheim, A. V. and Schafer, R. W. (1989). Discrete-Time Signal Processing. Englewood Cliffs, N.J.: Prentice-Hall.
11. Parks, T. W. and Burrus, C. S. (1987). Digital Filter Design. John Wiley and Sons.
12. Rice, John R. and Usow, Karl H. (1968). The Lawson Algorithm and Extensions. Mathematics of Computation, 22, 118-127.
13. Vargas, Ricardo A. and Burrus, Charles S. (1999). Adaptive Iterative Reweighted Least Squares Design of FIR Filters. Proc. ICASSP, 3, Signal Processing, Theory and Methods, 1129-32.

## Content actions

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

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

##### What is in a lens?

Lens makers point to 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 member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks