# Improved Computing-Efficiency Least-Squares Algorithm with Application to All-Pass Filter Design

Authors: Lo-Chyuan Su; Yue-Dar Jou; Fu-Kun Chen
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

All-pass filter design can be generally achieved by solving a system of linear equations. The associated matrices involved in the set of linear equations can be further formulated as a Toeplitz-plus-Hankel form such that a matrix inversion is avoided. Consequently, the optimal filter coefficients can be solved by using computationally efficient Levinson algorithms or Cholesky decomposition technique. In this paper, based on trigonometric identities and sampling the frequency band of interest uniformly, the authors proposed closed-form expressions to compute the elements of the Toeplitz-plus-Hankel matrix required in the least-squares design of IIR all-pass filters. Simulation results confirm that the proposed method achieves good performance as well as effectiveness.


---

## Body

## 1. Introduction

Digital all-pass filters have received much attention in many signal processing applications such as notch filtering, phase equalization in communication systems, multichannel filter bank, construction of a wavelet, and image compression [1–9]. The most commonly used design methods can be classified into four categories: (1) maximally flat design [10], (2) least-squares approximation [11–15], (3) minimax approximation [16–18], and (4) cepstral coefficients method [19, 20]. The aforementioned approaches [11–18] include solving a set of linear equations, using a linear programming method, or utilizing a generalized exchange algorithm. Several different optimization methods have been developed based on efficiency or applications. No matter what methods are used, the computation increases significantly as the filter length increases.In the literature [13], an efficient and robust weighted least-squares (WLS) method is proposed to design IIR all-pass filters that have a least-squares or an equiripple phase error response. The method is based on formulating a weighted square error reflecting the difference between the desired phase and the phase response of the design IIR all-pass filter in a quadratic form. The filter coefficients are then obtained by solving a system of linear equations which involves a Toeplitz-plus-Hankel matrix. Consequently, the optimal solution of this system of linear equations can be obtained efficiently using Cholesky decomposition or Levinson algorithm [21–23]. These efficient algorithms require only O(N2) complexity [21–24] to compute the optimal solution. It is computationally efficient as compared to solving the linear equations by directly computing a matrix inversion which involves O(N3) complexity. Furthermore, the computation of matrix inversion may cause numerical problems as the filter length is large. Therefore, the method proposed by Kidambi [13, 24] is not only efficient but also robust as applied to the filter design problems. It has been shown in [24] that the efficient algorithms do not result in numerical problem as the filter increases to 2000.As observed from the elements of Toeplitz-plus-Hankel matrix, they are associated with the filter specifications when the design problem is specified. Therefore, trigonometric properties can be further used to reduce the computation. In this paper, based on trigonometric identities and uniformly sampling the band of interest, the elements of Toeplitz-plus-Hankel matrix can be further simplified as closed-form expressions. Simulation results indicate that the computations can be greatly reduced as compared to the efficient algorithm presented by [13].This paper is organized as follows. In Section2, the problem formulation of IIR all-pass filter is briefly reviewed. In Section 3, the closed-form expressions for Toeplitz-plus-Hankel matrix are investigated in detail. The analysis of computational requirements is stated in Section 4. The simulation and design examples are shown in Section 5 to verify the advantages of the proposed method. Finally, the conclusion is stated in Section 6.


## 2. Problem Formulation

