- Email: [email protected]
Fast design of 2D fully oversampled DFT modulated filter bank using Toeplitz-block Toeplitz matrix inversion
Signal Processing 111 (2015) 194–198
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Fast communication
Fast design of 2D fully oversampled DFT modulated filter bank using Toeplitz-block Toeplitz matrix inversion$ Fang Zhou a, Jun-Zheng Jiang b,n, Peng-Lang Shui a a b
National Lab. of Radar Signal Processing, Xidian University, Xi’an 710071, China School of Information and Communication, Guilin University of Electronic Technology, Guilin 541004, China
a r t i c l e in f o
abstract
Article history: Received 15 September 2014 Received in revised form 19 November 2014 Accepted 17 December 2014 Available online 24 December 2014
This letter presents a low-complexity algorithm to design two-dimensional (2D) DFT modulated filter bank (DMFB) with nearly perfect reconstruction (NPR). The design problem is formulated into an unconstrained optimization problem whose cost function consists of the overall distortion and stopband energy. By exploiting the gradient information, the prototype filter (PF) coefficients are iteratively optimized with closedform formula. At each iteration, the computational complexity is dramatically reduced by evoking Toeplitz-block Toeplitz matrix inversion and matrix inversion lemma, which makes the algorithm suitable for design 2D DMFB with a large number of channels and PF coefficients. Numerical example and comparison are included to show the effectiveness of the proposed design algorithm. & 2014 Elsevier B.V. All rights reserved.
Keywords: Fully oversampled Non-convex optimization Toeplitz-block Toeplitz matrix inversion Two-dimensional discrete Fourier transform modulated filter bank Two-dimensional filter
1. Introduction Design of two-dimensional (2D) filter banks presents several difficulties, which do not arise in one-dimensional (1D) ones [1]. Among 2D filter banks, 2D DFT modulated filter banks (DMFBs) have achieved more and more interests recently [2–7]. Compared to that of 1D counterpart [8,9], the design of 2D DMFB is much more challenging, particularly when requiring a large number of channels and PF coefficients (known as large-scale). In previous work, bi-iterative second-order cone program (BI-SOCP) algorithm [3] was proposed to design 2D DMFB with double-prototype, in which the design problem is converted into constrained problem. Restricted by the scale of the problems for the SOCP to be able to solve, the BI-SOCP algorithm cannot design 2D large-scale DMFBs. To overcome the deficiency,
☆ This work was supported by the National Natural Science Foundation of China under Grant no. 61261032 and Guangxi Natural Science Foundation under Grant no. 2013GXNSFBA019264. n Corresponding author. Tel./fax: þ 86 773 2290203. E-mail addresses: [email protected] (F. Zhou), [email protected] (J.-Z. Jiang), [email protected] (P.-L. Shui).
http://dx.doi.org/10.1016/j.sigpro.2014.12.021 0165-1684/& 2014 Elsevier B.V. All rights reserved.
the modified Newton’s method [5] and conjugate gradient method [7] were presented to design 2D DMFBs with singleprototype. These two methods can lead to filter banks with larger scale than the BI-SOCP algorithm. Nevertheless, both the methods still suffer from heavy computation cost in large-scale case. For the modified Newton’s method, a matrix with very large dimension was required to be inverted at each iteration. And the conjugate gradient method needed large number of iterations. Efficient design of 2D large-scale DMFB remains an intractable problem. Similar to that in [3,5,7], we focus on the design of 2D fully oversampled DMFB, owing to that only fully oversampled one can provide good overall performance [3]. A computationally efficient iterative algorithm based on unconstrained optimization is presented for the design. In the optimization problem, a quadratic objective function is formulated. With the gradient vector, the PF is optimized with closed-form formula. By jointly using fast inversion of the Toeplitz-block Toeplitz matrix [10] and matrix inverse lemma [11], the size of matrix to be inverted at each iteration is greatly reduced, thus resulting in considerable reduction in computational complexity. Furthermore, unlike the method [5], the proposed iterative algorithm utilizes unit step size, avoiding the
F. Zhou et al. / Signal Processing 111 (2015) 194–198
heavy complexity of calculating exact step size in [5]. Consequently, the proposed algorithm possesses much lower complexity than the method [5]. Numerical example and comparison with existing methods are included to demonstrate the efficiency and good performance of the proposed algorithm. 2. Prototype filter design Fig. 1 illustrates the structure of 2D DMFB with modulation matrix D1 and decimation matrix D2 . Assume hðnÞ; n A ½ L; L2 as the impulse response of the PF. Denote the 2D frequency variable as ω ¼ ½ωx ; ωy T , thenthe frequency response is HðωÞ ¼ cT ðω; LÞh with cðω; LÞ ¼ e j½ L; Lω ; …; e j½ L;Lω ; …; e j½L;Lω gT and h ¼ hð L; LÞ; …; hð L; LÞ; T …; hðL; LÞg . The analysis filters H i ðωÞ and synthesis ones Gi ðωÞ are generated by Gi ðωÞ ¼ H i ðωÞ ¼ Hðω 2πD1 T ui Þ;
ui A NðDT1 Þ;
i ¼ 0; 1; …; jD1 j 1:
ð1Þ where NðDT1 Þ ¼ fui DT1 xi A Z2 : xi A ½0; 1Þ2 ; i ¼ 0; 1; …; jD1 j 1g and jD1 j denotes the absolution of determinant of D1 . The subband and reconstructed signals are expressed as jD2 j 1 1 X Y i ðωÞ ¼ H ðω 2πD2 T vl ÞXðω 2πD2 T vl Þ jD2 j l ¼ 0 i
vl A NðDT2 Þ ¼ fvl DT2 xl A Z2 : xl A ½0; 1Þ2 ;
195
with fully oversampled [3], where the passband and stopband regions of the PF have no intersection. Here, we present several properties of the filter bank, which is very favor for design task. Forconvenience, denote ST ¼ njn mod Dn1 ¼ 0; n A ½ 2L; 2L2 and let set S T correspond o to the set lr A Z; 0 rlr rð4L þ 1Þ2 1; r ¼ 0; …; K 1 in row stack, where K is the size of ST . It can be found that K r ð16L2 =jD1 jÞ þ 8L þ 1 by invoking the configurations of the fundamental parallelepiped (FPD) of matrix D1 and region ½ 2L; 2L2 [12]. Property 2.1. [7]: The overall transfer function and aliasing transfer functions can be fast calculated with the following formula jD1 j T e ðωÞ AGðhÞΛk h ; k ¼ 0; …; jD2 j 1 T k ðωÞ ¼ ð5Þ jD2 j where GðhÞ is a block matrix defined by formula (11) in [5] and T eðωÞ ¼ expð jωT nl0 Þ; …; exp jωT nlK 1 ( A½r; m ¼
1;
m ¼ lr ; r ¼ 0; …; K 1; m ¼ 0; …; ð4L þ 1Þ2 1
0;
otherwise ð6bÞ
l ¼ 0; 1; …; jD2 j 1g
n o T T T T 1 T 1 T 1 Λk ¼ diag ej2πvk D2 ½ L; L ; …; ej2πvk D2 ½ L;L ; …; ej2πvk D2 ½L;L
ð2Þ jDX 2j 1
ð6aÞ
ð6cÞ T
ð3Þ
where v0 ¼ ½0; 0 . Based on (4a) and Property 2.1, we can derive the following property.
where overall transfer function T 0 ðωÞ and aliasing transfer functions T k ðωÞ are given by
Property 2.2. The 2D DMFB has no overall distortion if and only if
X~ ðωÞ ¼ T 0 ðωÞXðωÞ þ
T k ðωÞXðω 2πD2 T vk Þ
k¼1
T 0 ðωÞ ¼
jD1 j 1 1 X H ðωÞGi ðωÞ jD2 j i ¼ 0 i
T k ðωÞ ¼
jD1 j 1 1 X H ðω 2πD2 T vk ÞGi ðωÞ; jD2 j i ¼ 0 i
AGðhÞh ¼ b ð4aÞ
k ¼ 1; …; jD2 j 1 ð4bÞ
It is well-known perfect-reconstruction (PR) filter bank indicates that T 0 ðωÞ is identical to a pure delay and functions T k ðωÞ are equal to zero. And if the condition is satisfied with sufficient accuracy, the filter bank is nearly PR (NPR). In this letter, we consider to design 2D NPR DMFB
ð7Þ
where b is K 1 zero vector except its ðK 1Þ=2 element is jD2 j=jD1 j. The proof is similar to that in [5] and neglected. Similar to that in [5], overall performance can be controlled by minimizing the overall distortion and stopband energy. Consequently, the PF design problem can be formulated into the following unconstrained optimization problem with respect to the PF coefficients h. min ΦðhÞ ¼ ‖AGðhÞh b‖22 þ αEs ðhÞ ð8Þ h
where α is the weighted factor to get compromise between the terms in objective function and Es ðhÞ is the stopband energy of the PF, which is expressed as n o 2 T T Es ðhÞ ¼ ∬ΩS HðωÞ dω ¼ h ∬ΩS cðω; LÞc† ðω; LÞdω h ¼ h Rs h ð9Þ where stopband region is defined as Ωs ¼ fω A ½ π; πÞ2 and ω2 = SPDðπD2 T Þg with SPDðπD T Þ defined in [12] and superscript † denotes transpose-conjugate. Matrix Rs can be analytically calculated, shown in Appendix A. The gradient vector of the objective function ΦðhÞ is given by [5]
Fig. 1. Structure of 2D DFT modulated filter bank with the modulation matrix D1 and decimation matrix D2 .
∇ΦðhÞ ¼ 4ðAGðhÞÞT AGðhÞh b þ2αRs h
ð10Þ
196
F. Zhou et al. / Signal Processing 111 (2015) 194–198
By setting the gradient vector to zero, we obtain
1 ðAGðhÞÞT b h ¼ 4ðAGðhÞÞT AGðhÞ þ 2αRs
ð11Þ
Starting from an appropriate initial PF h0 , we can iteratively optimize the PF by using (11), where each iteration updates the PF as follows.
1 h ¼ 4ðAGðh0 ÞÞT AGðh0 Þ þ 2αRs ðAGðh0 ÞÞT b ð12Þ In terms of formula (150) in [11] (known as matrix inversion lemma), (12) can be rewritten as h ¼ QBðIK þ BT QBÞ 1 b=2; Q ¼ ð2αR s Þ 1 ; B ¼ 2ðAGðh0 ÞÞT ð13Þ The dimension of matrix inversion is changed from ð2L þ1Þ2 ð2L þ 1Þ2 to K K. In summary, the proposed algorithm is described as follows. Algorithm 2.3. Fast design of 2D fully oversampled DFT modulated filter bank Step 1. Specify the modulation matrix D1 and decimation matrix D2 (under fully oversampled condition), and PF support parameter L. Step 2. Design an initial PF h0 . Step 3. Given h0 , update h using formula (13). Step 4. Evaluate ‖h h0 ‖2 o η (η ¼ 1 10 3 in numerical example) or the iterative number exceeds a fixed bound C (C ¼ 20 in numerical example). If yes, terminate the iteration and output h as the final solution; Otherwise, set h0 ¼ 0:5ðh0 þ hÞ and go back to Step 3. In algorithm, the initial PF can be designed by the least squares method [5] or other conventional methods. For efficiency, the initial PF is fast designed by the MATLAB routine ‘fwind2’ in simulation. Some key features of Algorithm 2.3 are noteworthy. (1) Matrix Rs is a Toeplitz-block Toeplitz matrix, its inversion can therefore be efficiently computed by using the Wax–Kailath algorithm [10]. The complexity of inversion is Οðð2L þ 1Þ5 Þ, as compared to Οðð2L þ 1Þ6 Þ when conventional matrix inversion is used. Also, the inversion is independent on the PF and needs to be calculated only once. (2) In (11), 4ðAGðh0 ÞÞT AGðh0 Þ þ2αRs is a ð2L þ1Þ2 ð2L þ 1Þ2 symmetric positive matrix, whose order is equal to the number of PF coefficients, whereas IK þBT QB is a K K matrix whose order is the number of equations in (7). Recall that K r ð16L2 =jD1 jÞ þ 8L þ1, inversion is at least reduced the complexity of matrix 3 from Οðð2L þ 1Þ6 Þ to Ο ð16L2 =jD1 jÞ þ 8L þ 1 . Evidently, this computation saving is very substantial, particularly for large L and jD1 j. (3) According to the discussion in [6], Algorithm 2.3 is one type of modified Newton’s method and the iteration is at least convergent to a stationary point of the objective function ΦðhÞ, i.e., lim ∇Φðhk Þ ¼ 0
k-1
ð14Þ
(4) The proposed algorithm looks like the modified Newton’s method with exact line search [5], there are distinct differences, including the step size and the size of matrix to be inverted at each iteration. On one hand, at each step of the algorithm in [5], a ð2L2 þ2L þ 1Þ ð2L2 þ 2L þ1Þ is required to be inverted, whose complexity is Οðð2L2 þ 2L þ 1Þ3 Þ. However, only a K K be inverted whose complexity matrix is required to 3 is no more than Ο ð16L2 =jD1 jÞ þ8L þ 1 . On the other hand, for the method [5], the exact step is the roots of a cubic polynomial and the calculation of its coefficients involves 8Kð2L þ 1Þ2 þ ð10L2 þ 10L þ 8Þ ð2L2 þ 2L þ1Þ multiplications. Nevertheless, at each iteration of the proposed method, unit step is utilized to avoid high computational cost in calculating exact step. In summary, the complexity of matrix inversion at each iteration is considerably reduced.
3. Numerical examples In this section, Algorithm 2.3 is utilized to design 2D DMFB and comparison with other methods is also presented. Generally, the performance of the filter bank is measured by overall distortion (denoted by εt ), aliasing error (εa ), and stopband attenuation (SA). The reconstruction error is neglected since it is jointly determined by the overall distortion and aliasing error. All examples are running on a PC with Intel i7-core CPU and 4 GB memory. Example: In this example, a 2D large-scale DFT modulated filter bank is designed satisfying the specifications as follows:
D1 ¼
20
20
20
20
;
D2 ¼
10
10
10
10
;
L ¼ 50;
α ¼ 1 10 3
The filter bank has 800 subbands and the number of PF coefficients is up to 10,201. The modified Newton’s method, the conjugate gradient method, and the proposed algorithm are utilized to design the filter bank. For fair comparison, the three methods use the same initial PF and weighted factor in optimization problem. The magnitude responses of the resultant PFs are shown in Fig. 2. Table 1 presents the design performance, iterative times, and CPU time of the methods. It is observed that the three methods have comparable performance, while the proposed algorithm possesses much lower computational cost than the other ones. The low cost is attributed to the fast inversion of Toeplitz-block Toeplitz matrix and the considerable dimensional reduction of inverted matrices. On one hand, by using the Wax–Kailath algorithm, the inversion of Rs only costs 18 s, as compared to 188 s with traditional method. On the other hand, only a 61 61 matrix is inverted for the proposed algorithm, as compared to a 5101 5101 matrix inversion is required to be calculated at modified Newton’s iteration [5]. Also, the conjugate gradient method requires much more iterative times than the other two methods, although its single iteration cost is relatively low. In summary, the proposed algorithm is more suitable for design 2D large-scale DMFB than the other two methods. In addition, the BI-SOCP algorithm is
F. Zhou et al. / Signal Processing 111 (2015) 194–198
not suitable for designing the filter bank due to its huge computational cost. In optimization problem (8), α serves to achieve trade-off between the overall distortion and stopband attenuation. In order to verify this viewpoint, the filter bank is redesigned by the proposed algorithm with α ¼ 0:1. The iterative times is equivalent to that of α ¼ 1 10 3 , while the CPU time of the
197
two cases are slight different due to the negligible systematic residual. The performance measures of the generated filter bank are listed in Table 1. It is observed that the overall distortion gets worse with increased α, while stopband attenuation is improved. Generally, parameter α can be adjusted according to desired trade-off between the reconstruction property and stopband characteristic.
Fig. 2. Normalized magnitude responses of the PFs. (a) PF designed the modified Newton’s method. (b) PF designed by the conjugate gradient method. (c) PF designed by the proposed algorithm with α ¼ 1 10 3 . (d) PF designed by the proposed algorithm with α ¼ 0:1.
Table 1 Performance comparison in example. Performance measures
SA (dB)
εt (dB)
εa (dB)
CPU time per iteration (s)
Total time (s)
Iteration times
Modified Newton’s method Conjugate gradient method
46.63 45.56 48.43
54.98 57.43 51.63
71.13 66.94 67.06
65.03 2.23 4.02
2233 376 68.02
8 168 12
49.58
40.71
63.52
4.03
68.17
12
Proposed method (α ¼ 1 10 3 ) Proposed method (α ¼ 0:1)
Table 2 Complexity comparison between the method [5] and the proposed algorithm. Computational terms
Matrix inversion
Methods
Method [5] Proposed
Step size calculation
Method [5] Proposed
20 20
20 20
40 40
L ¼ 50
40 40 L ¼ 120
Ο 51053
Ο 613
Ο 1:35 108
Ο 290413
Ο 853
Ο 2:06 109
Ο 804013
Ο 1153
Ο 3:25 1010
0
0
0
D1 ¼
D1 ¼
40 80 L ¼ 200 D1 ¼
40 80
198
F. Zhou et al. / Signal Processing 111 (2015) 194–198
Finally, Table 2 lists the computational cost of the modified Newton’s method [5] and the proposed algorithm in large L cases, as well as the case in above example. It is observed that the proposed algorithm possesses much low computational cost than the method [5]. In addition, the complexity of the conjugate gradient method [7] is ΟðL4 Þ at each iteration. 4. Conclusion The letter proposes an iterative algorithm to design 2D DMFB with large-scale. By using Toeplitz-block Toeplitz matrix inversion and matrix inversion lemma, the proposed algorithm possesses much lower computational cost than the existing methods. Additionally, it should be noticed that the Wax-Kailath algorithm is only valid for Toeplitz-block Toeplitz matrix with good condition number, as claimed by the authors. Appendix A. Analytical calculation of matrix Rs The stopband energy matrix R s ¼ ∬ΩS cðω; LÞc† ðω; LÞdω ¼ ∬ð π;π Þ2 cðω; LÞc† ðω; LÞdω ∬SPDðπD T Þ cðω; LÞc† ðω; LÞdω 2
Recall cðω; LÞ ¼ e its kth element is
j½ L; Lω
; …; e
j½ L;Lω
; …; e
k L; ck ðω; LÞ ¼ expð jðkx ωx þ ky ωy ÞÞ; kx ¼ 2L þ 1
k L ky ¼ k ð2L þ 1Þ 2L þ1
¼ ð D2 ð0; 1Þu þ D2 ð0; 0ÞvÞ=jD2 j . Using Jacobi transformation, R2 is calculated as R2 ¼ ∬½ π;πÞ2 e jð½m=ð2L þ 1Þ ½½n=ð2L þ 1ÞÞ ðD2 ð1; 1Þu D2 ð1; 0ÞvÞ=jD2 j e jfðm nÞ ð2L þ1Þ m=ð2L þ 1Þ n=ð2L þ 1Þ 1 dudv ð D2 ð0; 1Þu þ D2 ð0; 0ÞvÞ=jD2 j jD2 j ¼
4π 2 sincðaÞ sincðbÞ: jD2 j
ðA5Þ
where m n D2 ð1; 1Þ a¼ jD2 j 2L þ1 2L þ 1 m n D2 ð0; 1Þ ðm nÞ ð2L þ 1Þ jD2 j 2L þ1 2L þ 1 m n D2 ð1; 0Þ b¼ jD2 j 2L þ 1 2L þ1 m n D2 ð0; 0Þ þ ðm nÞ ð2L þ 1Þ jD2 j 2L þ1 2L þ 1 ðA6Þ
,
With (A4) and (A5), the matrix Rs can be analytically constructed. Moreover, it can be found that R s is a Toeplitz-block Toeplitz matrix. References
ðA2Þ
Consequently, the ðm; nÞth element of matrix Rs is given R s ½m; n ¼ ∬ð π;πÞ2 cm ðω; LÞcnn ðω; LÞdω 2
ωx ¼ ðD2 ð1; 1Þu D2 ð1; 0ÞvÞ=jD2 j; ωy
ðA1Þ
j½L;Lω T
by ∬SPDðπD T Þ cm ðω; LÞcnn ðω; LÞdω ¼ R1 R2
Let ½u; vT ¼ DT2 ½ωx ; ωy T , i.e.,
ðA3Þ
The double integral R1 is separable with respect to the two frequency components ωx ; ωy and can be directly calculated as Z π Z π e jðmx nx Þωx dωx e jðmy ny Þωy dωy R1 ¼ π Zππ ¼ e jð½m=ð2L þ 1Þ ½n=ð2L þ 1ÞÞωx dωx Z ππ e j ðm nÞ ð2L þ 1Þð½m=ð2L þ 1Þ ½n=ð2L þ 1ÞÞ ωy dωy π m n ¼ 4π 2 sinc 2L þ 1 2L þ 1 m n ðA4Þ sinc ðm nÞ ð2L þ 1Þ 2L þ 1 2L þ1 The calculation of R2 is more difficult due to the integral domain is not separable. Here, we use Jacobi transformation to convert nonseparable domain into separable one.
[1] Y.P. Lin, P.P. Vaidyanathan, Theory and design of two-dimensional filter banks: a review, Multidimension. Syst. Signal Process. 7 (1996) 263–330. [2] P.L. Shui, Image denoising using 2-D separable oversampled DFT modulated filter banks, IET Image Process. 3 (3) (2009) 163–173. [3] J.Z. Jiang, P.L. Shui, Design of 2D linear phase DFT modulated filter banks using bi-iterative second-order cone program, Signal Process. 90 (2010) 3065–3077. [4] J.Z. Jiang, P.L. Shui, Two-dimensional 2 oversampled DFT modulated filter banks and critically sampled modified DFT modulated filter banks, IEEE Trans. Signal Process. 58 (11) (2010) 5597–5611. [5] J.Z. Jiang, P.L. Shui, Design of 2D oversampled linear phase DFT modulated filter banks via modified Newton’s method, Signal Process. 92 (2012) 1411–1421. [6] J.Z. Jiang, F. Zhou, Iterative design of two-dimensional critically sampled MDFT modulated filter banks, Signal Process. 93 (2013) 3124–3132. [7] J.Z. Jiang, F. Zhou, S. Ouyang, Design of two-dimensional large-scale DFT-modulated filter bank, IET Signal Process. 7 (9) (2013) 807–813. [8] H.H. Dam, S. Nordholm, Iterative method for the design of DFT filter bank, IEEE Trans. Circuits Syst. II: Express Briefs 51 (11) (2004) 581–586. [9] M.R. Wilbur, T.N. Davidson, J.P. Reilly, Efficient design of oversampled NPR GDFT filterbanks, IEEE Trans. Signal Process. 52 (7) (2004) 1947–1963. [10] M. Wax, T. Kailath, Efficient inversion of Toeplitz-block Toeplitz matrix, IEEE Trans. Acoust. Speech Signal Process. 31 (5) (1983) 1218–1221. [11] K.B. Petersen, M.S. Pedersen, The Matrix Cookbook (2008). [12] P.P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice-Hall, Englewood Cliffs, NJ, 1993.
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Fast communication
Fast design of 2D fully oversampled DFT modulated filter bank using Toeplitz-block Toeplitz matrix inversion$ Fang Zhou a, Jun-Zheng Jiang b,n, Peng-Lang Shui a a b
National Lab. of Radar Signal Processing, Xidian University, Xi’an 710071, China School of Information and Communication, Guilin University of Electronic Technology, Guilin 541004, China
a r t i c l e in f o
abstract
Article history: Received 15 September 2014 Received in revised form 19 November 2014 Accepted 17 December 2014 Available online 24 December 2014
This letter presents a low-complexity algorithm to design two-dimensional (2D) DFT modulated filter bank (DMFB) with nearly perfect reconstruction (NPR). The design problem is formulated into an unconstrained optimization problem whose cost function consists of the overall distortion and stopband energy. By exploiting the gradient information, the prototype filter (PF) coefficients are iteratively optimized with closedform formula. At each iteration, the computational complexity is dramatically reduced by evoking Toeplitz-block Toeplitz matrix inversion and matrix inversion lemma, which makes the algorithm suitable for design 2D DMFB with a large number of channels and PF coefficients. Numerical example and comparison are included to show the effectiveness of the proposed design algorithm. & 2014 Elsevier B.V. All rights reserved.
Keywords: Fully oversampled Non-convex optimization Toeplitz-block Toeplitz matrix inversion Two-dimensional discrete Fourier transform modulated filter bank Two-dimensional filter
1. Introduction Design of two-dimensional (2D) filter banks presents several difficulties, which do not arise in one-dimensional (1D) ones [1]. Among 2D filter banks, 2D DFT modulated filter banks (DMFBs) have achieved more and more interests recently [2–7]. Compared to that of 1D counterpart [8,9], the design of 2D DMFB is much more challenging, particularly when requiring a large number of channels and PF coefficients (known as large-scale). In previous work, bi-iterative second-order cone program (BI-SOCP) algorithm [3] was proposed to design 2D DMFB with double-prototype, in which the design problem is converted into constrained problem. Restricted by the scale of the problems for the SOCP to be able to solve, the BI-SOCP algorithm cannot design 2D large-scale DMFBs. To overcome the deficiency,
☆ This work was supported by the National Natural Science Foundation of China under Grant no. 61261032 and Guangxi Natural Science Foundation under Grant no. 2013GXNSFBA019264. n Corresponding author. Tel./fax: þ 86 773 2290203. E-mail addresses: [email protected] (F. Zhou), [email protected] (J.-Z. Jiang), [email protected] (P.-L. Shui).
http://dx.doi.org/10.1016/j.sigpro.2014.12.021 0165-1684/& 2014 Elsevier B.V. All rights reserved.
the modified Newton’s method [5] and conjugate gradient method [7] were presented to design 2D DMFBs with singleprototype. These two methods can lead to filter banks with larger scale than the BI-SOCP algorithm. Nevertheless, both the methods still suffer from heavy computation cost in large-scale case. For the modified Newton’s method, a matrix with very large dimension was required to be inverted at each iteration. And the conjugate gradient method needed large number of iterations. Efficient design of 2D large-scale DMFB remains an intractable problem. Similar to that in [3,5,7], we focus on the design of 2D fully oversampled DMFB, owing to that only fully oversampled one can provide good overall performance [3]. A computationally efficient iterative algorithm based on unconstrained optimization is presented for the design. In the optimization problem, a quadratic objective function is formulated. With the gradient vector, the PF is optimized with closed-form formula. By jointly using fast inversion of the Toeplitz-block Toeplitz matrix [10] and matrix inverse lemma [11], the size of matrix to be inverted at each iteration is greatly reduced, thus resulting in considerable reduction in computational complexity. Furthermore, unlike the method [5], the proposed iterative algorithm utilizes unit step size, avoiding the
F. Zhou et al. / Signal Processing 111 (2015) 194–198
heavy complexity of calculating exact step size in [5]. Consequently, the proposed algorithm possesses much lower complexity than the method [5]. Numerical example and comparison with existing methods are included to demonstrate the efficiency and good performance of the proposed algorithm. 2. Prototype filter design Fig. 1 illustrates the structure of 2D DMFB with modulation matrix D1 and decimation matrix D2 . Assume hðnÞ; n A ½ L; L2 as the impulse response of the PF. Denote the 2D frequency variable as ω ¼ ½ωx ; ωy T , thenthe frequency response is HðωÞ ¼ cT ðω; LÞh with cðω; LÞ ¼ e j½ L; Lω ; …; e j½ L;Lω ; …; e j½L;Lω gT and h ¼ hð L; LÞ; …; hð L; LÞ; T …; hðL; LÞg . The analysis filters H i ðωÞ and synthesis ones Gi ðωÞ are generated by Gi ðωÞ ¼ H i ðωÞ ¼ Hðω 2πD1 T ui Þ;
ui A NðDT1 Þ;
i ¼ 0; 1; …; jD1 j 1:
ð1Þ where NðDT1 Þ ¼ fui DT1 xi A Z2 : xi A ½0; 1Þ2 ; i ¼ 0; 1; …; jD1 j 1g and jD1 j denotes the absolution of determinant of D1 . The subband and reconstructed signals are expressed as jD2 j 1 1 X Y i ðωÞ ¼ H ðω 2πD2 T vl ÞXðω 2πD2 T vl Þ jD2 j l ¼ 0 i
vl A NðDT2 Þ ¼ fvl DT2 xl A Z2 : xl A ½0; 1Þ2 ;
195
with fully oversampled [3], where the passband and stopband regions of the PF have no intersection. Here, we present several properties of the filter bank, which is very favor for design task. Forconvenience, denote ST ¼ njn mod Dn1 ¼ 0; n A ½ 2L; 2L2 and let set S T correspond o to the set lr A Z; 0 rlr rð4L þ 1Þ2 1; r ¼ 0; …; K 1 in row stack, where K is the size of ST . It can be found that K r ð16L2 =jD1 jÞ þ 8L þ 1 by invoking the configurations of the fundamental parallelepiped (FPD) of matrix D1 and region ½ 2L; 2L2 [12]. Property 2.1. [7]: The overall transfer function and aliasing transfer functions can be fast calculated with the following formula jD1 j T e ðωÞ AGðhÞΛk h ; k ¼ 0; …; jD2 j 1 T k ðωÞ ¼ ð5Þ jD2 j where GðhÞ is a block matrix defined by formula (11) in [5] and T eðωÞ ¼ expð jωT nl0 Þ; …; exp jωT nlK 1 ( A½r; m ¼
1;
m ¼ lr ; r ¼ 0; …; K 1; m ¼ 0; …; ð4L þ 1Þ2 1
0;
otherwise ð6bÞ
l ¼ 0; 1; …; jD2 j 1g
n o T T T T 1 T 1 T 1 Λk ¼ diag ej2πvk D2 ½ L; L ; …; ej2πvk D2 ½ L;L ; …; ej2πvk D2 ½L;L
ð2Þ jDX 2j 1
ð6aÞ
ð6cÞ T
ð3Þ
where v0 ¼ ½0; 0 . Based on (4a) and Property 2.1, we can derive the following property.
where overall transfer function T 0 ðωÞ and aliasing transfer functions T k ðωÞ are given by
Property 2.2. The 2D DMFB has no overall distortion if and only if
X~ ðωÞ ¼ T 0 ðωÞXðωÞ þ
T k ðωÞXðω 2πD2 T vk Þ
k¼1
T 0 ðωÞ ¼
jD1 j 1 1 X H ðωÞGi ðωÞ jD2 j i ¼ 0 i
T k ðωÞ ¼
jD1 j 1 1 X H ðω 2πD2 T vk ÞGi ðωÞ; jD2 j i ¼ 0 i
AGðhÞh ¼ b ð4aÞ
k ¼ 1; …; jD2 j 1 ð4bÞ
It is well-known perfect-reconstruction (PR) filter bank indicates that T 0 ðωÞ is identical to a pure delay and functions T k ðωÞ are equal to zero. And if the condition is satisfied with sufficient accuracy, the filter bank is nearly PR (NPR). In this letter, we consider to design 2D NPR DMFB
ð7Þ
where b is K 1 zero vector except its ðK 1Þ=2 element is jD2 j=jD1 j. The proof is similar to that in [5] and neglected. Similar to that in [5], overall performance can be controlled by minimizing the overall distortion and stopband energy. Consequently, the PF design problem can be formulated into the following unconstrained optimization problem with respect to the PF coefficients h. min ΦðhÞ ¼ ‖AGðhÞh b‖22 þ αEs ðhÞ ð8Þ h
where α is the weighted factor to get compromise between the terms in objective function and Es ðhÞ is the stopband energy of the PF, which is expressed as n o 2 T T Es ðhÞ ¼ ∬ΩS HðωÞ dω ¼ h ∬ΩS cðω; LÞc† ðω; LÞdω h ¼ h Rs h ð9Þ where stopband region is defined as Ωs ¼ fω A ½ π; πÞ2 and ω2 = SPDðπD2 T Þg with SPDðπD T Þ defined in [12] and superscript † denotes transpose-conjugate. Matrix Rs can be analytically calculated, shown in Appendix A. The gradient vector of the objective function ΦðhÞ is given by [5]
Fig. 1. Structure of 2D DFT modulated filter bank with the modulation matrix D1 and decimation matrix D2 .
∇ΦðhÞ ¼ 4ðAGðhÞÞT AGðhÞh b þ2αRs h
ð10Þ
196
F. Zhou et al. / Signal Processing 111 (2015) 194–198
By setting the gradient vector to zero, we obtain
1 ðAGðhÞÞT b h ¼ 4ðAGðhÞÞT AGðhÞ þ 2αRs
ð11Þ
Starting from an appropriate initial PF h0 , we can iteratively optimize the PF by using (11), where each iteration updates the PF as follows.
1 h ¼ 4ðAGðh0 ÞÞT AGðh0 Þ þ 2αRs ðAGðh0 ÞÞT b ð12Þ In terms of formula (150) in [11] (known as matrix inversion lemma), (12) can be rewritten as h ¼ QBðIK þ BT QBÞ 1 b=2; Q ¼ ð2αR s Þ 1 ; B ¼ 2ðAGðh0 ÞÞT ð13Þ The dimension of matrix inversion is changed from ð2L þ1Þ2 ð2L þ 1Þ2 to K K. In summary, the proposed algorithm is described as follows. Algorithm 2.3. Fast design of 2D fully oversampled DFT modulated filter bank Step 1. Specify the modulation matrix D1 and decimation matrix D2 (under fully oversampled condition), and PF support parameter L. Step 2. Design an initial PF h0 . Step 3. Given h0 , update h using formula (13). Step 4. Evaluate ‖h h0 ‖2 o η (η ¼ 1 10 3 in numerical example) or the iterative number exceeds a fixed bound C (C ¼ 20 in numerical example). If yes, terminate the iteration and output h as the final solution; Otherwise, set h0 ¼ 0:5ðh0 þ hÞ and go back to Step 3. In algorithm, the initial PF can be designed by the least squares method [5] or other conventional methods. For efficiency, the initial PF is fast designed by the MATLAB routine ‘fwind2’ in simulation. Some key features of Algorithm 2.3 are noteworthy. (1) Matrix Rs is a Toeplitz-block Toeplitz matrix, its inversion can therefore be efficiently computed by using the Wax–Kailath algorithm [10]. The complexity of inversion is Οðð2L þ 1Þ5 Þ, as compared to Οðð2L þ 1Þ6 Þ when conventional matrix inversion is used. Also, the inversion is independent on the PF and needs to be calculated only once. (2) In (11), 4ðAGðh0 ÞÞT AGðh0 Þ þ2αRs is a ð2L þ1Þ2 ð2L þ 1Þ2 symmetric positive matrix, whose order is equal to the number of PF coefficients, whereas IK þBT QB is a K K matrix whose order is the number of equations in (7). Recall that K r ð16L2 =jD1 jÞ þ 8L þ1, inversion is at least reduced the complexity of matrix 3 from Οðð2L þ 1Þ6 Þ to Ο ð16L2 =jD1 jÞ þ 8L þ 1 . Evidently, this computation saving is very substantial, particularly for large L and jD1 j. (3) According to the discussion in [6], Algorithm 2.3 is one type of modified Newton’s method and the iteration is at least convergent to a stationary point of the objective function ΦðhÞ, i.e., lim ∇Φðhk Þ ¼ 0
k-1
ð14Þ
(4) The proposed algorithm looks like the modified Newton’s method with exact line search [5], there are distinct differences, including the step size and the size of matrix to be inverted at each iteration. On one hand, at each step of the algorithm in [5], a ð2L2 þ2L þ 1Þ ð2L2 þ 2L þ1Þ is required to be inverted, whose complexity is Οðð2L2 þ 2L þ 1Þ3 Þ. However, only a K K be inverted whose complexity matrix is required to 3 is no more than Ο ð16L2 =jD1 jÞ þ8L þ 1 . On the other hand, for the method [5], the exact step is the roots of a cubic polynomial and the calculation of its coefficients involves 8Kð2L þ 1Þ2 þ ð10L2 þ 10L þ 8Þ ð2L2 þ 2L þ1Þ multiplications. Nevertheless, at each iteration of the proposed method, unit step is utilized to avoid high computational cost in calculating exact step. In summary, the complexity of matrix inversion at each iteration is considerably reduced.
3. Numerical examples In this section, Algorithm 2.3 is utilized to design 2D DMFB and comparison with other methods is also presented. Generally, the performance of the filter bank is measured by overall distortion (denoted by εt ), aliasing error (εa ), and stopband attenuation (SA). The reconstruction error is neglected since it is jointly determined by the overall distortion and aliasing error. All examples are running on a PC with Intel i7-core CPU and 4 GB memory. Example: In this example, a 2D large-scale DFT modulated filter bank is designed satisfying the specifications as follows:
D1 ¼
20
20
20
20
;
D2 ¼
10
10
10
10
;
L ¼ 50;
α ¼ 1 10 3
The filter bank has 800 subbands and the number of PF coefficients is up to 10,201. The modified Newton’s method, the conjugate gradient method, and the proposed algorithm are utilized to design the filter bank. For fair comparison, the three methods use the same initial PF and weighted factor in optimization problem. The magnitude responses of the resultant PFs are shown in Fig. 2. Table 1 presents the design performance, iterative times, and CPU time of the methods. It is observed that the three methods have comparable performance, while the proposed algorithm possesses much lower computational cost than the other ones. The low cost is attributed to the fast inversion of Toeplitz-block Toeplitz matrix and the considerable dimensional reduction of inverted matrices. On one hand, by using the Wax–Kailath algorithm, the inversion of Rs only costs 18 s, as compared to 188 s with traditional method. On the other hand, only a 61 61 matrix is inverted for the proposed algorithm, as compared to a 5101 5101 matrix inversion is required to be calculated at modified Newton’s iteration [5]. Also, the conjugate gradient method requires much more iterative times than the other two methods, although its single iteration cost is relatively low. In summary, the proposed algorithm is more suitable for design 2D large-scale DMFB than the other two methods. In addition, the BI-SOCP algorithm is
F. Zhou et al. / Signal Processing 111 (2015) 194–198
not suitable for designing the filter bank due to its huge computational cost. In optimization problem (8), α serves to achieve trade-off between the overall distortion and stopband attenuation. In order to verify this viewpoint, the filter bank is redesigned by the proposed algorithm with α ¼ 0:1. The iterative times is equivalent to that of α ¼ 1 10 3 , while the CPU time of the
197
two cases are slight different due to the negligible systematic residual. The performance measures of the generated filter bank are listed in Table 1. It is observed that the overall distortion gets worse with increased α, while stopband attenuation is improved. Generally, parameter α can be adjusted according to desired trade-off between the reconstruction property and stopband characteristic.
Fig. 2. Normalized magnitude responses of the PFs. (a) PF designed the modified Newton’s method. (b) PF designed by the conjugate gradient method. (c) PF designed by the proposed algorithm with α ¼ 1 10 3 . (d) PF designed by the proposed algorithm with α ¼ 0:1.
Table 1 Performance comparison in example. Performance measures
SA (dB)
εt (dB)
εa (dB)
CPU time per iteration (s)
Total time (s)
Iteration times
Modified Newton’s method Conjugate gradient method
46.63 45.56 48.43
54.98 57.43 51.63
71.13 66.94 67.06
65.03 2.23 4.02
2233 376 68.02
8 168 12
49.58
40.71
63.52
4.03
68.17
12
Proposed method (α ¼ 1 10 3 ) Proposed method (α ¼ 0:1)
Table 2 Complexity comparison between the method [5] and the proposed algorithm. Computational terms
Matrix inversion
Methods
Method [5] Proposed
Step size calculation
Method [5] Proposed
20 20
20 20
40 40
L ¼ 50
40 40 L ¼ 120
Ο 51053
Ο 613
Ο 1:35 108
Ο 290413
Ο 853
Ο 2:06 109
Ο 804013
Ο 1153
Ο 3:25 1010
0
0
0
D1 ¼
D1 ¼
40 80 L ¼ 200 D1 ¼
40 80
198
F. Zhou et al. / Signal Processing 111 (2015) 194–198
Finally, Table 2 lists the computational cost of the modified Newton’s method [5] and the proposed algorithm in large L cases, as well as the case in above example. It is observed that the proposed algorithm possesses much low computational cost than the method [5]. In addition, the complexity of the conjugate gradient method [7] is ΟðL4 Þ at each iteration. 4. Conclusion The letter proposes an iterative algorithm to design 2D DMFB with large-scale. By using Toeplitz-block Toeplitz matrix inversion and matrix inversion lemma, the proposed algorithm possesses much lower computational cost than the existing methods. Additionally, it should be noticed that the Wax-Kailath algorithm is only valid for Toeplitz-block Toeplitz matrix with good condition number, as claimed by the authors. Appendix A. Analytical calculation of matrix Rs The stopband energy matrix R s ¼ ∬ΩS cðω; LÞc† ðω; LÞdω ¼ ∬ð π;π Þ2 cðω; LÞc† ðω; LÞdω ∬SPDðπD T Þ cðω; LÞc† ðω; LÞdω 2
Recall cðω; LÞ ¼ e its kth element is
j½ L; Lω
; …; e
j½ L;Lω
; …; e
k L; ck ðω; LÞ ¼ expð jðkx ωx þ ky ωy ÞÞ; kx ¼ 2L þ 1
k L ky ¼ k ð2L þ 1Þ 2L þ1
¼ ð D2 ð0; 1Þu þ D2 ð0; 0ÞvÞ=jD2 j . Using Jacobi transformation, R2 is calculated as R2 ¼ ∬½ π;πÞ2 e jð½m=ð2L þ 1Þ ½½n=ð2L þ 1ÞÞ ðD2 ð1; 1Þu D2 ð1; 0ÞvÞ=jD2 j e jfðm nÞ ð2L þ1Þ m=ð2L þ 1Þ n=ð2L þ 1Þ 1 dudv ð D2 ð0; 1Þu þ D2 ð0; 0ÞvÞ=jD2 j jD2 j ¼
4π 2 sincðaÞ sincðbÞ: jD2 j
ðA5Þ
where m n D2 ð1; 1Þ a¼ jD2 j 2L þ1 2L þ 1 m n D2 ð0; 1Þ ðm nÞ ð2L þ 1Þ jD2 j 2L þ1 2L þ 1 m n D2 ð1; 0Þ b¼ jD2 j 2L þ 1 2L þ1 m n D2 ð0; 0Þ þ ðm nÞ ð2L þ 1Þ jD2 j 2L þ1 2L þ 1 ðA6Þ
,
With (A4) and (A5), the matrix Rs can be analytically constructed. Moreover, it can be found that R s is a Toeplitz-block Toeplitz matrix. References
ðA2Þ
Consequently, the ðm; nÞth element of matrix Rs is given R s ½m; n ¼ ∬ð π;πÞ2 cm ðω; LÞcnn ðω; LÞdω 2
ωx ¼ ðD2 ð1; 1Þu D2 ð1; 0ÞvÞ=jD2 j; ωy
ðA1Þ
j½L;Lω T
by ∬SPDðπD T Þ cm ðω; LÞcnn ðω; LÞdω ¼ R1 R2
Let ½u; vT ¼ DT2 ½ωx ; ωy T , i.e.,
ðA3Þ
The double integral R1 is separable with respect to the two frequency components ωx ; ωy and can be directly calculated as Z π Z π e jðmx nx Þωx dωx e jðmy ny Þωy dωy R1 ¼ π Zππ ¼ e jð½m=ð2L þ 1Þ ½n=ð2L þ 1ÞÞωx dωx Z ππ e j ðm nÞ ð2L þ 1Þð½m=ð2L þ 1Þ ½n=ð2L þ 1ÞÞ ωy dωy π m n ¼ 4π 2 sinc 2L þ 1 2L þ 1 m n ðA4Þ sinc ðm nÞ ð2L þ 1Þ 2L þ 1 2L þ1 The calculation of R2 is more difficult due to the integral domain is not separable. Here, we use Jacobi transformation to convert nonseparable domain into separable one.
[1] Y.P. Lin, P.P. Vaidyanathan, Theory and design of two-dimensional filter banks: a review, Multidimension. Syst. Signal Process. 7 (1996) 263–330. [2] P.L. Shui, Image denoising using 2-D separable oversampled DFT modulated filter banks, IET Image Process. 3 (3) (2009) 163–173. [3] J.Z. Jiang, P.L. Shui, Design of 2D linear phase DFT modulated filter banks using bi-iterative second-order cone program, Signal Process. 90 (2010) 3065–3077. [4] J.Z. Jiang, P.L. Shui, Two-dimensional 2 oversampled DFT modulated filter banks and critically sampled modified DFT modulated filter banks, IEEE Trans. Signal Process. 58 (11) (2010) 5597–5611. [5] J.Z. Jiang, P.L. Shui, Design of 2D oversampled linear phase DFT modulated filter banks via modified Newton’s method, Signal Process. 92 (2012) 1411–1421. [6] J.Z. Jiang, F. Zhou, Iterative design of two-dimensional critically sampled MDFT modulated filter banks, Signal Process. 93 (2013) 3124–3132. [7] J.Z. Jiang, F. Zhou, S. Ouyang, Design of two-dimensional large-scale DFT-modulated filter bank, IET Signal Process. 7 (9) (2013) 807–813. [8] H.H. Dam, S. Nordholm, Iterative method for the design of DFT filter bank, IEEE Trans. Circuits Syst. II: Express Briefs 51 (11) (2004) 581–586. [9] M.R. Wilbur, T.N. Davidson, J.P. Reilly, Efficient design of oversampled NPR GDFT filterbanks, IEEE Trans. Signal Process. 52 (7) (2004) 1947–1963. [10] M. Wax, T. Kailath, Efficient inversion of Toeplitz-block Toeplitz matrix, IEEE Trans. Acoust. Speech Signal Process. 31 (5) (1983) 1218–1221. [11] K.B. Petersen, M.S. Pedersen, The Matrix Cookbook (2008). [12] P.P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice-Hall, Englewood Cliffs, NJ, 1993.