WLS Design of Sparse FIR Digital Filters - Semantic Scholar
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
125
WLS Design of Sparse FIR Digital Filters Aimin Jiang, Member, IEEE, and Hon Keung Kwan, Senior Member, IEEE
Abstract—In this paper, we propose a novel algorithm for sparse finite impulse response (FIR) filter designs. The objective of the sparse digital filter design problem considered in this paper is to reduce the number of nonzero-valued filter coefficients, subject to a weighted least-squares (WLS) approximation error constraint imposed on the frequency domain. The proposed design method is inspired by the iterative shrinkage/thresholding (IST) algorithms, which are used in sparse and redundant representation for signals. The basic idea of the proposed design algorithm is to successively transform the original nonconvex problem to a series of constrained subproblems in a simpler form. Despite of their nonconvexity, these subproblems can be efficiently and reliably solved in each iterative step by a numerical approach developed in this paper. Furthermore, it can be demonstrated that the obtained solutions are essentially optimal to their respective subproblems. Since its major part only involves scalar operations, the proposed algorithm is computationally efficient. Three sets of numerical examples are presented in this paper to illustrate the effectiveness of the proposed design algorithm. Index Terms—Finite impulse response (FIR) digital filter, iterative shrinkage/thresholding (IST), sparse filter design, weighted least-squares (WLS).
I. INTRODUCTION
D
GITAL finite impulse response filters (FIR) are widely used in various applications of communication and signal processing. Classical FIR filter design focuses on developing numerical methods, which can achieve optimal solutions under some criterion. However, the implementation of obtained FIR filters may need a large number of multipliers, especially when design specifications are extremely stringent. One of potential ways to efficiently realize FIR filters is to reduce the number of arithmetic operations, which can be taken into account during the stage of designing filters. Many researchers prefer to prefilter-equalizer design approaches [1]–[4], where a prefilter is designed first, which is generally composed of one or more cascaded FIR filters with simple structures, and an equalizer is cascaded after the prefilter to further improve the performance of Manuscript received October 14, 2011; revised February 24, 2012 and April 06, 2012; accepted April 27, 2012. Date of publication September 19, 2012; date of current version January 04, 2013. This work was supported in part by the Fundamental Research Funds for the Central Universities of China under Grant 2011B11214, and the National Nature Science Foundation of China under Grant 61101158. This paper was recommended by Associate Editor C.-C. Tseng. A. Jiang is with the College of Computer and Information Engineering, Hohai University, Changzhou Campus, Changzhou 213022, China (e-mail: [email protected]). H. K. Kwan is with the Department of Electrical and Computer Engineering, University of Windsor, Windsor, ON N9B 3P4, Canada (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TCSI.2012.2215742
the whole system. In [1] and [2], prefilters are implemented by using the first 104 cyclotomic polynomials, whose coefficients are only 0 or 1, such that only equalizers require multipliers. In [3] and [4], the efficient structure of recursive running sum (RSS), which is essentially an FIR filter with zeros evenly placed on the unit circle, is utilized as a prefilter. Another way for efficient realization of FIR filters is to constrain filter coefficients in signed power-of-two (SPT) space, such that the implementation of an FIR filter only needs shifters and adders. In [5]–[8], tree search techniques are combined with other optimization approaches to locate optimal solutions. Some other algorithms [9], [10] apply convex relaxation techniques to design FIR filters with SPT coefficients, where some nonconvex constraints are relaxed such that original nonconvex design problems can be recast as convex optimization problems. The obvious advantage of this design strategy is that many efficient numerical methods (e.g., interior-point methods [11]) can be deployed to achieve approximate solutions, which may provide useful information about global solutions [12], [13]. In [14], the FIR filter design problem using SPT coefficients is directly formulated as a polynomial programming problem, which is then resolved by use of the corresponding numerical solvers. For sharp FIR filter design problems, the frequency masking technique [15]–[17] can be utilized to reduce the complexity of resulting filters. In this paper, we concern sparse FIR filter designs, which aim to find FIR filters with as few nonzero-valued coefficients as possible. If a sparse FIR filter is attained, the multipliers corresponding to the zero-valued coefficients can be omitted, which consequently reduce the implementation complexity. Theoretically speaking, we can locate optimal positions of zero-valued impulse responses by exhaustive search. However, this is impossible due to its huge computational complexity when the filter order is large. In practice, it is more realistic to find locally optimal designs. The major difficulty of local search is to determine which filter coefficients could impact the performance of an FIR filter less than some others. In [18], orthogonal matching pursuit (OMP) algorithm is employed to determine sparse solutions to equiripple filter design problems, which are recast as a linear system of equations. This algorithm is belonged to greedy methods. Compared with exhaustive search, the computational complexity of the OMP design algorithm is dramatically reduced. In [19], the branch-and-bound algorithm is applied to the design of sparse FIR filters. Following the depth-first search, filter coefficients to be zeros are gradually incorporated in subproblems until the given specifications are violated. A half-band FIR filter design algorithm is presented in [20], where by introducing auxiliary variables the sparse design problem is formulated as a mixed integer linear programming (MILP) problem and then solved by branch-and-bound or branch-and-cut technique. Although the design algorithms in
1549-8328/$31.00 © 2012 IEEE
126
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
[19] and [20] employ intelligent search to locate optimal solutions, the computational complexity is still high especially when designing FIR filters with high orders. Using linear programming (LP) as building blocks, two heuristic approaches are proposed in [21] to design sparse FIR filters. The first approach iteratively thins nonzero-valued coefficients. In each iterative step, one or more coefficients are forced to zeros using either the minimum-increase or smallest-coefficient rule. The second approach is to determine the sparsest filter by solving an approximate problem with an -norm objective function. Based on the similar idea of the second approach in [21], a two-phase design algorithm is proposed in [22]. A convex design problem with an -norm regularization term is first solved. Then, the hard thresholding operation is applied to force coefficients with small magnitudes equal to 0. Finally, a post-processing is employed to further refine the design result. Similar design strategy is also applied in sparse IIR filter design []. Many of design approaches described above only consider linear-phase FIR filter designs subject to peak error constraints imposed on magnitude response. However, some of them can be readily extended to handle general FIR filter designs. In this paper, we present an iterative shrinkage/thresholding (IST)-based algorithm for sparse FIR filter designs subject to a WLS approximation error constraint. Both linear- and nonlinear-phase FIR filter designs are handled under a unified framework. The basic idea behind the proposed design algorithm is to successively transform the original nonconvex design problem to a series of subproblems in a simpler form. Although these subproblems are still nonconvex, they can be efficiently solved by a numerical approach developed in this paper. Furthermore, by means of the Lagrange dual problems of such subproblems, it will be demonstrated in Section III-B that the obtained solutions are optimal to their respective subproblems. The paper is organized as follows. The sparse FIR filter design problem is first presented in Section II. A novel iterative design method is then proposed in Section III. Several design examples are presented in Section IV. Finally, conclusions are drawn in Section V. II. PROBLEM FORMULATION Let be a given desirable frequency response and represent the union of frequency bands of interest. A dense grid of frequency points are sampled over . The original WLS design problem can be expressed by (1) where the filter coefficient vector is defined by (2) and is a given weighting function. In (2), the superscript represents the transpose operation of a matrix or vector. In this paper, we only consider design problems of FIR filters with real-valued coefficients. However, the proposed algorithm
can be readily extended to handle design problems of FIR filters with complex-valued coefficients. Using , the frequency response of an FIR filter can be computed by (3) where (4) It is known that design problem (1) is convex and can be equivalently cast as (5) where
.. .
.. .
(6)
(7)
and denotes the Euclidean norm of a vector . In (6) and (7), and are operators used to retrieve real and imaginary parts of a complex value, respectively. The optimal solution of (5) can be obtained by solving a set of linear equations or using some convex optimization solver (e.g., SeDuMi [26] or SDPT3 [27]) if additional convex constraints are incorporated. For convenience of latter discussion, we designate as the optimal objective value of the WLS design problem (5). Based on (5), the sparse FIR filter design problem under consideration is expressed by (8a) (8b) where denotes the norm of , which is equal to the number of nonzero-valued elements of . It should be mentioned here that the norm is essentially a semi-norm as it does not satisfy the property of positive scalability. In (8), is a specified upper bound of the WLS approximation error. Without loss of generalization, we assume that (9) We can achieve a sparser design by a larger , while a lower WLS approximation error by a smaller . The design problem
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
127
(8) is highly nonconvex and, in general, it is difficult to definitely obtain its optimal solutions. Note that if we appropriately choose , , and , (8) can also be applied to formulate sparse linear-phase FIR filter design problems. For instance, in order to design an th-order Class I linear-phase FIR filter, and can be chosen, respectively, as
.. .
..
.
.. .
.. .
obtained in each iterative step by a numerical method presented in Section III-B. A. Iterative Design Procedure The proposed design algorithm adopts an iterative procedure to design sparse FIR filters. The initial point can be obtained by solving the original WLS design problem (5). In the -th iterative step, we construct a constrained subproblem other than the regularized one. Its objective function is in the same form as in (8). Constraint (8b) is modified as (11)
Correspondingly, is composed by the desired magnitude responses sampled on the frequencies , . For linear-phase FIR filter designs, the WLS approximation error constraint (8b) only needs to be imposed on the magnitude response instead of the frequency response, as the group delay of a linear-phase FIR filter is constant over the whole frequency band. The nonconvexity of (8) can be overcome by replacing the norm of by its norm in the objective function. Although it is a good approximation of (8), the -norm design problem does not directly lead to a real sparse solution. In the next section, we shall introduce an efficient design algorithm to solve (8), which can yield sparse filters directly.
where is the design result in the last iterative step and the additional term is defined by
(12) In (12), is chosen such that further implies
is always convex, which (13)
One promising way to efficiently solve (8) is to mix its -norm objective function and the -norm constraint in the form
where represents the maximal eigenvalue of a symmetric matrix . It should be mentioned that under (13) the term is always nonnegative. Thus, the feasibility domain defined by (11) is contained within the one defined by (8b). In order to make the restricted feasibility domain as large as possible, should be chosen as . After some manipulations, constraint (11) can be rewritten by
(10)
(14)
III. ITERATIVE DESIGN METHOD
where the positive regularization parameter controls the tradeoff between the approximation error and the sparsity of filter coefficients. Many numerical approaches have been proposed to tackle this unconstrained optimization problem. Among them, the IST algorithms have drawn much attention due to their computational efficiency. The IST algorithms were proposed independently by several authors. The major advantage of the IST algorithms is that in each iterative step of the IST algorithms the optimization can be decomposed into a set of independent scalar optimization problems, which generally have closed-form solutions, such that large-scale problems can be efficiently resolved. Readers are referred to [24] and [25] for a more thorough review of the IST algorithms and their variants. The proposed algorithm is based on a design strategy similar to that of the IST algorithms. However, instead of minimizing the regularized objective function , the proposed algorithm directly deals with (8) so as to avoid the dilemma of choosing the regularization parameter . The proposed algorithm contains an iterative design procedure, which is developed in Section III-A. In each iteration, a subproblem is constructed using the solution obtained in the previous iterative step. It will be shown that although this subproblem is still nonconvex, its globally optimal solution can be efficiently
where (15)
(16) In (16), represents an identity matrix whose size can be determined in context. Using (11) or, equivalently, (14), we can reformulate the subproblem in each iterative step as (17a) (17b) Note that (17) cannot be further decomposed to a set of scalar optimization problems. We shall introduce an efficient design method in the next subsection to solve (17). Let denote an optimal solution to (17). The following proposition shows that by appropriately selecting an initial point for the iterative design procedure, the feasibility domain of (17) is nonempty. Then, if in each iteration a global solution to (17)
128
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
can be always achieved, we can find a local solution of the original design problem (8). Proposition 1: If
where is obtained by extracting from the columns corresponding to the complement set of . Note that if is an empty set, (22) is also the optimal solution to (5).
(18) B. Successive Activation Algorithm the feasibility domain defined by (14) is nonempty for any . Proof: First of all, we show that for some specific , if is satisfied, is nonnegative. Since (19) and
we have
The remaining task is to solve the subproblem (17). Although (17) is still nonconvex due to the existence of the -norm objective function, compared with (8), it has been simplified. Constraint (17b) defines a sphere of radius and center on . If the origin lies inside the ball, which means , the optimal solution to (17) is a zero vector . In other cases, we employ an iterative algorithm to search the optimal solution to (17). The searching procedure starts from the origin of , and activates elements of one by one until satisfies constraint (17b). We refer to this numerical approach as the successive activation algorithm, whose main steps are shown below. Successive Activation Algorithm:
and, consequently
If , there exists at least one feasible point (e.g., ) for the design problem (17). Combined with the fact that the optimal solution of (17) also satisfies (19), in all the successive iterative steps, the feasibility domain defined by (14) is always nonempty. If the iterative procedure starts from a given which satisfies (18), then the feasibility domain defined by (14) is nonempty for any . In our designs, the iterative design procedure continues until the following condition is satisfied (20)
Inputs:
and
Output:
.
.
. If Step 1: Choose ; Otherwise, go to Step 2;
, return
Step 2: According to their magnitudes, rearrange elements of in a nonascending order, i.e., , and set
; , choose
Step 3: If
(23) as any value
and go to Step 4; Otherwise, specify
where
within the interval or the number of iterations is larger than a predefined number MaxIter. Although the convergence of the proposed design algorithm has not been strictly guaranteed, in all the designs we have conducted, the iterative procedure can always converge to final solutions. After the iterative procedure described above, we may further reduce the WLS approximation error by solving a WLS design problem similar to (5) (21a)
and return Step 4: Set
;
and go to Step 3.
Several remarks regarding the successive activation algorithm are made below. 1) In Step 3, if , is the last element of situation, we can verify that
to be activated. In this satisfies (17b) if
(21b) In (21b), is a subset of indices at which is identified to be zero by the iterative procedure described before. The nonzerovalued coefficients of the optimal solution to (21) can be computed by
, the equality of (17b) is When achieved. In our designs, we always adopt
(22)
(24)
along with
and
.
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
129
to determine the last nonzero-valued coefficient. In (24), is defined by if , otherwise.
(25)
Among all the feasible values, (24) represents the one with the smallest magnitude. This updating strategy is inspired by the -norm objective function of (17). 2) In the successive activation algorithm, most of nonzerovalued coefficients are set equal to the corresponding elements of . When the indices of zero-valued coefficients are fixed, (15) suggests that the successive activation algorithm works like the classical steepest descent method. 3) In Step 3, if several elements of have the same value ), we can proceed (e.g., with any one of them. In our designs, we always choose the one which has the smallest index number in . 4) In Step 4, if is larger than , the successive activation algorithm fails to find a solution feasible to (17). However, this could happen only when . In view of Proposition 1, it can be concluded that successive activation algorithm can converge to final solutions of (17). The successive activation algorithm can be derived from another perspective. The Lagrangian associated with (17) is defined by
one extra coefficient will be turned to a nonzero value. Furthermore, the dual objective function increases or decreases monotonically with respect to within the interval , , where denotes the subset conif sisting of indices of all the zero-valued coefficients when applying to (27). Since the objective function of the Lagrange dual problem (28) is always a concave function, its optimal solution can be found among s. Based on this fact, the searching procedure continues until begins to decrease, which can be written by or, equivalently,
(29)
Using the identity simplified as
(30)
(26) where is a Lagrange multiplier. Without the term independent of , we can decompose to a set of scalar functions with respect to . By minimizing such scalar functions with respect to , we obtain if otherwise.
,
, (29) can be further
Apparently, (30) is valid only when the following condition is satisfied (31)
(27)
be a subset of indices at which s are specified Let equal to 0 according to (27). Taking (27) into (26), we achieve the Lagrange dual problem associated with (17) (28a)
which is essentially the termination condition employed in Step 3 of the successive activation algorithm. Note that the searching procedure described above to solve (28) adopts the same rule to locate the indices of nonzero-valued coefficients as the successive activation algorithm. Thereby, it can find the same indices of nonzero-valued coefficients as the successive activation algorithm, if they start from the same initial point . Moreover, the termination condition (31) implies
(28b) The optimal solution to (28) can be achieved by a searching procedure similar to the successive activation algorithm described before. The idea behind this searching procedure is that by adjusting we can control the size of as well as the dual objective value of (28). It can be observed from (27) that when all the elements of are equal to 0, which essentially corresponds to the initial point of the successive activation algorithm. In order to locate an optimal value of for , we first rearrange elements of as in Step 2 of the successive activation algorithm. Then, starting from 0, we can augment in the subsequent steps to turn more elements of to nonzero values. Suppose that s are different from each in each step, other. Then, if we set
(32) In view of (31) and (32), we have (33) which means that the gap between the primal and dual objective values of (17) is less than 1. Since the primal objective function only takes nonnegative integer values, (33) further indicates that
130
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
the obtained primal solution is optimal to (17). It is worth noting that the only difference between the searching procedure mentioned above and the successive activation algorithm is on the value of the last nonzero coefficient. If (24) other than the corresponding element of is specified to the last nonzero coefficient in (27), it can be verified that the gap between and the primal objective value goes to 0. The Lagrange dual problem (28) can also be expressed as an LP problem by using the identity (34a) (34b) (34c) (34d) where represents a vector with all the entries equal to 1, and s are a set of auxiliary variables introduced to simplify the dual problem formulation. When we obtain an optimal solution , the corresponding primal solution can be achieved by taking into (27). In practice, however, we prefer the successive activation algorithm to solving (34) due to its computational efficiency. To conclude this section, we plot the flowchart of the complete design algorithm in Fig. 1. For the ease of notation, we refer to the successive activation algorithm as SAA. It can be observed from Fig. 1 that since the successive activation algorithm only involves scalar operations, the major computation is expended to solve the WLS design problems (5) and (21). Thereby, compared with some other design approaches, the proposed design algorithm is computationally efficient. In practice, the design procedure depicted in Fig. 1 can be repeated for several times to further improve the design performance. All stages start from the design output of the previous stage except the first one. This multi-stage design strategy will be used in each numerical design presented in Section IV. IV. SIMULATIONS In this section, three sets of examples are presented to demonstrate the effectiveness of our proposed design method. In our designs, used in (20) is always set equal to . The predefined maximum iteration number MaxIter is specified as 200 in each design presented in this section. Without any explicit declaration, the initial design of the iterative design procedure is chosen as the optimal solution to the WLS design problem (5). Similarly, if it is not explicitly defined in the specifications, the weighting function is always set equal to 1 over and 0 otherwise. The total number of frequency points uniformly sampled over the whole frequency band is chosen as in all designs. Since the sparsity (i.e., the number of zero-valued filter coefficients) is the principle objective of (8), it is adopted as the major measure of design results. The (weighted) and errors ( and , respectively, for short) are also measured on frequency response (FR), magnitude response (MR) and group delay (GD) of nonlinear-phase FIR filters. Since the group delay
Fig. 1. Flowchart of the complete design algorithm.
TABLE I SPECIFICATIONS OF EXAMPLE I
of a linear-phase FIR filter is constant over the whole frequency band, we only measure the (weighted) and errors on magnitude response of linear-phase FIR filters. All designs presented in this section are conducted on a desktop computer with an AMD Phenom II quad core processor of 3.00 GHz. A. Example 1 The first example is taken from [21], where the objective is to determine optimal weights for a uniform linear beamformer. In essence, this design problem can be equivalently translated to a linear-phase FIR filter design problem, whose specifications are given in Table I. For comparison, we also employ the successive thinning and minimum 1-norm algorithms of [21] to design FIR filters. Although these two algorithms are originally proposed for sparse FIR filter designs with peak error constraints, they can be readily extended to deal with WLS designs. Designs with various within [1, 5] are conducted in this example. Note that when the obtained solution is optimal to the WLS design problem (5). All the design results are summarized in Table II. For , magnitude and impulse responses of the linear-phase FIR filters designed, respectively, by the proposed algorithm, and the successive thinning and minimum 1-norm algorithms of [21] are depicted in Fig. 2. Due to the symmetric
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
131
TABLE II NUMBER OF ZERO-VALUED COEFFICIENTS AND APPROXIMATION ERRORS FOR DIFFERENT
structure of a linear-phase FIR filter, only half of impulse responses are shown in Fig. 2(c). In each design, the minimum-increase and smallest-coefficient selection rules of the successive thinning algorithm yield the same design result. Thereby, we only plot in Fig. 2 the magnitude and impulse responses of the FIR filter obtained by the minimum-increase rule of the successive thinning algorithm. It can be observed that the proposed algorithm can achieve design results comparable to or better than those of the successive thinning algorithm in most of designs. The proposed algorithm attains sparser designs than the successive thinning algorithm when , 4.0, and 5.0 at a cost of higher magnitude approximation errors. For and 3.5, FIR filters designed by these two algorithms have the same number of zero-valued coefficients. However, the design results obtained by the proposed algorithm have lower peak errors but higher approximation errors. In the designs of , 2.5 and 3.5, both the proposed algorithm and the successive thinning algorithm can achieve the same results. Compared with the minimm 1-norm algorithm, the proposed algorithm yields much better designs in terms of the sparsity of obtained FIR filters. For the design with , we also plot in Fig. 3 the variation of with respect to the iteration number . It can be observed that the proposed algorithm depicted in Fig. 1 is repeated for two times. During the second stage, the curve decreases monotonically with respect to . Nevertheless, several leaps appear in the curve of the first stage. A large number of experiments reveal that these leaps generally happen with the change of the number or indices of zero-valued coefficients. B. Example 2 In this example, the proposed algorithm is applied to design variable fractional delay (VFD) FIR filters. The transfer function of a VFD FIR filter is defined by
where denotes a fractional delay and polynomial of degree
is a
IN
EXAMPLE 1
By adjusting , we can control the group delay of the VFD FIR filter. The objective of VFD FIR filter designs is to find a set of polynomial coefficients , such that the frequency response of the resulting VFD filter can maximally approximate the ideal frequency response , where is a fixed integer delay. In [28], the VFD FIR filter design problem in the WLS sense is to minimize the objective function
where and are nonnegative weighting functions, is the cutoff frequency close to , is composed by , and . The optimal solution to the design problem above has been given in [28]. Since is a convex quadratic function with respect to , it can also be cast in the form of the objective function of (5), where the filter coefficient vector is obtained by stacking columns of on each other. To find sparse VFD FIR filters, we first employ the WLS design algorithm of [28] to achieve a benchmark VFD filter with , , , for , and
Then, the WLS approximation error of the benchmark VFD FIR filter is chosen as the upper bound of (8b). Finally, we apply the proposed algorithm to solve sparse VFD FIR filter design problems using the same set of specifications except and/or . The sparsity and approximation errors of VFD filters obtained by the proposed algorithm are given in Table III. Note that when and , the proposed algorithm also leads to the benchmark filter designed by the algorithm of [28]. The magnitude response and group delays of the VFD FIR filter designed by the proposed algorithm with and are depicted in Fig. 4. The variation of obtained in the same design is shown in Fig. 5. Note that although the design procedure illustrated in Fig. 1 is invoked for two times, the design procedure of the second stage converges
132
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
Fig. 3. Variation of number obtained from the design with
with respect to the iteration in Example 1.
Simulation results show that by increasing and/or we can achieve VFD FIR filters with a large number of zero-valued polynomial coefficients. Large and/or mean that the resulting VFD FIR filter has a large number of polynomial coefficients. However, this does not imply that the implementation complexity increases accordingly. For example, the benchmark filter has totally 544 polynomial coefficients, while the VFD FIR filter designed by the proposed algorithm with and has 8 extra coefficients. However, by slightly increasing the VFD FIR filter order, we can gain 114 zero-valued polynomial coefficients under the same approximation error constraint. The arithmetic operations or circuit components corresponding to these coefficients can be eliminated or deactivated. C. Example 3
Fig. 2. Magnitude and impulse responses of linear-phase FIR filters obtained by the proposed algorithm, the successive thinning algorithm [21], and the minimum 1-norm algorithm [21] in Example 1. (a) Magnitude response; (b) magnitude response in passband; (c) impulse response.
in two iterations. Thereby, we only plot in Fig. 5 the variation of obtained in the first stage.
In this example, we employ the proposed algorithm to design a linear-phase bandpass FIR filter. The detailed specifications presented in Table IV are similar to those given in Example 3 of [22]. We also adopt the design algorithm proposed in [22] to design FIR filters under the same set of specifications. The design algorithm of [22] is based on the -norm optimization. The WLS approximation error of an FIR filter is mixed with the norm of the filter coefficient vector, i.e., replacing the norm of by its norm in (10). In this way, the original nonconvex optimization problem can be reformulated in a convex form. In the subsequent discussion, we refer to the design algorithm of [22] as the hard thresholding algorithm. In our designs, we use a prescribed sparsity level (given in Table IV) on the prototype filter to achieve a sparse design, which is equivalent to using a thresholding parameter to guide the hard thresholding operation. For a fair comparison, the squared WLS approximation error of each FIR filter obtained by the hard thresholding algorithm is adopted as in (9) by the proposed algorithm. In [22], the regularization parameter has to be pre-specified by designers. Since there is no systematical way to choose , we conduct four sets of designs with different . The sparsity and approximation errors of each FIR filter are summarized in Table V for comparison. The magnitude responses of FIR filters obtained in the design with
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
133
TABLE III NUMBER OF ZERO-VALUED COEFFICIENTS AND APPROXIMATION ERRORS IN EXAMPLE 2
Fig. 5. Variation of number obtained from the design with
with respect to the iteration and in Example 2.
TABLE IV SPECIFICATIONS OF EXAMPLE 3
Fig. 4. Magnitude response and group delay of the VFD FIR filters designed by and in Example 2. (a) Magnitude the proposed algorithm with response; (b) group delay in passband.
and 20 zero-valued coefficients specified for the hard thresholding algorithm are plotted in Fig. 6. The variation of with respect to the iteration number , which is obtained in the same design, is shown in Fig. 7.
Simulation results presented in Table V show that the proposed algorithm outperforms the hard thresholding algorithm [22] in each design. Since the norm is a convex approximation of the norm, for a large , the performance of the hard thresholding algorithm is close to that of the proposed algorithm. However, when takes a small value, the prototype filter obtained by the hard thresholding algorithm is close to the optimal solution to (5). In this case, the proposed algorithm can achieve much better designs than the hard thresholding algorithm, as the original WLS design problem (5) does not take the sparsity requirement into account. V. CONCLUSIONS A novel sparse FIR filter design algorithm in the WLS sense has been presented in this paper. In order to tackle the nonconvexity of the original design problem incurred by the -norm
134
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
TABLE V NUMBER OF ZERO-VALUED COEFFICIENTS AND APPROXIMATION ERROR FOR DIFFERENT SPARSITY LEVEL IN EXAMPLE 3
Fig. 7. Variation of with respect to the iteration and 20 zero-valued coeffinumber obtained from the design with cients set for the hard thresholding algorithm in Example 3.
Fig. 6. Magnitude responses of linear-phase FIR filters obtained by the proposed algorithm and the hard thresholding algorithm [22] in Example 3. (a) Magnitude response; (b) magnitude response in passband.
greatly simplified and can be solved by the successive activation algorithm developed in Section III-B. The successive activation algorithm is of computational efficiency as it only involves scalar operations. Furthermore, it is demonstrated that the design result obtained by the successive activation algorithm is optimal to each subproblem. Design results obtained by the iterative procedure is finally refined by solving (22). The design procedure can be repeated for several times to further improve the design performance. It is worth noting that the proposed algorithm is essentially independent on the problem formulation. It can be applied to solve any design problem, which can be formulated in the form of (8). Simulation results reveal that the proposed algorithm can achieve design results comparable to or better than those of some other popular sparse filter design approaches. REFERENCES
objective function, an iterative design procedure is developed in Section III. In each iterative step, by use of the design result attained in the last iteration, a constrained subproblem is first constructed. Although it is still nonconvex, this subproblem is
[1] W. J. Oh and Y. H. Lee, “Design of efficient FIR filters with cyclotomic polynomial prefilters using mixed integer linear programming,” IEEE Signal Process. Lett., vol. 3, pp. 239–241, Aug. 1996. [2] R. J. Hartnett and G. F. Boudreaux-Bartels, “On the use of cyclotomic polynomial prefilters for efficient FIR filter design,” IEEE Trans. Signal Process., vol. 41, pp. 1766–1779, May 1993.
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
[3] J. Adams and A. Willson, Jr., “A new approach to FIR digital filters with fewer multipliers and reduced sensitivity,” IEEE Trans. Circuits Syst., vol. CAS-30, pp. 277–283, May 1983. [4] J. Adams and A. Willson, Jr., “Some efficient digital prefilter structures,” IEEE Trans. Circuits Syst., vol. CAS-31, pp. 260–266, Mar. 1984. [5] Y. C. Lim and S. R. Parker, “FIR filter design over a discrete power-of-two coefficient space,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-31, pp. 583–591, Jun. 1983. [6] Y. C. Lim and S. R. Parker, “Discrete coefficient FIR digital filter design based upon an LMS criterion,” IEEE Trans. Circuits Syst., vol. CAS-30, pp. 723–739, Oct. 1983. [7] J.-H. Lee, C.-K. Chen, and Y. C. Lim, “Design of discrete coefficient FIR digital filters with arbitrary amplitude and phase responses,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 40, pp. 444–448, Jul. 1993. [8] J.-J. Shyu and Y.-C. Lin, “A new approach to the design of discrete coefficient FIR digital filters,” IEEE Trans. Signal Process., vol. 43, pp. 310–314, Jan. 1995. [9] W.-S. Lu, “Design of FIR filters with discrete coefficients: A semidefinite programming relaxation approach,” in Proc. IEEE Int. Symp. Circuits Syst., Sydney, Australia, May 2001, pp. 297–300. [10] W.-S. Lu, “Design of FIR digital filters with discrete coefficients via convex relaxation,” in Proc. IEEE Int. Symp. Circuits Syst., Kobe, Japan, May 2005, pp. 1831–1834. [11] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge, U.K.: Cambridge Univ. Press, 2004. [12] A. Jiang and H. K. Kwan, “Minimax design of IIR digital filters using iterative SOCP,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 57, no. 6, pp. 1326–1337, Jun. 2010. [13] A. Jiang and H. K. Kwan, “Minimax design of IIR digital filters using SDP relaxation technique,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 57, pp. 378–390, Feb. 2010. [14] W.-S. Lu and T. Hinamoto, “Design of FIR digital filters with discrete coefficients via polynomial programming: Towards the global solution,” in Proc. IEEE Int. Symp. Circuits Syst., May 2007, pp. 2048–2051. [15] Y. C. Lim, “Frequency-response masking approach for the synthesis of sharp linear phase digital filters,” IEEE Trans. Circuits Syst., vol. CAS-33, no. 4, pp. 357–364, Apr. 1986. [16] Y. C. Lim and R. Lian, “Frequency-response masking approach for digital filter design: Complexity reduction via masking filter factorization,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 41, no. 8, pp. 518–525, Aug. 1994. [17] W.-S. Lu and T. Hinamoto, “Optimal design of frequency-responsemasking filters using semidefinite programming,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 50, no. 4, pp. 557–568, Apr. 2003. [18] D. Mattera, F. Palmieri, and S. Haykin, “Efficient sparse FIR filter design,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2002, vol. 2, pp. 1537–1540. [19] Y.-S. Song and Y. H. Lee, “Design of sparse FIR filters based on branch-and-bound algorithm,” in Proc. 40th Midwest Symp. Circuits Syst., Aug. 1997, vol. 2, pp. 1445–1448. [20] O. Gustafsson, L. S. DeBrunner, V. DeBrunner, and H. Johansson, “On the design of sparse half-band like FIR filters,” in Proc. 41st Asilomar Conf. Signal, Syst., Comp., Nov. 2007, pp. 1098–1102. [21] T. Baran, D. Wei, and A. V. Oppenheim, “Linear programming algorithms for sparse filter design,” IEEE Trans. Signal Process., vol. 58, pp. 1605–1617, Mar. 2010. [22] W.-S. Lu and T. Hinamto, “Digital filters with sparse coefficients,” in Proc. IEEE Int. Symp. Circuits Syst., Paris, France, May 2010, pp. 169–172.
135
[23] W.-S. Lu and T. Hinamto, “Minimax design of stable IIR filters with sparse coefficients,” in Proc. IEEE Int. Symp. Circuits Syst., Rio de Janeiro, Brazil, May 2011, pp. 398–401. [24] M. Zibulevsky and M. Elad, “L1–L2 optimization in signal and image processing,” IEEE Signal Process. Mag., vol. 27, pp. 76–88, May 2010. [25] M. Elad, M. A. T. Figueiredo, and Y. Ma, “On the role of sparse and redundant representations in image processing,” Pro. IEEE, vol. 98, pp. 972–982, Jun. 2010. [26] J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones,” Optim. Methods Softw., vol. 11/12, pp. 625–653, 1999. [27] K. C. Toh, M. J. Todd, and R. H. Tutuncu, “SDPT3-a MATLAB software package for semidefinite programming,” Optim. Methods Softw., vol. 11, pp. 545–581, 1999. [28] W.-S. Lu and T.-B. Deng, “An improved weighted least-squares design for variable fractional delay FIR filters,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 46, no. 8, pp. 1035–1040, Aug. 1999.
Aimin Jiang (S’07–M’11) received the B.E. and M.E. degrees from Nanjing University of Aeronautics and Astronautics, Nanjing, China, in 2001 and 2004, respectively, and the Ph.D. degree from University of Windsor, Windsor, ON, Canada, in 2010, all in electrical engineering. He was an Associate Professor since 2011 and was promoted to a Professor in June, 2012 at the College of Computer and Information Engineering, Hohai University, Changzhou Campus, Changzhou. His research interests include mathematical optimization and its applications to digital signal processing and communications.
Hon Keung Kwan (M’81–SM’86) received the D.I.C. and Ph.D. degree in electrical engineering (signal processing) from Imperial College London, U.K., in 1981. His previous experiences include working as a Design Engineer in the electronics and computer memory industry during 1977–1978 and serving as a faculty member in the Department of Electronic and Information Engineering in The Hong Kong Polytechnic University during 1981, and then in the Department of Electrical and Electronic Engineering in The University of Hong Kong. He subsequently joined the University of Windsor and has been a Professor in Electrical and Computer Engineering since 1989. His research interests are in the area of neural systems, digital filters, and their applications to intelligent signal processing. He has a variety of teaching, administrative, and professional experiences. Dr. Kwan had served as the Chair and various officers in each of the Digital Signal Processing Technical Committee and the Neural Systems and Applications Technical Committee of the IEEE Circuits and Systems Society. Dr. Kwan is a licensed Professional Engineer (ON, Canada), a Chartered Electrical Engineer (U.K.), and a Fellow of the Institution of Engineering and Technology (U.K.).
125
WLS Design of Sparse FIR Digital Filters Aimin Jiang, Member, IEEE, and Hon Keung Kwan, Senior Member, IEEE
Abstract—In this paper, we propose a novel algorithm for sparse finite impulse response (FIR) filter designs. The objective of the sparse digital filter design problem considered in this paper is to reduce the number of nonzero-valued filter coefficients, subject to a weighted least-squares (WLS) approximation error constraint imposed on the frequency domain. The proposed design method is inspired by the iterative shrinkage/thresholding (IST) algorithms, which are used in sparse and redundant representation for signals. The basic idea of the proposed design algorithm is to successively transform the original nonconvex problem to a series of constrained subproblems in a simpler form. Despite of their nonconvexity, these subproblems can be efficiently and reliably solved in each iterative step by a numerical approach developed in this paper. Furthermore, it can be demonstrated that the obtained solutions are essentially optimal to their respective subproblems. Since its major part only involves scalar operations, the proposed algorithm is computationally efficient. Three sets of numerical examples are presented in this paper to illustrate the effectiveness of the proposed design algorithm. Index Terms—Finite impulse response (FIR) digital filter, iterative shrinkage/thresholding (IST), sparse filter design, weighted least-squares (WLS).
I. INTRODUCTION
D
GITAL finite impulse response filters (FIR) are widely used in various applications of communication and signal processing. Classical FIR filter design focuses on developing numerical methods, which can achieve optimal solutions under some criterion. However, the implementation of obtained FIR filters may need a large number of multipliers, especially when design specifications are extremely stringent. One of potential ways to efficiently realize FIR filters is to reduce the number of arithmetic operations, which can be taken into account during the stage of designing filters. Many researchers prefer to prefilter-equalizer design approaches [1]–[4], where a prefilter is designed first, which is generally composed of one or more cascaded FIR filters with simple structures, and an equalizer is cascaded after the prefilter to further improve the performance of Manuscript received October 14, 2011; revised February 24, 2012 and April 06, 2012; accepted April 27, 2012. Date of publication September 19, 2012; date of current version January 04, 2013. This work was supported in part by the Fundamental Research Funds for the Central Universities of China under Grant 2011B11214, and the National Nature Science Foundation of China under Grant 61101158. This paper was recommended by Associate Editor C.-C. Tseng. A. Jiang is with the College of Computer and Information Engineering, Hohai University, Changzhou Campus, Changzhou 213022, China (e-mail: [email protected]). H. K. Kwan is with the Department of Electrical and Computer Engineering, University of Windsor, Windsor, ON N9B 3P4, Canada (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TCSI.2012.2215742
the whole system. In [1] and [2], prefilters are implemented by using the first 104 cyclotomic polynomials, whose coefficients are only 0 or 1, such that only equalizers require multipliers. In [3] and [4], the efficient structure of recursive running sum (RSS), which is essentially an FIR filter with zeros evenly placed on the unit circle, is utilized as a prefilter. Another way for efficient realization of FIR filters is to constrain filter coefficients in signed power-of-two (SPT) space, such that the implementation of an FIR filter only needs shifters and adders. In [5]–[8], tree search techniques are combined with other optimization approaches to locate optimal solutions. Some other algorithms [9], [10] apply convex relaxation techniques to design FIR filters with SPT coefficients, where some nonconvex constraints are relaxed such that original nonconvex design problems can be recast as convex optimization problems. The obvious advantage of this design strategy is that many efficient numerical methods (e.g., interior-point methods [11]) can be deployed to achieve approximate solutions, which may provide useful information about global solutions [12], [13]. In [14], the FIR filter design problem using SPT coefficients is directly formulated as a polynomial programming problem, which is then resolved by use of the corresponding numerical solvers. For sharp FIR filter design problems, the frequency masking technique [15]–[17] can be utilized to reduce the complexity of resulting filters. In this paper, we concern sparse FIR filter designs, which aim to find FIR filters with as few nonzero-valued coefficients as possible. If a sparse FIR filter is attained, the multipliers corresponding to the zero-valued coefficients can be omitted, which consequently reduce the implementation complexity. Theoretically speaking, we can locate optimal positions of zero-valued impulse responses by exhaustive search. However, this is impossible due to its huge computational complexity when the filter order is large. In practice, it is more realistic to find locally optimal designs. The major difficulty of local search is to determine which filter coefficients could impact the performance of an FIR filter less than some others. In [18], orthogonal matching pursuit (OMP) algorithm is employed to determine sparse solutions to equiripple filter design problems, which are recast as a linear system of equations. This algorithm is belonged to greedy methods. Compared with exhaustive search, the computational complexity of the OMP design algorithm is dramatically reduced. In [19], the branch-and-bound algorithm is applied to the design of sparse FIR filters. Following the depth-first search, filter coefficients to be zeros are gradually incorporated in subproblems until the given specifications are violated. A half-band FIR filter design algorithm is presented in [20], where by introducing auxiliary variables the sparse design problem is formulated as a mixed integer linear programming (MILP) problem and then solved by branch-and-bound or branch-and-cut technique. Although the design algorithms in
1549-8328/$31.00 © 2012 IEEE
126
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
[19] and [20] employ intelligent search to locate optimal solutions, the computational complexity is still high especially when designing FIR filters with high orders. Using linear programming (LP) as building blocks, two heuristic approaches are proposed in [21] to design sparse FIR filters. The first approach iteratively thins nonzero-valued coefficients. In each iterative step, one or more coefficients are forced to zeros using either the minimum-increase or smallest-coefficient rule. The second approach is to determine the sparsest filter by solving an approximate problem with an -norm objective function. Based on the similar idea of the second approach in [21], a two-phase design algorithm is proposed in [22]. A convex design problem with an -norm regularization term is first solved. Then, the hard thresholding operation is applied to force coefficients with small magnitudes equal to 0. Finally, a post-processing is employed to further refine the design result. Similar design strategy is also applied in sparse IIR filter design []. Many of design approaches described above only consider linear-phase FIR filter designs subject to peak error constraints imposed on magnitude response. However, some of them can be readily extended to handle general FIR filter designs. In this paper, we present an iterative shrinkage/thresholding (IST)-based algorithm for sparse FIR filter designs subject to a WLS approximation error constraint. Both linear- and nonlinear-phase FIR filter designs are handled under a unified framework. The basic idea behind the proposed design algorithm is to successively transform the original nonconvex design problem to a series of subproblems in a simpler form. Although these subproblems are still nonconvex, they can be efficiently solved by a numerical approach developed in this paper. Furthermore, by means of the Lagrange dual problems of such subproblems, it will be demonstrated in Section III-B that the obtained solutions are optimal to their respective subproblems. The paper is organized as follows. The sparse FIR filter design problem is first presented in Section II. A novel iterative design method is then proposed in Section III. Several design examples are presented in Section IV. Finally, conclusions are drawn in Section V. II. PROBLEM FORMULATION Let be a given desirable frequency response and represent the union of frequency bands of interest. A dense grid of frequency points are sampled over . The original WLS design problem can be expressed by (1) where the filter coefficient vector is defined by (2) and is a given weighting function. In (2), the superscript represents the transpose operation of a matrix or vector. In this paper, we only consider design problems of FIR filters with real-valued coefficients. However, the proposed algorithm
can be readily extended to handle design problems of FIR filters with complex-valued coefficients. Using , the frequency response of an FIR filter can be computed by (3) where (4) It is known that design problem (1) is convex and can be equivalently cast as (5) where
.. .
.. .
(6)
(7)
and denotes the Euclidean norm of a vector . In (6) and (7), and are operators used to retrieve real and imaginary parts of a complex value, respectively. The optimal solution of (5) can be obtained by solving a set of linear equations or using some convex optimization solver (e.g., SeDuMi [26] or SDPT3 [27]) if additional convex constraints are incorporated. For convenience of latter discussion, we designate as the optimal objective value of the WLS design problem (5). Based on (5), the sparse FIR filter design problem under consideration is expressed by (8a) (8b) where denotes the norm of , which is equal to the number of nonzero-valued elements of . It should be mentioned here that the norm is essentially a semi-norm as it does not satisfy the property of positive scalability. In (8), is a specified upper bound of the WLS approximation error. Without loss of generalization, we assume that (9) We can achieve a sparser design by a larger , while a lower WLS approximation error by a smaller . The design problem
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
127
(8) is highly nonconvex and, in general, it is difficult to definitely obtain its optimal solutions. Note that if we appropriately choose , , and , (8) can also be applied to formulate sparse linear-phase FIR filter design problems. For instance, in order to design an th-order Class I linear-phase FIR filter, and can be chosen, respectively, as
.. .
..
.
.. .
.. .
obtained in each iterative step by a numerical method presented in Section III-B. A. Iterative Design Procedure The proposed design algorithm adopts an iterative procedure to design sparse FIR filters. The initial point can be obtained by solving the original WLS design problem (5). In the -th iterative step, we construct a constrained subproblem other than the regularized one. Its objective function is in the same form as in (8). Constraint (8b) is modified as (11)
Correspondingly, is composed by the desired magnitude responses sampled on the frequencies , . For linear-phase FIR filter designs, the WLS approximation error constraint (8b) only needs to be imposed on the magnitude response instead of the frequency response, as the group delay of a linear-phase FIR filter is constant over the whole frequency band. The nonconvexity of (8) can be overcome by replacing the norm of by its norm in the objective function. Although it is a good approximation of (8), the -norm design problem does not directly lead to a real sparse solution. In the next section, we shall introduce an efficient design algorithm to solve (8), which can yield sparse filters directly.
where is the design result in the last iterative step and the additional term is defined by
(12) In (12), is chosen such that further implies
is always convex, which (13)
One promising way to efficiently solve (8) is to mix its -norm objective function and the -norm constraint in the form
where represents the maximal eigenvalue of a symmetric matrix . It should be mentioned that under (13) the term is always nonnegative. Thus, the feasibility domain defined by (11) is contained within the one defined by (8b). In order to make the restricted feasibility domain as large as possible, should be chosen as . After some manipulations, constraint (11) can be rewritten by
(10)
(14)
III. ITERATIVE DESIGN METHOD
where the positive regularization parameter controls the tradeoff between the approximation error and the sparsity of filter coefficients. Many numerical approaches have been proposed to tackle this unconstrained optimization problem. Among them, the IST algorithms have drawn much attention due to their computational efficiency. The IST algorithms were proposed independently by several authors. The major advantage of the IST algorithms is that in each iterative step of the IST algorithms the optimization can be decomposed into a set of independent scalar optimization problems, which generally have closed-form solutions, such that large-scale problems can be efficiently resolved. Readers are referred to [24] and [25] for a more thorough review of the IST algorithms and their variants. The proposed algorithm is based on a design strategy similar to that of the IST algorithms. However, instead of minimizing the regularized objective function , the proposed algorithm directly deals with (8) so as to avoid the dilemma of choosing the regularization parameter . The proposed algorithm contains an iterative design procedure, which is developed in Section III-A. In each iteration, a subproblem is constructed using the solution obtained in the previous iterative step. It will be shown that although this subproblem is still nonconvex, its globally optimal solution can be efficiently
where (15)
(16) In (16), represents an identity matrix whose size can be determined in context. Using (11) or, equivalently, (14), we can reformulate the subproblem in each iterative step as (17a) (17b) Note that (17) cannot be further decomposed to a set of scalar optimization problems. We shall introduce an efficient design method in the next subsection to solve (17). Let denote an optimal solution to (17). The following proposition shows that by appropriately selecting an initial point for the iterative design procedure, the feasibility domain of (17) is nonempty. Then, if in each iteration a global solution to (17)
128
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
can be always achieved, we can find a local solution of the original design problem (8). Proposition 1: If
where is obtained by extracting from the columns corresponding to the complement set of . Note that if is an empty set, (22) is also the optimal solution to (5).
(18) B. Successive Activation Algorithm the feasibility domain defined by (14) is nonempty for any . Proof: First of all, we show that for some specific , if is satisfied, is nonnegative. Since (19) and
we have
The remaining task is to solve the subproblem (17). Although (17) is still nonconvex due to the existence of the -norm objective function, compared with (8), it has been simplified. Constraint (17b) defines a sphere of radius and center on . If the origin lies inside the ball, which means , the optimal solution to (17) is a zero vector . In other cases, we employ an iterative algorithm to search the optimal solution to (17). The searching procedure starts from the origin of , and activates elements of one by one until satisfies constraint (17b). We refer to this numerical approach as the successive activation algorithm, whose main steps are shown below. Successive Activation Algorithm:
and, consequently
If , there exists at least one feasible point (e.g., ) for the design problem (17). Combined with the fact that the optimal solution of (17) also satisfies (19), in all the successive iterative steps, the feasibility domain defined by (14) is always nonempty. If the iterative procedure starts from a given which satisfies (18), then the feasibility domain defined by (14) is nonempty for any . In our designs, the iterative design procedure continues until the following condition is satisfied (20)
Inputs:
and
Output:
.
.
. If Step 1: Choose ; Otherwise, go to Step 2;
, return
Step 2: According to their magnitudes, rearrange elements of in a nonascending order, i.e., , and set
; , choose
Step 3: If
(23) as any value
and go to Step 4; Otherwise, specify
where
within the interval or the number of iterations is larger than a predefined number MaxIter. Although the convergence of the proposed design algorithm has not been strictly guaranteed, in all the designs we have conducted, the iterative procedure can always converge to final solutions. After the iterative procedure described above, we may further reduce the WLS approximation error by solving a WLS design problem similar to (5) (21a)
and return Step 4: Set
;
and go to Step 3.
Several remarks regarding the successive activation algorithm are made below. 1) In Step 3, if , is the last element of situation, we can verify that
to be activated. In this satisfies (17b) if
(21b) In (21b), is a subset of indices at which is identified to be zero by the iterative procedure described before. The nonzerovalued coefficients of the optimal solution to (21) can be computed by
, the equality of (17b) is When achieved. In our designs, we always adopt
(22)
(24)
along with
and
.
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
129
to determine the last nonzero-valued coefficient. In (24), is defined by if , otherwise.
(25)
Among all the feasible values, (24) represents the one with the smallest magnitude. This updating strategy is inspired by the -norm objective function of (17). 2) In the successive activation algorithm, most of nonzerovalued coefficients are set equal to the corresponding elements of . When the indices of zero-valued coefficients are fixed, (15) suggests that the successive activation algorithm works like the classical steepest descent method. 3) In Step 3, if several elements of have the same value ), we can proceed (e.g., with any one of them. In our designs, we always choose the one which has the smallest index number in . 4) In Step 4, if is larger than , the successive activation algorithm fails to find a solution feasible to (17). However, this could happen only when . In view of Proposition 1, it can be concluded that successive activation algorithm can converge to final solutions of (17). The successive activation algorithm can be derived from another perspective. The Lagrangian associated with (17) is defined by
one extra coefficient will be turned to a nonzero value. Furthermore, the dual objective function increases or decreases monotonically with respect to within the interval , , where denotes the subset conif sisting of indices of all the zero-valued coefficients when applying to (27). Since the objective function of the Lagrange dual problem (28) is always a concave function, its optimal solution can be found among s. Based on this fact, the searching procedure continues until begins to decrease, which can be written by or, equivalently,
(29)
Using the identity simplified as
(30)
(26) where is a Lagrange multiplier. Without the term independent of , we can decompose to a set of scalar functions with respect to . By minimizing such scalar functions with respect to , we obtain if otherwise.
,
, (29) can be further
Apparently, (30) is valid only when the following condition is satisfied (31)
(27)
be a subset of indices at which s are specified Let equal to 0 according to (27). Taking (27) into (26), we achieve the Lagrange dual problem associated with (17) (28a)
which is essentially the termination condition employed in Step 3 of the successive activation algorithm. Note that the searching procedure described above to solve (28) adopts the same rule to locate the indices of nonzero-valued coefficients as the successive activation algorithm. Thereby, it can find the same indices of nonzero-valued coefficients as the successive activation algorithm, if they start from the same initial point . Moreover, the termination condition (31) implies
(28b) The optimal solution to (28) can be achieved by a searching procedure similar to the successive activation algorithm described before. The idea behind this searching procedure is that by adjusting we can control the size of as well as the dual objective value of (28). It can be observed from (27) that when all the elements of are equal to 0, which essentially corresponds to the initial point of the successive activation algorithm. In order to locate an optimal value of for , we first rearrange elements of as in Step 2 of the successive activation algorithm. Then, starting from 0, we can augment in the subsequent steps to turn more elements of to nonzero values. Suppose that s are different from each in each step, other. Then, if we set
(32) In view of (31) and (32), we have (33) which means that the gap between the primal and dual objective values of (17) is less than 1. Since the primal objective function only takes nonnegative integer values, (33) further indicates that
130
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
the obtained primal solution is optimal to (17). It is worth noting that the only difference between the searching procedure mentioned above and the successive activation algorithm is on the value of the last nonzero coefficient. If (24) other than the corresponding element of is specified to the last nonzero coefficient in (27), it can be verified that the gap between and the primal objective value goes to 0. The Lagrange dual problem (28) can also be expressed as an LP problem by using the identity (34a) (34b) (34c) (34d) where represents a vector with all the entries equal to 1, and s are a set of auxiliary variables introduced to simplify the dual problem formulation. When we obtain an optimal solution , the corresponding primal solution can be achieved by taking into (27). In practice, however, we prefer the successive activation algorithm to solving (34) due to its computational efficiency. To conclude this section, we plot the flowchart of the complete design algorithm in Fig. 1. For the ease of notation, we refer to the successive activation algorithm as SAA. It can be observed from Fig. 1 that since the successive activation algorithm only involves scalar operations, the major computation is expended to solve the WLS design problems (5) and (21). Thereby, compared with some other design approaches, the proposed design algorithm is computationally efficient. In practice, the design procedure depicted in Fig. 1 can be repeated for several times to further improve the design performance. All stages start from the design output of the previous stage except the first one. This multi-stage design strategy will be used in each numerical design presented in Section IV. IV. SIMULATIONS In this section, three sets of examples are presented to demonstrate the effectiveness of our proposed design method. In our designs, used in (20) is always set equal to . The predefined maximum iteration number MaxIter is specified as 200 in each design presented in this section. Without any explicit declaration, the initial design of the iterative design procedure is chosen as the optimal solution to the WLS design problem (5). Similarly, if it is not explicitly defined in the specifications, the weighting function is always set equal to 1 over and 0 otherwise. The total number of frequency points uniformly sampled over the whole frequency band is chosen as in all designs. Since the sparsity (i.e., the number of zero-valued filter coefficients) is the principle objective of (8), it is adopted as the major measure of design results. The (weighted) and errors ( and , respectively, for short) are also measured on frequency response (FR), magnitude response (MR) and group delay (GD) of nonlinear-phase FIR filters. Since the group delay
Fig. 1. Flowchart of the complete design algorithm.
TABLE I SPECIFICATIONS OF EXAMPLE I
of a linear-phase FIR filter is constant over the whole frequency band, we only measure the (weighted) and errors on magnitude response of linear-phase FIR filters. All designs presented in this section are conducted on a desktop computer with an AMD Phenom II quad core processor of 3.00 GHz. A. Example 1 The first example is taken from [21], where the objective is to determine optimal weights for a uniform linear beamformer. In essence, this design problem can be equivalently translated to a linear-phase FIR filter design problem, whose specifications are given in Table I. For comparison, we also employ the successive thinning and minimum 1-norm algorithms of [21] to design FIR filters. Although these two algorithms are originally proposed for sparse FIR filter designs with peak error constraints, they can be readily extended to deal with WLS designs. Designs with various within [1, 5] are conducted in this example. Note that when the obtained solution is optimal to the WLS design problem (5). All the design results are summarized in Table II. For , magnitude and impulse responses of the linear-phase FIR filters designed, respectively, by the proposed algorithm, and the successive thinning and minimum 1-norm algorithms of [21] are depicted in Fig. 2. Due to the symmetric
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
131
TABLE II NUMBER OF ZERO-VALUED COEFFICIENTS AND APPROXIMATION ERRORS FOR DIFFERENT
structure of a linear-phase FIR filter, only half of impulse responses are shown in Fig. 2(c). In each design, the minimum-increase and smallest-coefficient selection rules of the successive thinning algorithm yield the same design result. Thereby, we only plot in Fig. 2 the magnitude and impulse responses of the FIR filter obtained by the minimum-increase rule of the successive thinning algorithm. It can be observed that the proposed algorithm can achieve design results comparable to or better than those of the successive thinning algorithm in most of designs. The proposed algorithm attains sparser designs than the successive thinning algorithm when , 4.0, and 5.0 at a cost of higher magnitude approximation errors. For and 3.5, FIR filters designed by these two algorithms have the same number of zero-valued coefficients. However, the design results obtained by the proposed algorithm have lower peak errors but higher approximation errors. In the designs of , 2.5 and 3.5, both the proposed algorithm and the successive thinning algorithm can achieve the same results. Compared with the minimm 1-norm algorithm, the proposed algorithm yields much better designs in terms of the sparsity of obtained FIR filters. For the design with , we also plot in Fig. 3 the variation of with respect to the iteration number . It can be observed that the proposed algorithm depicted in Fig. 1 is repeated for two times. During the second stage, the curve decreases monotonically with respect to . Nevertheless, several leaps appear in the curve of the first stage. A large number of experiments reveal that these leaps generally happen with the change of the number or indices of zero-valued coefficients. B. Example 2 In this example, the proposed algorithm is applied to design variable fractional delay (VFD) FIR filters. The transfer function of a VFD FIR filter is defined by
where denotes a fractional delay and polynomial of degree
is a
IN
EXAMPLE 1
By adjusting , we can control the group delay of the VFD FIR filter. The objective of VFD FIR filter designs is to find a set of polynomial coefficients , such that the frequency response of the resulting VFD filter can maximally approximate the ideal frequency response , where is a fixed integer delay. In [28], the VFD FIR filter design problem in the WLS sense is to minimize the objective function
where and are nonnegative weighting functions, is the cutoff frequency close to , is composed by , and . The optimal solution to the design problem above has been given in [28]. Since is a convex quadratic function with respect to , it can also be cast in the form of the objective function of (5), where the filter coefficient vector is obtained by stacking columns of on each other. To find sparse VFD FIR filters, we first employ the WLS design algorithm of [28] to achieve a benchmark VFD filter with , , , for , and
Then, the WLS approximation error of the benchmark VFD FIR filter is chosen as the upper bound of (8b). Finally, we apply the proposed algorithm to solve sparse VFD FIR filter design problems using the same set of specifications except and/or . The sparsity and approximation errors of VFD filters obtained by the proposed algorithm are given in Table III. Note that when and , the proposed algorithm also leads to the benchmark filter designed by the algorithm of [28]. The magnitude response and group delays of the VFD FIR filter designed by the proposed algorithm with and are depicted in Fig. 4. The variation of obtained in the same design is shown in Fig. 5. Note that although the design procedure illustrated in Fig. 1 is invoked for two times, the design procedure of the second stage converges
132
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
Fig. 3. Variation of number obtained from the design with
with respect to the iteration in Example 1.
Simulation results show that by increasing and/or we can achieve VFD FIR filters with a large number of zero-valued polynomial coefficients. Large and/or mean that the resulting VFD FIR filter has a large number of polynomial coefficients. However, this does not imply that the implementation complexity increases accordingly. For example, the benchmark filter has totally 544 polynomial coefficients, while the VFD FIR filter designed by the proposed algorithm with and has 8 extra coefficients. However, by slightly increasing the VFD FIR filter order, we can gain 114 zero-valued polynomial coefficients under the same approximation error constraint. The arithmetic operations or circuit components corresponding to these coefficients can be eliminated or deactivated. C. Example 3
Fig. 2. Magnitude and impulse responses of linear-phase FIR filters obtained by the proposed algorithm, the successive thinning algorithm [21], and the minimum 1-norm algorithm [21] in Example 1. (a) Magnitude response; (b) magnitude response in passband; (c) impulse response.
in two iterations. Thereby, we only plot in Fig. 5 the variation of obtained in the first stage.
In this example, we employ the proposed algorithm to design a linear-phase bandpass FIR filter. The detailed specifications presented in Table IV are similar to those given in Example 3 of [22]. We also adopt the design algorithm proposed in [22] to design FIR filters under the same set of specifications. The design algorithm of [22] is based on the -norm optimization. The WLS approximation error of an FIR filter is mixed with the norm of the filter coefficient vector, i.e., replacing the norm of by its norm in (10). In this way, the original nonconvex optimization problem can be reformulated in a convex form. In the subsequent discussion, we refer to the design algorithm of [22] as the hard thresholding algorithm. In our designs, we use a prescribed sparsity level (given in Table IV) on the prototype filter to achieve a sparse design, which is equivalent to using a thresholding parameter to guide the hard thresholding operation. For a fair comparison, the squared WLS approximation error of each FIR filter obtained by the hard thresholding algorithm is adopted as in (9) by the proposed algorithm. In [22], the regularization parameter has to be pre-specified by designers. Since there is no systematical way to choose , we conduct four sets of designs with different . The sparsity and approximation errors of each FIR filter are summarized in Table V for comparison. The magnitude responses of FIR filters obtained in the design with
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
133
TABLE III NUMBER OF ZERO-VALUED COEFFICIENTS AND APPROXIMATION ERRORS IN EXAMPLE 2
Fig. 5. Variation of number obtained from the design with
with respect to the iteration and in Example 2.
TABLE IV SPECIFICATIONS OF EXAMPLE 3
Fig. 4. Magnitude response and group delay of the VFD FIR filters designed by and in Example 2. (a) Magnitude the proposed algorithm with response; (b) group delay in passband.
and 20 zero-valued coefficients specified for the hard thresholding algorithm are plotted in Fig. 6. The variation of with respect to the iteration number , which is obtained in the same design, is shown in Fig. 7.
Simulation results presented in Table V show that the proposed algorithm outperforms the hard thresholding algorithm [22] in each design. Since the norm is a convex approximation of the norm, for a large , the performance of the hard thresholding algorithm is close to that of the proposed algorithm. However, when takes a small value, the prototype filter obtained by the hard thresholding algorithm is close to the optimal solution to (5). In this case, the proposed algorithm can achieve much better designs than the hard thresholding algorithm, as the original WLS design problem (5) does not take the sparsity requirement into account. V. CONCLUSIONS A novel sparse FIR filter design algorithm in the WLS sense has been presented in this paper. In order to tackle the nonconvexity of the original design problem incurred by the -norm
134
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 60, NO. 1, JANUARY 2013
TABLE V NUMBER OF ZERO-VALUED COEFFICIENTS AND APPROXIMATION ERROR FOR DIFFERENT SPARSITY LEVEL IN EXAMPLE 3
Fig. 7. Variation of with respect to the iteration and 20 zero-valued coeffinumber obtained from the design with cients set for the hard thresholding algorithm in Example 3.
Fig. 6. Magnitude responses of linear-phase FIR filters obtained by the proposed algorithm and the hard thresholding algorithm [22] in Example 3. (a) Magnitude response; (b) magnitude response in passband.
greatly simplified and can be solved by the successive activation algorithm developed in Section III-B. The successive activation algorithm is of computational efficiency as it only involves scalar operations. Furthermore, it is demonstrated that the design result obtained by the successive activation algorithm is optimal to each subproblem. Design results obtained by the iterative procedure is finally refined by solving (22). The design procedure can be repeated for several times to further improve the design performance. It is worth noting that the proposed algorithm is essentially independent on the problem formulation. It can be applied to solve any design problem, which can be formulated in the form of (8). Simulation results reveal that the proposed algorithm can achieve design results comparable to or better than those of some other popular sparse filter design approaches. REFERENCES
objective function, an iterative design procedure is developed in Section III. In each iterative step, by use of the design result attained in the last iteration, a constrained subproblem is first constructed. Although it is still nonconvex, this subproblem is
[1] W. J. Oh and Y. H. Lee, “Design of efficient FIR filters with cyclotomic polynomial prefilters using mixed integer linear programming,” IEEE Signal Process. Lett., vol. 3, pp. 239–241, Aug. 1996. [2] R. J. Hartnett and G. F. Boudreaux-Bartels, “On the use of cyclotomic polynomial prefilters for efficient FIR filter design,” IEEE Trans. Signal Process., vol. 41, pp. 1766–1779, May 1993.
JIANG AND KWAN: WLS DESIGN OF SPARSE FIR DIGITAL FILTERS
[3] J. Adams and A. Willson, Jr., “A new approach to FIR digital filters with fewer multipliers and reduced sensitivity,” IEEE Trans. Circuits Syst., vol. CAS-30, pp. 277–283, May 1983. [4] J. Adams and A. Willson, Jr., “Some efficient digital prefilter structures,” IEEE Trans. Circuits Syst., vol. CAS-31, pp. 260–266, Mar. 1984. [5] Y. C. Lim and S. R. Parker, “FIR filter design over a discrete power-of-two coefficient space,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-31, pp. 583–591, Jun. 1983. [6] Y. C. Lim and S. R. Parker, “Discrete coefficient FIR digital filter design based upon an LMS criterion,” IEEE Trans. Circuits Syst., vol. CAS-30, pp. 723–739, Oct. 1983. [7] J.-H. Lee, C.-K. Chen, and Y. C. Lim, “Design of discrete coefficient FIR digital filters with arbitrary amplitude and phase responses,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 40, pp. 444–448, Jul. 1993. [8] J.-J. Shyu and Y.-C. Lin, “A new approach to the design of discrete coefficient FIR digital filters,” IEEE Trans. Signal Process., vol. 43, pp. 310–314, Jan. 1995. [9] W.-S. Lu, “Design of FIR filters with discrete coefficients: A semidefinite programming relaxation approach,” in Proc. IEEE Int. Symp. Circuits Syst., Sydney, Australia, May 2001, pp. 297–300. [10] W.-S. Lu, “Design of FIR digital filters with discrete coefficients via convex relaxation,” in Proc. IEEE Int. Symp. Circuits Syst., Kobe, Japan, May 2005, pp. 1831–1834. [11] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge, U.K.: Cambridge Univ. Press, 2004. [12] A. Jiang and H. K. Kwan, “Minimax design of IIR digital filters using iterative SOCP,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 57, no. 6, pp. 1326–1337, Jun. 2010. [13] A. Jiang and H. K. Kwan, “Minimax design of IIR digital filters using SDP relaxation technique,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 57, pp. 378–390, Feb. 2010. [14] W.-S. Lu and T. Hinamoto, “Design of FIR digital filters with discrete coefficients via polynomial programming: Towards the global solution,” in Proc. IEEE Int. Symp. Circuits Syst., May 2007, pp. 2048–2051. [15] Y. C. Lim, “Frequency-response masking approach for the synthesis of sharp linear phase digital filters,” IEEE Trans. Circuits Syst., vol. CAS-33, no. 4, pp. 357–364, Apr. 1986. [16] Y. C. Lim and R. Lian, “Frequency-response masking approach for digital filter design: Complexity reduction via masking filter factorization,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 41, no. 8, pp. 518–525, Aug. 1994. [17] W.-S. Lu and T. Hinamoto, “Optimal design of frequency-responsemasking filters using semidefinite programming,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 50, no. 4, pp. 557–568, Apr. 2003. [18] D. Mattera, F. Palmieri, and S. Haykin, “Efficient sparse FIR filter design,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2002, vol. 2, pp. 1537–1540. [19] Y.-S. Song and Y. H. Lee, “Design of sparse FIR filters based on branch-and-bound algorithm,” in Proc. 40th Midwest Symp. Circuits Syst., Aug. 1997, vol. 2, pp. 1445–1448. [20] O. Gustafsson, L. S. DeBrunner, V. DeBrunner, and H. Johansson, “On the design of sparse half-band like FIR filters,” in Proc. 41st Asilomar Conf. Signal, Syst., Comp., Nov. 2007, pp. 1098–1102. [21] T. Baran, D. Wei, and A. V. Oppenheim, “Linear programming algorithms for sparse filter design,” IEEE Trans. Signal Process., vol. 58, pp. 1605–1617, Mar. 2010. [22] W.-S. Lu and T. Hinamto, “Digital filters with sparse coefficients,” in Proc. IEEE Int. Symp. Circuits Syst., Paris, France, May 2010, pp. 169–172.
135
[23] W.-S. Lu and T. Hinamto, “Minimax design of stable IIR filters with sparse coefficients,” in Proc. IEEE Int. Symp. Circuits Syst., Rio de Janeiro, Brazil, May 2011, pp. 398–401. [24] M. Zibulevsky and M. Elad, “L1–L2 optimization in signal and image processing,” IEEE Signal Process. Mag., vol. 27, pp. 76–88, May 2010. [25] M. Elad, M. A. T. Figueiredo, and Y. Ma, “On the role of sparse and redundant representations in image processing,” Pro. IEEE, vol. 98, pp. 972–982, Jun. 2010. [26] J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones,” Optim. Methods Softw., vol. 11/12, pp. 625–653, 1999. [27] K. C. Toh, M. J. Todd, and R. H. Tutuncu, “SDPT3-a MATLAB software package for semidefinite programming,” Optim. Methods Softw., vol. 11, pp. 545–581, 1999. [28] W.-S. Lu and T.-B. Deng, “An improved weighted least-squares design for variable fractional delay FIR filters,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 46, no. 8, pp. 1035–1040, Aug. 1999.
Aimin Jiang (S’07–M’11) received the B.E. and M.E. degrees from Nanjing University of Aeronautics and Astronautics, Nanjing, China, in 2001 and 2004, respectively, and the Ph.D. degree from University of Windsor, Windsor, ON, Canada, in 2010, all in electrical engineering. He was an Associate Professor since 2011 and was promoted to a Professor in June, 2012 at the College of Computer and Information Engineering, Hohai University, Changzhou Campus, Changzhou. His research interests include mathematical optimization and its applications to digital signal processing and communications.
Hon Keung Kwan (M’81–SM’86) received the D.I.C. and Ph.D. degree in electrical engineering (signal processing) from Imperial College London, U.K., in 1981. His previous experiences include working as a Design Engineer in the electronics and computer memory industry during 1977–1978 and serving as a faculty member in the Department of Electronic and Information Engineering in The Hong Kong Polytechnic University during 1981, and then in the Department of Electrical and Electronic Engineering in The University of Hong Kong. He subsequently joined the University of Windsor and has been a Professor in Electrical and Computer Engineering since 1989. His research interests are in the area of neural systems, digital filters, and their applications to intelligent signal processing. He has a variety of teaching, administrative, and professional experiences. Dr. Kwan had served as the Chair and various officers in each of the Digital Signal Processing Technical Committee and the Neural Systems and Applications Technical Committee of the IEEE Circuits and Systems Society. Dr. Kwan is a licensed Professional Engineer (ON, Canada), a Chartered Electrical Engineer (U.K.), and a Fellow of the Institution of Engineering and Technology (U.K.).