The frequency response of an IIR all-pass filter withN real coefficients, h(n),n=0,…,N-1, can be expressed as
(1)H(ejω)=∑n=0Nh(N-n)e-jωn∑n=0Nh(n)e-jωn=e-jωND(e-jω)D(ejω),
where the denominator is an FIR filter with frequency response D(ejω)=∑n=0Nh(n)e-jωn and h(0)=1. It can be found that the magnitude response of H(ejω) is unit gain at all frequencies. The phase response of H(ejω) can be written as
(2)θ(ω)=-Nω+2tan-1(∑n=1Nh(n)sin(nω)1+∑n=1Nh(n)cos(nω)).
Given a desired phase response, θd(ω),0≤ω≤π, the purpose is to find the real coefficients h(n) such that phase errors e(ω) can be minimized in the least-squares sense. The phase error response between θ(ω) and θd(ω) can be represented as
(3)e(ω)=θd(ω)-θ(ω)=θd(ω)+Nω-2tan-1(∑n=1Nh(n)sin(nω)1+∑n=1Nh(n)cos(nω)).It is clear that the phase error response in (3) is a nonlinear representation of filter coefficients. Minimizing the phase error response directly will result in a highly nonlinear optimization problem. Based on the assumption of phase approximation presented by Kidambi [13], (2) can be alternatively rewritten as
(4)tan(θ(ω)+Nω2)=hTs(ω)1+hTc(ω),
where
(5)h=[h(1),h(2),…,h(N)]T,c(ω)=[cos(ω),cos(2ω),…,cos(Nω)]T,s(ω)=[sin(ω),sin(2ω),…,sin(Nω)]T.
If the phase response of the design IIR all-pass filter is very close to the desired phase, that is, θd(ω)≈θ(ω), then the expression in (4) can be approximated as [13]
(6)sin[ρd(ω)]cos[ρd(ω)]≈hTs(ω)1+hTc(ω),
where ρd(ω)=(θd(ω)+Nω)/2. After performing some mathematical treatment, (6) can be simplified as
(7)hT[sin[ρd(ω)]c(ω)-cos[ρd(ω)]s(ω)]≈-sin[ρd(ω)].
As a consequence, (7) can be expressed in a more compact form
(8)hTs1(ω)≈-sin[ρd(ω)],
where
(9)s1(ω)=sin[ρd(ω)]c(ω)-cos[ρd(ω)]s(ω),=[sin(ρd(ω)-ω)sin(ρd(ω)-2ω)⋯sin(ρd(ω)-Nω)]T.
Therefore, the objective function to be minimized in the least-squares sense can be formulated in a quadratic form as
(10)E=∑l=1L{hTs1(ωl)+sin(ρd(ωl))}2,
where Lis the number of frequency points at which the desired phase is sampled. By setting ∂E/∂h(i)=0, for i=1,2,…,N, a system of linear equations given by Qh=d is easily obtained, where
(11)Q=∑l=1Ls1(ωl)·s1T(ωl),d=-∑l=1Lsin[ρd(ωl)]s1(ωl).SinceQ is a real, symmetric, and positive-definite matrix, a unique solution is guaranteed. In solving the system of linear equations, a matrix inversion Q-1 which involves O(N3) complexity can be used directly. As the filter length is long, the computational complexity increases significantly. Furthermore, large filter length may cause the numerical problem when the matrix is ill conditioned. Therefore, algorithms with efficient to design and robust to avoid numerical problem should be used.


## 3. Proposed Closed-Form Expressions for Toeplitz-Plus-Hankel Matrix

