Trade-off of Error Measures and Design Specifications In many filter design problems, more than one criterion is important. For example, both L2 and L∞ may be of interest in one filter. Often one is posed as a constraint and the other as an optimized variable. Indeed, because L2 approximation minimizes the error energy and because Parseval's theorem states that an optimal L2 frequency domain approximation is also an optimal L2 time domain approximation, an L∞ constrained minimum L2 error approximation seems a good practical approach. To see how this might have advantages, it is informative to examine the relationship of the L2 error to the L∞ error as the constraint is varied from tight to loose , in . From this one can see just how sensitive one error is to the other and how the traditional designs are extremes on this diagram.

The Squared Error vs. the Chebyshev Error for the Constrained Least Squared Error FIR Filter Another trade-off is the error in a Chebyshev design as a function of the transition band location. There are certain locations of transition band or band edges that give much lower ripple size than others. Rabiner has examined that relation , .

Constrained Least Squares Design There are problems where the peak error or Chebyshev error is important. This can be minimized directly using the Remez exchange algorithm but, in many cases, is better controlled by use of a peak error constraint on the basic least squared error formulation of the problem , , , . An efficient algorithm for minimizing the constrained least squared error uses Lagrange multipliers , and the Kuhn-Tucker conditions , . Similar to the Chebyshev design problem, there are two formulations of the problem: one where there is a well defined transition band separating the desired signal spectrum (passband) from the noise or interfering signal spectrum (stopband) and the second where there is a well defined frequency that separates the pass and stopband but no well defined transition band. The first case would include situations with signals residing in specified bands separated by “guard bands" such as commercial radio and TV transmissions. It also includes cases where due to multirate sampling, certain well defined bands are aliased into other well defined bands. The Parks-McClellan and Shpak-Antoniou Chebyshev designs address this case for the Chebyshev error. Adams' method , , , , , described below applies to the constrained least squares design with a specified transition band. The second case would include signals with known spectral support with additive white or broad-band noise. In these cases there is no obvious transition band or “don't care" band. The Hoffstetter-Oppenheim-Siegel and the method of address this case for a Chebyshev design. The method in section below applies to the constrained least squares design without a specified transition band.

