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

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,

where the magnitude response is given by

and the phase response is

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

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

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

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

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

Equation 8 becomes

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

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

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

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.

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

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

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

or

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:

The main issue becomes iteratively finding suitable weights ww for Equation 18 so that the algorithm converges to the optimal solution a⋆ a⋆ 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

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

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 l∞l∞ 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

where

and

and follow the basic algorithm.

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

and uses a^a^ in

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 2≤p<∞2≤p<∞). 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

resulting in

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 a2⋆ a2⋆ is the optimal l2l2 solution, there exists a pp for which the optimal lplp solution ap⋆ ap⋆ is arbitrarily close to a2⋆ a2⋆ ; that is, for a given δ>0δ>0

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

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,

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

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

(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,

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 l∞l∞ 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 l∞l∞ 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 l∞l∞ 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.