Kidambi [13, 24] expanded the matrix Q into a sum of symmetric Toeplitz and Hankel matrix. Consequently, the system of linear equations can be written as (T+H)h=d. There exist several iterative and efficient algorithms [21–23] which involve O(N2) complexity that can be used to obtain the optimal filter. In this paper, some trigonometric identities [25] and uniformly sampled frequency points are exploited such that the Toeplitz and Hankel matrices can be further simplified as closed-form expressions. As a result, the computational requirements for the associated matrices can be greatly reduced especially for filter length is large.Evidently, each element of matrixQ is the sum of cosine functions and can be written as a Toeplitz-plus-Hankel matrix form:
(12)Q(i,j)=∑l=1Lsin[ρd(ωl)-iωl]·sin[ρd(ωl)-jωl]=12{∑l=1Lcos(i-j)ωl-∑l=1Lcos[(i+j)ωl-2ρd(ωl)]}=T(i,j)+H(i,j).
These elements can be further expanded as the matrix forms(13)T=[T(i,j)]=12[∑l=1L1∑l=1Lcos(ωl)⋯∑l=1Lcos(N-1)ωl∑l=1Lcos(ωl)∑l=1L1⋯∑l=1Lcos(N-2)ωl⋮⋮⋱⋮∑l=1Lcos(N-1)ωl∑l=1Lcos(N-2)ωl⋯∑l=1L1],H=[H(i,j)]=-12[∑l=1Lcos(2ωl-2ρd(ωl))⋯∑l=1Lcos((N+1)ωl-2ρd(ωl))⋮⋱⋮∑l=1Lcos((N+1)ωl-2ρd(ωl))⋯∑l=1Lcos(2Nωl-2ρd(ωl))].Therefore, only the first row (or column) of T has to be evaluated. There are N distinct elements required to be computed for an N×N Toeplitz matrix. Similarly, 2N-1 distinct elements at the first row (or column) and the last row (or column) are required to be computed for an N×N Hankel matrix. Using the Toeplitz-plus-Hankel matrix, the computations can be significantly reduced. As observed from (12), the summation is performed on the discrete frequency points, so the computations are affected by the frequency sampling points. If L points are  uniformly sampled in the band of interest [ωo,ωL], ωl can be formulated as
(14)ωl=ωo+(ωL-ωo)·lL=ωo+Δl,l=1,2,…,L,
where Δ=(ωL-ωo)/L is the spacing between two adjacent frequency points. To obtain the closed-form expressions, the following useful formulae which have been shown in [25] can be used:
(15)∑l=mncos(α+lβ)=cos[α+(m+n)β/2]sin[(n-m+1)β/2]sin(β/2),∑l=mnsin(α+lβ)=sin[α+(m+n)β/2]sin[(n-m+1)β/2]sin(β/2).
Consequently, only the closed-form expressions for Toeplitz matrix can be readily achieved for the case of i=j,
(16)T(i,j)=12∑l=1Lcos[0·ωl]=12∑l=1L1=L2.
And for the case of i≠j,
(17)T(i,j)=12∑l=1Lcos[(i-j)(ωo+lΔ)]=12∑l=1Lcos[(i-j)ωo+lΔ(i-j)]=12cos[(i-j)ωo+(L+1)(i-j)Δ/2]·sin[(i-j)LΔ/2]sin[(i-j)Δ/2].
Similarly, the closed-form expression for the elements of the Hankel matrix can be written as
(18)H(i,j)=-12∑l=1Lcos[(i+j)ωl-2ρd(ωl)].
It is clear that the Hankel matrix is not only dependent on cosine function but also dependent on the desired phase response. Considering for a Hilbert transformer with phase response θd(ω)=-Nω-π/2 is designed, it implies that ρd(ω)=-π/4. Accordingly, (18) can be further simplified as
(19)H(i,j)=-12∑l=1Lcos[(i+j)ωl+π2]=12∑l=1Lsin[(i+j)(ωo+lΔ)]=12sin[(i+j)ωo+(L+1)(i+j)Δ/2]·sin[(i+j)LΔ/2]sin[(i+j)Δ/2].
Performing the analogous procedure, the closed-form expression of d can be shown below:
(20)d(i)=-∑l=1Lsin[ρd(ωl)]·sin[ρd(ωl)-iωl]=-12∑l=1L{cos(iωl)-cos[iωl-2ρd(ωl)]}=-12∑l=1Lcos(iωl)+12∑l=1Lcos(iωl+π2)=-12∑l=1Lcos(iωl)-12∑l=1Lsin(iωl)=-12cos[iωo+(L+1)iΔ/2]sin(iLΔ/2)sin(i·Δ/2)-12sin[iωo+(L+1)iΔ/2]sin(iLΔ/2)sin(i·Δ/2)=-22sin(iLΔ/2)sin(iΔ/2)·sin[iωo+(L+1)iΔ2+π4].
Using these closed-form expressions for the associated matrices in the system of linear equations, the computational complexity can be significantly reduced.In order to verify the accuracy of the closed-form expressions for Toeplitz-plus-Hankel matrix calculated, a Hilbert transformer with lengthN=6, ωo=0.06π, ωL=0.94π, and L=60 is designed [14]. Our results agree with those obtained from the method in [13]. It can be found that the difference occurs at approximately 15 decimal places of accuracy for the computations of Toeplitz-plus-Hankel matrix and vector d. Therefore, the proposed method is not only efficient but also accurate for the computation of associated matrices. As the closed-form expressions of the associated matrices have been determined, the efficient Levinson algorithm [21–23] or Cholesky decomposition can be used to solve the optimal solution.The problem of Toeplitz-plus-Hankel system of equations can be converted to a block Toeplitz matrix using a block Levinson algorithm [21–23]. The system (T+H)h=d is equivalent to the following symmetric system:
(21)[TJHJHH][Jhh]=[Jdd],
where the N×N matrix J denotes the reverse matrix with unities on the antidiagonal and zeros elsewhere as shown below:
(22)J=[00⋯10⋮10⋮1⋱⋮10⋯0].
It is clear that the coefficient matrix in (21) is Toeplitz and symmetric. Although, the size of the problem has increased by a factor of two, the desired solution is symmetric. As a consequence, an efficient block Levinson algorithm [21, 23] can be applied to solve this system. Furthermore, the symmetric split Levinson algorithm [22] that uses only symmetric quantities with lower complexity can also be exploited to solve this problem. Both algorithms have O(N2) complexity to solve the Toeplitz-plus-Hankel system, as compared to the general algorithm which requires O(N3) complexity.