The Lagrangian To pose the constrained least squared error optimization problem, we use a Lagrange multiplier formulation. First define the Lagrangian as L = P ∫ 0 π ( A ( ω ) - A d ( ω ) ) 2 d ω + ∑ i μ i A ( ω i ) - [ A d ( ω i ) ± T ( ω i ) ] where the μi are the necessary number of Langrange multipliers and P is a scale factor that can be chosen for simplicity later. The first term in () is the integral squared error of the frequency response to be minimized and the second term will be zero when the equality constraints are satisfied at the frequencies, ωi. The function T(ω) is the constraint function in that A(ω) must satisfy A d ( ω ) + T ( ω ) ≥ A ( ω ) ≥ A d ( ω ) - T ( ω ) . Necessary conditions for the minimization of the integral squared error are that the derivative of the Lagrangian with respect to the filter parameters a(n) defined in () and to the Lagrange multipliers μi be zero . The derivatives of the Lagrangian with respect to a(n) are d L d a ( n ) = P ∫ 0 π 2 ( A ( ω ) - A d ( ω ) ) d A d a d ω + ∑ i μ i d A d a ω i where from () we have for n=1,2,⋯,M d A ( ω ) d a ( n ) = cos ( ω n ) and for n=0 d A ( ω ) d a ( 0 ) = K . For n=1,2,⋯,M this gives d L d a ( n ) = 2 P ∫ A ( ω ) cos ( ω n ) d ω - ∫ A d ( ω ) cos ( ω n ) d ω + ∑ i μ i cos ( ω i n ) and for n=0 gives d L d a ( 0 ) = 2 P K ∫ A ( ω ) d ω - ∫ A d ( ω ) d ω + ∑ i μ i K . Using () for n=1,2,⋯,M, we have d L d a ( n ) = π P a ( n ) - a d ( n ) + ∑ i μ i cos ( ω i n ) = 0 and for n=0 d L d a ( 0 ) = 2 π P K 2 a ( 0 ) - a d ( 0 ) + K ∑ i μ i = 0 . Choosing P=1/π gives a ( n ) = a d ( n ) - ∑ i μ i cos ( ω i n ) and a ( 0 ) = a d ( 0 ) - 1 2 K ∑ i μ i Writing () and () in matrix form gives a = a d - H μ . where H is a matrix with elements h ( n , i ) = cos ( ω i n ) except for the first row which is h ( 0 , i ) = 1 2 K because of the normalization of the a(0) term. The ad(n) are the cosine coefficients for the unconstrained approximation to the ideal filter which result from truncating the inverse DTFT of Ad(ω). The derivative of the Lagrangian in () with respect to the Lagrange multipliers μi, when set to zero, gives A ( ω i ) = A d ( ω i ) ± T ( ω i ) = A c ( ω i ) which is simply a statement of the equality constraints. In terms of the filter's cosine coefficients a(n), from (), this can be written A c ( ω i ) = ∑ n a ( n ) cos ( ω i n ) + K a ( 0 ) and as matrices A c = G a where Ac is the vector of frequency response values which are the desired response plus or minus the constraints evaluated at the frequencies in the constraint set. The frequency response must interpolate these values. The matrix G is g ( i , n ) = cos ( ω i n ) except for the first column which is g ( i , 0 ) = K . Notice that if K=1/2, the first rows and columns are such that we have GT=H. The two equations () and () that must be satisfied can be written as a single matrix equation of the form I H G 0 a μ = a d A c or, if K=1/2, as I G T G 0 a μ = a d A c which have as solutions μ = ( G H ) - 1 ( G a d - A c ) a = a d - H μ The filter corresponding to the cosine coefficients a(n) minimize the L2 error norm subject the equality conditions in (). Notice that the term in () of the form Gad is the frequency response of the optimal unconstrained filter evaluated at the constraint set frequencies. Equation () could, therefore, be written μ = ( G H ) - 1 ( A u - A c )