## 4. Computational Requirements Analysis

As observed from the above mentioned matrices, the required computations include addition, multiplication, and trigonometric operation. To obtain eachQ(i,j) from (12) requires 3L-1 additions, 3L multiplications, and 2L trigonometric operations. There are N2 distinct elements that need to be computed. Therefore, the total numbers of additions, multiplications, and trigonometric operations required to obtain Q are (3L-1)N2, 3LN2, and 2LN2, respectively. On the other hand, there are N distinct elements in d as shown in (20). Similarly, the total numbers of additions, multiplications, and trigonometric operations required to obtain d are (2L-1)N, 2LN, and 2LN, respectively.To obtain theN distinct elements of the Toeplitz matrix shown in (12), the total numbers of additions, multiplications, and trigonometric operations required to compute T are (2L-1)N, LN, and LN, respectively, for Kidambi’s method [13, 24]. In the proposed closed-form expressions of Toeplitz matrix shown in (17), the value of (i-j)Δ/2 needs to be computed only once and stored for further use. Evidently, each T(i,j) requires 3 additions, 6 multiplications, and 3 trigonometric operations. Therefore, the total numbers of additions, multiplications, and trigonometric operations required to obtain Q are 3N, 6N, and 3N, respectively.Similar computational analysis can be found on Hankel matrix. It is concluded that by using the closed-form expressions for the Toeplitz-plus-Hankel matrix, the computational cost can be reduced significantly. Table1 shows the comparison of the computations required for the matrices Q, H, T, and d. The total computation of the associated matrices is approximately 24LN-4N-7L and 49N operations for Kidambi approach [13] and the proposed technique, respectively. That is the closed-form expressions of the matrices involved here depend only on the length of the filter while those of the matrices involved in [13] depend on both the length of the filter and the number of sampling grid. In general, the number of the sampling grid is approximately L=kN with k≥4 for the compromise of integral squared error and computational requirements. It is notable that the required operations for the Toeplitz-plus-Hankel matrix of the proposed technique are O(N), whereas the Kidambi method [13] requires computational complexity of O(N2). On the other hand, the conventional least-squares method requires O(N3) complexity.Table 1
Comparison of computational requirements for associated matrices.
Matrix or vector
Operation
Algorithm
Conventional
Kidambi [13]
Proposed
Q
Addition
(
3
L
-
1
)
N
2
Multiplication
3
L
N
2
Sine or cosine
2
L
N
2
T
Addition
(
2
L
-
1
)
N
3
N
Multiplication
L
N
6
N
Sine or cosine
L
N
3
N
H
Addition
(
3
L
-
1
)
(
2
N
-
1
)
3
(
2
N
-
1
)
Multiplication
2
L
(
2
N
-
1
)
7
(
2
N
-
1
)
Sine or cosine
L
(
2
N
-
1
)
3
(
2
N
-
1
)
d
Addition
(
2
L
-
1
)
N
(
2
L
-
1
)
N
3
N
Multiplication
2
L
N
2
L
N
5
N
Sine or cosine
2
L
N
2
L
N
3
N
Total
Addition
3
L
N
2
-
N
2
+
2
L
N
-
N
10
L
N
-
3
L
-
4
N
+
1
12
N
-
3
Multiplication
3
L
N
2
+
2
L
N
7
L
N
-
3
L
25
N
-
7
Sine or cosine
2
L
N
2
+
2
L
N
5
L
N
-
L
12
N
-
3


## 5. Simulation and Design Examples

In this section, MATLAB programming language is used to design an IIR all-pass filter having the same specifications as that of [13] for evaluating the performance of the proposed technique. The simulations are evaluated on the IBM PC with an Intel Pentium IV-2.67 GHz CPU and 512 MB RAM.Example(Design of a Hilbert Transformer). For an IIR all-pass filter, the desired phase response is given by
(23)θd(ω)=-Nω-π2,0.08π≤ω≤0.92π,
where N=30 and the frequency sampling grid is set to L=10N. The design results including designed impulse response, magnitude response, phase response, phase error response, and group delay response are shown in Figure 1. As observed from the analysis of computations and accuracy for Toeplitz-plus-Hankel matrix, the proposed technique achieves the same performance as those of the Kidambi method [13] and the least-squares method, whereas the computational complexity can be reduced considerably. Comparison of the design performance with method of [13] and the least-squares method is indistinguishable; we just show the proposed design results. Figure 1(a) shows the design impulse response. The design magnitude response of the IIR all-pass filter is plotted in Figure 1(b). It is obvious that the proposed method achieves unit magnitude response at the bands of interest. The error response of magnitude response can be seen in Figure 1(c). The maximum magnitude peak error response is approximately 6.6613×10-15, which is very small. Figures 1(d) and 1(e) illustrate the phase response and error phase response, respectively. The constant desired group delay response and the designed one are compared at Figure 1(f). The maximum peak phase error response is 0.022 rad and the maximum peak group delay response is 0.5544. Moreover, the maximum pole radius is 0.9301, so the designed IIR all-pass filter not only is stable but also achieves excellent performance with effectiveness.A 30th Hilbert transformer. (a) Impulse response. (b) Designed magnitude response. (c) Designed magnitude error response. (d) Phase response. (e) Phase error response. (f) Group delay response.
(a)
(b)
(c)
(d)
(e)
(f)In order to compare the computational efficiency of the proposed technique, a Hilbert transformer withωo=0.04π, ωL=0.94π, and varying filter lengths from N=40 to N=240 are designed. Figure 2 shows the comparison ofCPU design time required for the computations of associated matrices with efficient method in [13]. It is evident that the proposed closed-form expressions for the computations of the Toeplitz-plus-Hankel matrix grow very slowly as the filter length increases. However, the design CPU time of the robust and efficient method [13] increases rapidly as compared to the proposed method when the filter is long. This is because the proposed closed-form expressions depend only on the filter length while the Kidambi [13] method depends both on filter length and on frequency sampling grid. Therefore, it is concluded that the proposed method has the superiority in efficiency as well as design accuracy.Figure 2
Comparison of CPU time for a Hilbert transformer with varying filter length.Although the Hilbert transformer is only considered for the derivation of the closed-form expression, the proposed method can also be applied to the quadrature mirror filter design that generally with a constant phase responseθd(ω)=-kω is required [26].


## 6. Conclusion

The least-squares design of all-pass IIR filter can be formulated by solving a system of linear equations and the associated matrix can be further represented as a Toeplitz-plus-Hankel form. This system of linear equations can be solved using efficient algorithms withO(N2) complexity. Based on the trigonometric identities and sampling the frequency domain uniformly, this paper presents a computationally more efficient and closed-form expressions for the simplification of Toeplitz-plus-Hankel matrix. As a result, not only the design accuracy can be achieved, but also the computational requirements can be reduced greatly.



---

*Source: 249021-2013-06-18.xml*