The Constrained Weighted Least Squares Design of FIR Filters Combining the weighted least squared error formulation with the constrained least squared error gives the general formulation of this class of problems. We now modify the Lagrangian in () to allow a weighted squared error giving L = 1 π ∫ 0 π W ( ω ) ( A ( ω ) - A d ( ω ) ) 2 d ω + ∑ i μ i A ( ω i ) - A d ( ω i ) ± T ( ω i ) with a corresponding derivative of d L d a ( n ) = 2 π ∫ W ( ω ) ( A ( ω ) - A d ( ω ) d A d a d ω + ∑ i μ i d A d a ω i The integral cannot be carried out analytically for a general weighting function, but if the weight function is constant over each subband, () can be written d L d a ( n ) = 2 π ∑ k ∫ ω k ω k + 1 W k ( ∑ m = 1 M a ( m ) cos ( ω m ) + K a ( 0 ) - A d ( ω ) ) cos ( ω n ) d ω + ∑ i μ i d A d a ω i which after rearranging is = ∑ m = 1 M 2 π ∑ k W k ∫ ω k ω k + 1 ( cos ( ω m ) cos ( ω n ) ) d ω a ( m ) - 2 π ∑ k W k ∫ ω k ω k + 1 A d ( ω ) cos ( ω n ) d ω + ∑ i μ i cos ( ω i n ) = 0 where the integral in the first term can now be done analytically. In matrix notation () is R a - a d w + H μ = 0 This is a similar form to that in the multiband paper where the matrix R gives the effects of weighting with elements r ( n , m ) = 2 π ∑ k W k ∫ ω k ω k + 1 ( cos ( ω m ) cos ( ω n ) ) d ω except for the first row which should be divided by 2K because of the normalizing of the a(0) term in () and () and the first column which should be multiplied by K because of () and (). The matrix R is a sum of a Toeplitz matrix and a Hankel matrix and this fact might be used to advantage and adw is the vector of modified filter parameters with elements a d w ( n ) = 2 π ∑ W k ∫ ω k ω k + 1 A d ( ω ) cos ( ω n ) d ω and the matrix H is the same as used in () and defined in (). Equations () and () can be written together as a matrix equation R H G 0 a μ = a d w A c The solutions to () and () or to () are μ = ( G R - 1 H ) - 1 ( GR - 1 a d w - A c ) a = R - 1 ( a d w - H μ ) which are ideally suited to a language like Matlab and are implemented in the programs at the end of this book. Since the solution of Rau=adw is the optimal unconstrained weighted least squares filter, we can write () and () in the form μ = ( G R - 1 H ) - 1 ( G a u - A c ) = ( G R - 1 H ) - 1 ( A u - A c ) a = a u - R - 1 H μ a = a u - R - 1 H μ

The Exchange Algorithms This Lagrange multiplier formulation together with applying the Kuhn-Tucker conditions are used in an iterative multiple exchange algorithm similar to the Remez exchange algorithm to give the complete design method. One version of this exchange algorithm applies to the problem posed by Adams with specified pass and stopband edges and with zero error weighting in the transition band. This problem has the structure of a quadratic programming problem and could be solved using general QP methods but the multiple exchange algorithm suggested here is probably faster. The second version of this exchange algorithm applies to the problem where there is no explicitly specified transition band. This problem is not strictly a quadratic programming problem and our exchange algorithm has no proof of convergence (the HOS algorithm also has no proof of convergence). However, in practice, this program has proven to be robust and converges for a wide variety of lengths, constraints, weights, and band edges. The performance is completely independent of the normalizing parameter K. Notice that the inversion of the R matrix is done once and does not have to be done each iteration. The details of the program are included in the filter design paper and in the Matlab program at the end of this book. At mentioned earlier, this design problem might be addressed by general constrained quadratic programming methods , , , , , , , , , .

Examples and Observations on CLS Designs Here we show that the CLS FIR filter design approach is probably the best general FIR filter design method. For example, a length-31 linear phase lowpass FIR filter is designed for a band edge of 0.3 and the constraint that the response in the stop cannot be greater than 0.03 is illustrated in .

Response of a Constrained Least Squared Error Filter Design This filter was designed using the Matlab command: 'fircls1()' function.

Approximation and the Iterative Reweighted Least Squares Method We now consider the general Lp approximation which contains the least squares L2 and the Chebyshev L∞ cases. This approach is described in , , using the iterative reweighted least squared (IRLS) alorithm and looks attractive in that it can use different p in different frequency bands. This would allow, for example, a least squared error approximation in the passband and a Chebyshev approximation in the stopband. The IRLS method can also be used for complex Chebyshev approximations and constrained L2 approximatin.

Iterative Reweighted Least Squares Filter Design Methods There are cases where it is desirable to design an FIR filter that will minimize the Lp error norm. The error is defined by q = ∫ Ω | A ( ω ) - A d ( ω ) | p d ω but we usually work with Q2. For large p, the results are essentially the same as the Chebyshev filter and this gives a continuum of design between L2 and L∞. It also allows the very interesting important possibility of allowing p(ω) to be a function of frequency. This means one could have different error criteria in different frequency bands. It can be modified to give the same effects as a constraint. This approach is discussed in . It can be applied to complex approximation and to two-dimensional filter design , . The least squared error and the minimum Chebyshev error criteria are the two most commonly used linear-phase FIR filter design methods . There are many situations where better total performance would be obtained with a mixture of these two error measures or some compromise design that would give a trade-off between the two. We show how to design a filter with an L2 approximation in the passband and a Chebyshev approximation in the stopband. We also show that by formulating the Lp problem we can solve the constrained L2 approximation problem . This section first explores the minimization of the pth power of the error as a general approximation method and then shows how this allows L2 and L∞ approximations to be used simultaneous in different frequency bands of one filter and how the method can be used to impose constraints on the frequency response. There are no analytical methods to find this approximation, therefore, an iterative method is used over samples of the error in the frequency domain. The methods developed here , are based on what is called an iterative reweighted least squared (IRLS) error algorithm , , and they can solve certain FIR filter design problems that neither the Remez exchange algorithm nor analytical L2 methods can. The idea of using an IRLS algorithm to achieve a Chebyshev or L∞ approximation seems to have been first developed by Lawson and extended to Lp by Rice and Usow , . The basic IRLS method for Lp was given by Karlovitz and extended by Chalmers, et. al. , Bani and Chalmers , and Watson . Independently, Fletcher, Grant and Hebden developed a similar form of IRLS but based on Newton's method and Kahng did likewise as an extension of Lawson's algorithm. Others analyzed and extended this work , , , . Special analysis has been made for 1≤p<2 by , , , , , , and for p=∞ by , , , , , . Relations to the Remez exchange algorithm , were suggested by , to homotopy by , and to Karmarkar's linear programming algorithm by , . Applications of Lawson's algorithm to complex Chebyshev approximation in FIR filter design have been made in , , , and to 2-D filter design in . Reference indicates further results may be forthcoming. Application to array design can be found in and to statistics in . This paper unifies and extends the IRLS techniques and applies them to the design of FIR digital filters. It develops a framework that relates all of the above referenced work and shows them to be variations of a basic IRLS method modified so as to control convergence. In particular, we generalize the work of Rice and Usow on Lawson's algorithm and explain why its asymptotic convergence is slow. The main contribution here is a new robust IRLS method , that combines an improved convergence acceleration scheme and a Newton based method. This gives a very efficient and versatile filter design algorithm that performs significantly better than the Rice-Usow-Lawson algorithm or any of the other IRLS schemes. Both the initial and asymptotic convergence behavior of the new algorithm is examined and the reason for occasional slow convergence of this and all other IRLS methods is discovered. We then show that the new IRLS method allows the use of p as a function of frequency to achieve different error criteria in the pass and stopbands of a filter. Therefore, this algorithm can be applied to solve the constrained Lp approximation problem. Initial results of applications to the complex and two-dimensional filter design problem are presented. Although the traditional IRLS methods were sometimes slower than competing approaches, the results of this paper and the availability of fast modern desktop computers make them practical now and allow exploitation of their greater flexibility and generality.

Minimum Squared Error Approximations Various approximation methods can be developed by considering different definitions of norm or error measure. Commonly used definitions are L1, L2, and Chebyshev or L∞. Using the L2 norm, gives the scalar error to minimize q = ∑ k = 0 L - 1 | A ( ω k ) - A d ( ω k ) | 2 or in matrix notation using (), the error or residual vector is defined by q = C a - A d giving the scalar error of () as q = ϵ T ϵ . This can be minimized by solution of the normal equations , , C T C a = C T A d . The weighted squared error defined by q = ∑ k = 0 L - 1 w k 2 | A ( ω k ) - A d ( ω k ) | 2 . or, in matrix notation using () and () causes () to become q = ϵ T W T W ϵ which can be minimized by solving W C a = W A d with the normal equations C T W T W C a = C T W T W A d where W is an L by L diagonal matrix with the weights wk from () along the diagonal. A more general formulation of the approximation simply requires WTW to be positive definite. Some authors define the weighted error in () using wk rather than wk2. We use the latter to be consistent with the least squared error algorithms in Matlab . Solving () is a direct method of designing an FIR filter using a weighted least squared error approximation. To minimize the sum of the squared error and get approximately the same result as minimizing the integral of the squared error, one must choose L to be 3 to 10 or more times the length L of the filter being designed.

Iterative Algorithms to Minimize the Error There is no simple direct method for finding the optimal approximation for any error power other than two. However, if the weighting coefficients wk as elements of W in () could be set equal to the elements in |A-Ad|, minimizing () would minimize the fourth power of |A-Ad|. This cannot be done in one step because we need the solution to find the weights! We can, however, pose an iterative algorithm which will first solve the problem in () with no weights, then calculate the error vector ϵ from () which will then be used to calculate the weights in (). At each stage of the iteration, the weights are updated from the previous error and the problem solved again. This process of successive approximations is called the iterative reweighted least squared error algorithm (IRLS). The basic IRLS equations can also be derived by simply taking the gradient of the p-error with respect to the filter coefficients h or a and setting it equal to zero , . These equations form the basis for the iterative algorithm. If the algorithm is a contraction mapping , the successive approximations will converge and the limit is the solution of the minimum L4 approximation problem. If a general problem can be posed , , as the solution of an equation in the form x = f ( x ) , a successive approximation algorithm can be proposed which iteratively calculates x using x m + 1 = f ( x m ) starting with some x0. The function f(·) maps xm into xm+1 and, if limm→∞xm=x0 where x0=f(x0), x0 is the fixed point of the mapping and a solution to (). The trick is to find a mapping that solves the desired problem, converges, and converges fast. By setting the weights in () equal to w ( k ) = | A ( ω k ) - A d ( ω k ) | ( p - 2 ) / 2 , the fixed point of a convergent algorithm minimizes q = ∑ k = 0 L - 1 | A ( ω k ) - A d ( ω k ) | p . It has been shown that weights always exist such that minimizing () also minimizes (). The problem is to find those weights efficiently.

Basic Iterative Reweighted Least Squares The basic IRLS algorithm is started by initializing the weight matrix defined in () and () for unit weights with W0=I. Using these weights to start, the mth iteration solves () for the filter coefficients with a m = [ C T W m T W m C ] - 1 C T W m T W m A d This is a formal statement of the operation. In practice one should not invert a matrix, one should use a sophisticated numerical method to solve the overdetermined equations in () The error or residual vector () for the mth iteration is found by ϵ m = C a m - A d A new weighting vector is created from this error vector using () by w m + 1 = | ϵ m | ( p - 2 ) / 2 whose elements are the diagonal elements of the new weight matrix W m + 1 = d i a g [ w m + 1 ] . Using this weight matrix, we solve for the next vector of filter coefficients by going back to () and this defines the basic iterative process of the IRLS algorithm. It can easily be shown that the a that minimizes () is a fixed point of this iterative map. Unfortunately, applied directly, this basic IRLS algorithm does not converge and/or it has numerical problems for most practical cases . There are three aspects that must be addressed. First, the IRLS algorithm must theoretically converge. Second, the solution of () must be numerically stable. Finally, even if the algorithm converges and is numerically stable, it must converge fast enough to be practical. Both theory and experience indicate there are different convergence problems connected with several different ranges and values of p. In the range 2≤p<3, virtually all methods converge , , . In the range 3≤p<∞, the algorithm diverges and the various methods discussed in this paper must be used. As p becomes large compared to 2, the weights carry a larger contribution to the total minimization than the underlying least squared error minimization, the improvement at each iteration becomes smaller, and the likelihood of divergence becomes larger. For p=∞ we can use to advantage the fact that the optimal approximation solution to () is unique but the weights in () that give that solution are not. In other words, different matrices W give the same solution to () but will have different convergence properties. This allows certain alteration to the weights to improve convergence without harming the optimality of the results . In the range 1<p<2, both convergence and numerical problems exist as, in contrast to p>2, the IRLS iterations are undoing what the underlying least squares is doing. In particular, the weights near frequencies with small errors become very large. Indeed, if the error happens to be zero, the weight becomes infinite because of the negative exponent in (). For p=1 the solution to the optimization problem is not even unique. The various algorithms that are presented below are based on schemes to address these problems.

The Karlovitz Method In order to achieve convergence, a second order update is used which only partially changes the filter coefficients am in () each iteration. This is done by first calculating the unweighted L2 approximation filter coefficients using () as a 0 = [ C T C ] - 1 C T A d . The error or residual vector () for the mth iteration is found as () by ϵ m = C a m - A d and the new weighting vector is created from this error vector using () by w m + 1 = | ϵ m | ( p - 2 ) / 2 whose elements are the diagonal elements of the new weight matrix W m + 1 = d i a g [ w m + 1 ] . This weight matrix is then used to calculate a temporary filter coefficient vector by a ^ m + 1