Construction of Bivariate Nonseparable Compactly Supported Orthogonal Wavelets

Abstract

A method for constructing bivariate nonseparable compactly supported orthogonal scaling functions, and the corresponding wavelets, using the dilation matrix 𝑀∶=2^𝑛𝕀=2^𝑛10₀₁, (𝑑=det𝑀=2^2𝑛≥4,𝑛∈ℕ) is given. The accuracy and smoothness of the scaling functions are studied, thus showing that they have the same accuracy order as the univariate Daubechies low-pass filter 𝑚₀(𝜔), to be used in our method. There follows that the wavelets can be made arbitrarily smooth by properly choosing the accuracy parameter 𝑟.

1. Introduction

Wavelet analysis (see [1, 2]) has become a powerful tool in neural networks, harmonic analysis, numerical analysis, and signal processing, especially in the area of image compression [3]. The wavelet transform is a simple and practical mathematical tool that cuts up data or functions into different frequency components and then studies each component with a resolution matched to its scale. The main feature of the wavelet transform is to hierarchically decompose general functions, as a signal or a process, into a set of approximation functions with different scales [1, 2, 4, 5].

Wavelet theory is closely related to subband coding, and it provides a functional space structure for subband coding, often leading to a better understanding of signals.

Univariate wavelets have found successful applications in signal processing [5–7]. However, in order to apply wavelet methods to data processing and to analyze multidimensional signals, such as image processing, we have to construct multivariate wavelets. The most commonly used method to define multivariate wavelets is the tensor product of univariate wavelets [4, 8–10]. This construction leads to a separable wavelet function having, as a major drawback, the horizontal and vertical (obliged) directions, while usually in imaging process the information contents spread along all directions, in a nonisotropic way. Much efforts have been spent for constructing multivariate nonseparable compactly supported orthogonal wavelets (see, e.g., [11–26]); however, the generalization to higher dimensions is not a simple task, so that a general simple methodology is still missing. However, it is possible to characterize the construction of bivariate scaling functions by the two scale difference equation, which in 1 dimension is 𝜑(𝑥)=_𝑘∈ℤ𝑐_𝑘𝜑(2𝑥−𝑘),(𝑥∈ℝ) (1) and in 2 dimensions 𝜑(𝐱)=_𝐤∈ℤ2𝜑(𝑀𝐱−𝐤),𝐱∈ℝ² (2) depends on the coefficients of a suitable second-order dilation matrix 𝑀. Usually [9, 11–13] the determinant of this matrix is 2. For example, Ayache [11, 12] proposed two methods for constructing bivariate nonseparable compactly supported orthonormal wavelets. The wavelets constructed using one of the methods are called semiseparable wavelets. Belogay and Wang [13] constructed bivariate nonseparable compactly supported orthonormal wavelets using the dilation matrix 𝑀=⁰²₁₀. Belogay and Wang gave the sufficient condition for the low-pass filters constructed by their method to generate nonseparable compactly supported orthogonal accurate scaling function with arbitrarily high smoothness. However, for proving the accuracy these authors used the matrix 𝑀²=±2𝕀. Moreover, their method has been modified by Lai [18] to construct nonseparable wavelets using the dilation matrix 2002 with determinant equal to 4. He and Lai [16] constructed many examples of nonseparable orthogonal wavelets. Lai and Roach [19] provided a method to construct bivariate nonseparable compactly supported wavelets by using the method of symmetry. Stanhill and Zeevi [24, 25] provided other methods for constructing bivariate nonseparable orthogonal wavelets.

By combining the Belogay and Wang method, the Lai method, and the Ayache method, we present another method for constructing bivariate nonseparable orthogonal wavelets using the dilation matrix 𝑀=2𝕀=²⁰₀₂. Our method holds true also for dilation matrix of the form 𝑀=2^𝑛𝕀=2^𝑛()1001𝑛∈ℕ (3) being 𝕀∶=¹⁰₀₁ the identity matrix, so that 𝑑=det𝑀=2^2𝑛(𝑛∈ℕ). (4) It follows that the two-scale difference equation simplifies into 𝜑(𝐱)=_𝐤∈ℤ2𝜑2^𝑛𝐱−𝐤𝐱∈ℝ²,𝑛∈ℕ (5) which is more alike the corresponding equation in 1 dimension. The minor drawback is that for the multivariate scaling function there exist 𝑑−1 corresponding wavelets (see, e.g., [9]), so that for a matrix with determinant 4 there exist, together with a scaling function, 3 wavelet functions (instead of 1 wavelet when 𝑑=2). However, due to the special form of the dilation matrix, we can take into account the results of [11–13] and show that the scaling function constructed by our method has the same accuracy as the univariate Daubechies low-pass filter to be used in our method. Following [13] we will also show that the smoothness of our bivariate scaling function can be manipulated by properly choosing the accuracy parameter 𝑟.

The paper is organized as follows. In Section 2, we introduce some notations and briefly recall some basic elements of the theory of bivariate multiresolution analysis and bivariate orthogonal wavelets. In Section 3, we provide a method for constructing a class of bivariate nonseparable compactly supported orthogonal filter banks. Then we study the accuracy and the smoothness of the scaling function constructed by our method. Finally, we give an example of the bivariate nonseparable compactly supported scaling function.

2. Preliminaries

Let ℤ and ℝbe the sets of all integers and real numbers, respectively. Denote 𝐿²ℝ²∶=𝑓(𝐱)∶_ℝ2𝑓(𝐱)²𝑑𝐱<∞,𝐱∈ℝ²,ℓ²ℤ²=𝐶_{𝐧𝐧∈ℤ2}∶_𝐧∈ℤ2𝐶_𝐧²<∞. (6) The Fourier transform of 𝑓(𝐱) is defined by 𝑓(𝝎)=_ℝ2𝑓(𝐱)𝑒^{−𝑖𝝎⋅𝐱}𝑑𝐱, (7) where 𝝎,𝐱∈ℝ²,𝑖=√−1, and 𝝎⋅𝐱∶=𝜔₁𝑥₁+𝜔₂𝑥₂ is the dot product.

For 𝑓(𝐱),𝑔(𝐱)∈𝐿²(ℝ²), their inner product is defined by 𝑓(𝐱),𝑔(𝐱)∶=_ℝ2𝑓(𝐱)𝑔(𝐱)𝑑𝐱, (8) which is explicitly 𝑓𝑥₁,𝑥₂,𝑔𝑥₁,𝑥₂∶=_ℝ2𝑓𝑥₁,𝑥₂𝑔𝑥₁,𝑥₂𝑑𝑥₁𝑑𝑥₂. (9)

A function 𝑓(𝐱)∈𝐿²(ℝ²) is called orthogonal if the set of its translate {𝑓(𝐱−𝐤)}_𝐤∈ℤ² is orthogonal; that is, 𝑓(𝐱−𝐤),𝑓(𝐱−𝐡)=𝛿_𝐤𝐡, (10) where 𝐤,𝐡∈ℤ² and 𝛿_𝐤𝐡 is the Kronecker notation defined by 𝛿_𝐤𝐡=1 if 𝐤=𝐡 and 𝛿_𝐤𝐡=0 if 𝐤≠𝐡, which is explicitly 𝑓𝑥₁−𝑘₁,𝑥₂−𝑘₂,𝑓𝑥₁−ℎ₁,𝑥₂−ℎ₂=𝛿_𝑘1ℎ1𝛿_𝑘2ℎ2. (11)

The whole construction of wavelets is based on the multiresolution analysis (MRA) and its corresponding axioms, which can be easily extended to the bivariate multiresolution analysis as follows. Let 𝜑(𝐱)∈𝐿²(ℝ²) satisfy the following refinement equation: 𝜑(𝐱)=_𝐤∈ℤ2𝑝_𝐤𝜑(2𝐱−𝐤), (12) where {𝑝_𝐤}_𝐤∈ℤ² is a real constant sequence.

We assume that 𝑝_𝐤≠0 for only finitely many 𝐤∈ℤ² (ensuring that 𝜑(𝐱) has compact support).

Define a closed subspace 𝑉_𝑗⊂𝐿²(ℝ²) by 𝑉_𝑗=clos_𝐿²_(ℝ²₎span2^𝑗𝜑2^𝑗𝐱−𝐤∶𝐤∈ℤ²,𝑗∈ℤ. (13)

Definition 1. One says that 𝜑(𝐱) in (12) generates a bivariate multiresolution analysis (MRA) {𝑉_𝑗}_𝑗∈ℤ of 𝐿²(ℝ²) if the sequence {𝑉_𝑗}_𝑗∈ℤ defined by (13) satisfies(i)⋯⊂𝑉₋₁⊂𝑉₀⊂𝑉₁⋯; (ii)⋃_𝑗∈ℤ𝑉_𝑗is dense in𝐿²(ℝ²) and ⋂_𝑗∈ℤ𝑉_𝑗={0}, where 0 is the zero function; (iii)𝑓(𝐱)∈𝑉_𝑗⇔𝑓(2𝐱)∈𝑉_𝑗+1,𝑗∈ℤ; (iv){𝜑(𝐱−𝐤)}_𝐤∈ℤ²is an orthogonal basis of 𝑉₀.The function 𝜑(𝐱) is called an orthogonal scaling function of the bivariate multiresolution analysis {𝑉_𝑗}_𝑗∈ℤ.

By taking the Fourier transform for both sides of (12), we have 𝜑𝜉,𝜂=𝑃₀𝑧₁,𝑧₂𝜑𝜉2,𝜂2,𝜉,𝜂∈ℝ, (14) where 𝑃₀𝑧₁,𝑧₂=_{()𝑘,𝑙∈ℤ2}𝑝_𝑘,𝑙𝑧^𝑘₁𝑧^𝑙₂,𝑧₁∶=𝑒^𝑖𝜉/2,𝑧₂∶=𝑒^𝑖𝜂/2, (15) is a Laurent polynomial satisfying 𝑃₀(1,1)=1.

The orthogonality condition (10) implies that _{𝑘1,𝑘2∈{}×{}0,10,1}𝑃₀(−1)^𝑘1𝑧₁,(−1)^𝑘2𝑧₂²=1. (16) The function 𝑃₀(𝑧₁,𝑧₂) is called a low-pass filter (refinement mask). By (14) and 𝑃₀(1,1)=1, we get 𝜑𝜉,𝜂=^∞_𝑘=1𝑃₀𝑒^{𝑖𝜉/2𝑘},𝑒^{𝑖𝜂/2𝑘}, (17) so that the scaling function can be defined by a suitable choice of the low-pass filter.

Let 𝐶(𝜉,𝜂)=𝑃₀(𝑧₁,𝑧₂); the filter 𝑃₀(𝑧₁,𝑧₂) satisfies Cohen’s criterion [26] if there exists a compact fundamental domain Ω of the lattice 2𝜋ℤ² with the property 𝐶𝜉2^𝑗,𝜂2^𝑗≠0,∀𝑗≥1,∀𝜉,𝜂∈Ω. (18)

Let us note that (16) is only a necessary condition for 𝜑(𝐱) to be orthogonal. The condition (16) becomes also a sufficient condition if the filter 𝑃₀(𝑧₁,𝑧₂) satisfies Cohen’s criterion (18).

Let us now assume that there exist three closed subspaces 𝑊_0,𝑖⊂𝐿²(ℝ²),𝑖=1,2,3, such that 𝑉₁=𝑉₀⊕𝑊_0,1⊕𝑊_0,2⊕𝑊_0,3, (19) where ⊕ denotes orthogonal direct sum of spaces. We say that the compactly supported functions 𝜓_𝑖(𝐱),𝑖=1,2,3 are orthogonal wavelets if {𝜓_𝑖(𝐱−𝐤)}_𝐤∈ℤ² form orthogonal bases of 𝑊_0,𝑖,𝑖=1,2,3; that is, 𝑊_0,𝑖=clos_𝐿²_(ℝ²₎span𝜓_𝑖(𝐱−𝐤)∶𝐤∈ℤ²𝜓,𝑖=1,2,3,_𝑖(𝐱−𝐤),𝜓_𝑗(𝐱−𝐥)=𝛿_𝑖,𝑗𝛿_𝐤,𝑙, (20) where 𝑖,𝑗=1,2,3 and 𝐤,𝑙∈ℤ².

If there exist three Laurent polynomials 𝑃_𝑖(𝑧₁,𝑧₂),𝑖=1,2,3, such that 𝑃𝑄=⁰𝑧₁,𝑧₂𝑃₀−𝑧₁,𝑧₂𝑃₀𝑧₁,−𝑧₂𝑃₀−𝑧₁,−𝑧₂𝑃₁𝑧₁,𝑧₂𝑃₁−𝑧₁,𝑧₂𝑃₁𝑧₁,−𝑧₂𝑃₁−𝑧₁,−𝑧₂𝑃₂𝑧₁,𝑧₂𝑃₂−𝑧₁,𝑧₂𝑃₂𝑧₁,−𝑧₂𝑃₂−𝑧₁,−𝑧₂𝑃₃𝑧₁,𝑧₂𝑃₃−𝑧₁,𝑧₂𝑃₃𝑧₁,−𝑧₂𝑃₃−𝑧₁,−𝑧₂ (21) is unitary, that is, det𝑄=1, then the Fourier transforms of the functions 𝜓_𝑖(𝐱),𝑖=1,2,3, can be given by 𝜓_𝑖𝜉,𝜂=𝑃_𝑖𝑒^𝑖𝜉/2,𝑒^𝑖𝜂/2𝜑𝜉2,𝜂2, (22) and the three Laurent polynomials 𝑃_𝑖(𝑧₁,𝑧₂),𝑖=1,2,3, are called the high-pass filters (wavelet masks). The whole set of functions {𝑃_𝑖(𝑧₁,𝑧₂)}_0≤𝑖≤3 is called a filter bank.

By (17) and (22), the construction of bivariate compactly supported orthogonal scaling function and the corresponding wavelets is reduced to the construction of filter bank {𝑃_𝑖(𝑧₁,𝑧₂)}_0≤𝑖≤3 which satisfies (16), (18), and (21).

3. Construction of Bivariate Nonseparable Orthogonal Wavelets

3.1. Construction of the Filter Banks

According to Section 2, for constructing bivariate orthogonal scaling function and the corresponding wavelets, we mainly need to define the bivariate orthogonal filter bank {𝑃_𝑖(𝑧₁,𝑧₂)}_0≤𝑖≤3 satisfying (16), (18), and (21).

Let us shortly recall the Belogay and Wang method, the Lai method, and the Ayache method for constructing bivariate nonseparable compactly supported orthogonal wavelets and their corresponding filter banks.

Belogay and Wang [13] constructed the bivariate nonseparable orthogonal scaling functions using the dilation matrix 0210 with special coefficients: the coefficients are aligned along two adjacent rows. The corresponding low-pass filters have the form 𝑃₀𝑧₁,𝑧₂=𝐴𝑧₁+𝑧₂𝐵𝑧₁, (23) where 𝐴(𝑧),𝐵(𝑧) are two univariate polynomials of 𝑧=𝑒^𝑖𝜉 satisfying some conditions (see [13]). The high-pass filter coefficients 𝑑_𝐧 are given by 𝑑_𝐧=(−1)^𝑛1𝑐_𝐞−𝐧, (24) where 𝐧=(𝑛₁,𝑛₂)∈ℤ²,𝐞=(1,0), and 𝑐_𝐧 are low-pass filter coefficients.

Lai [18] modified the Belogay and Wang method to construct nonseparable wavelets using the dilation matrix 2002. According to the Lai method, their filter bank has the form 𝑃₀𝑧₁,𝑧₂=𝐴𝑧₁+𝑧₂𝐵𝑧₁𝐴𝑧₂+𝑧²₁𝐵𝑧₂,𝑃₁𝑧₁,𝑧₂=𝑧₂𝐴𝑧₁+𝑧₂𝐵𝑧₁𝐴−𝑧₂⁻¹+𝑧²₁𝐵−𝑧₂⁻¹,𝑃₂𝑧₁,𝑧₂=𝑧₁𝐴−𝑧₁⁻¹+𝑧₂⁻¹𝐵−𝑧₁⁻¹×𝐴𝑧₂+𝑧²₁𝐵𝑧₂,𝑃₃𝑧₁,𝑧₂=𝑧₁𝑧₂𝐴−𝑧₁⁻¹+𝑧₂⁻¹𝐵−𝑧₁⁻¹×𝐴−𝑧₂⁻¹+𝑧²₁𝐵−𝑧₂⁻¹, (25) where 𝐴(𝑧),𝐵(𝑧) are two univariate polynomials of 𝑧=𝑒^𝑖𝜉 satisfying some conditions, 𝑧₁=𝑒^𝑖𝜉, and 𝑧₂=𝑒^𝑖𝜂 (see [18]).

Ayache [11, 12] proposed two methods for constructing bivariate compactly supported orthogonal wavelets using the dilation matrix 2002. Let 𝛼(𝑥),𝛽(𝑥) be two univariate trigonometric polynomials, such that 𝛼(0)=1,|𝛼(𝑥)|²+𝛽(𝑥)²=1. (26) Then according to one of the two methods, the nonseparable filter banks have the form 𝑃₀𝑧₁,𝑧₂=𝑢₀𝑧₁𝛼2𝜂𝑣₀𝑧₂+𝛽2𝜂𝑣₁𝑧₂𝑧₁⁻²,𝑃₁𝑧₁,𝑧₂=𝑢₀𝑧₁𝛽2𝜂𝑣₀𝑧₂−𝛼2𝜂𝑣₁𝑧₂𝑧₁⁻²,𝑃₂𝑧₁,𝑧₂=𝑢₁𝑧₁𝛼2𝜂𝑣₀𝑧₂+𝛽2𝜂𝑣₁𝑧₂𝑧₁⁻²,𝑃₃𝑧₁,𝑧₂=𝑢₁𝑧₁𝛽2𝜂𝑣₀𝑧₂−𝛼2𝜂𝑣₁𝑧₂𝑧₁⁻², (27) where {𝑢_𝑘(𝑧)}_0≤𝑘≤1,{𝑣_𝑘(𝑧)}_0≤𝑘≤1 are two given univariate filter banks, 𝑧₁=𝑒^𝑖𝜉, and 𝑧₂=𝑒^𝑖𝜂.

By combining the previous three methods, we propose, in the following, a method based on a suitable choice of the univariate filter bank, together with the basic ideas of Belogay-Wang approach. In particular, we use the univariate Daubechies filter bank to define the bivariate nonseparable filter bank as follows.

Let 𝑟 be a nonnegative integer, and let {𝑚_𝑖(𝑧)}_0≤𝑖≤1 be the univariate Daubechies filter bank [2, 27] with accuracy 𝑟+1; that is, 𝑚₀(𝑧)=1+𝑧2^𝑟+1𝐻_𝑟+1(𝑧),𝑚₁(𝑧)=𝑒^−𝑖𝜉𝑚₀−𝑧⁻¹, (28) where 𝐻_𝑟+1(𝑧)²∶=^𝑟_𝑘=0𝑘𝑟+𝑘sin²𝜉2^𝑘 (29) and 𝑧=𝑒^𝑖𝜉. Let 𝐴(𝑧),𝐵(𝑧) be two univariate polynomials of 𝑧=𝑒^𝑖𝜉 with real coefficients, and let 𝐴(1)≠0,𝐵(0)≠0. We define four bivariate polynomials of 𝑧₁=𝑒^𝑖𝜉, 𝑧₂=𝑒^𝑖𝜂 which are the forms of the filter banks of this paper, as follows: 𝑃₀𝑧₁,𝑧₂∶=𝑚₀𝑧₂𝐴𝑧₁+𝑧²₂𝐵𝑧₁𝑃₁,𝑧₁,𝑧₂∶=𝑧₁𝑚₀𝑧₂𝐵−𝑧₁⁻¹+𝑧²₂𝐴−𝑧₁⁻¹𝑃₂,𝑧₁,𝑧₂∶=𝑚₁𝑧₂𝐴𝑧₁+𝑧²₂𝐵𝑧₁𝑃₃,𝑧₁,𝑧₂∶=𝑧₁𝑚₁𝑧₂𝐵−𝑧₁⁻¹+𝑧²₂𝐴−𝑧₁⁻¹. (30)(31)(32)(33)

Let 𝑟 be a nonnegative integer, and let 𝜈 be an odd integer with 𝜈≥deg𝐴(𝑧). Now, we define the polynomials 𝐴(𝑧) and 𝐵(𝑧) of 𝑧=𝑒^𝑖𝜉 as follows: 𝐴(𝑧)∶=𝑧^𝜈1+𝑧⁻¹2^𝑟𝐿𝑧⁻¹𝑆𝑧⁻²,𝐵(𝑧)∶=1+𝑧2^𝑟1−𝑧2^2𝑟𝐿(−𝑧)𝑄𝑧², (34) where 𝐿(𝑧), 𝑆(𝑧), and 𝑄(𝑧) are given by 𝑆𝑧²²=1−1−𝑧²4^𝑟1−𝑧⁻²4^𝑟𝑄𝑧²²,|𝐿(𝑧)|²=^𝑟−1_𝑘=0𝑘𝑟+𝑘−11−𝑢2^𝑘+(1−𝑢)^𝑟𝑢𝑅𝑢²,1𝑢∶=2𝑧+𝑧⁻¹,𝑆(1)=𝐿(1)=1,𝑄(1)=(−1)^𝑟𝐿(),𝐿(0)𝑄(0)−1≠0, (35)(36)(37) where 𝑅(𝑧), 𝑄(𝑧) are two arbitrary polynomials and 𝑄(𝑧) satisfies (37).

The following condition (to be used in the following) is given for the minimal degree of 𝐴(𝑧) and 𝐵(𝑧): 𝑄(𝑧)=constant=(−1)^𝑟𝐿(),𝑅(𝑧)−1=constant=0,𝜈=4𝑟−1. (38)

Remark 2. Obviously, the polynomial 𝐿(𝑧) defined by (36) satisfies 1+𝑧2^2𝑟|𝐿(𝑧)|²+1−𝑧2^2𝑟|𝐿(−𝑧)|²=1. (39) This equation was solved by Daubechies in [27]. If (38) holds, then 𝐿(𝑧)=𝐻_𝑟(𝑧).

Lemma 3. The polynomial 𝑃₀(𝑧₁,𝑧₂) given by (30), where 𝑚₀(𝑧) is the univariate Daubechies filer defined by (28), satisfies (16) if and only if the polynomials 𝐴(𝑧), 𝐵(𝑧) in (30) satisfy |𝐴(𝑧)|²+|𝐴(−𝑧)|²+|𝐵(𝑧)|²+|𝐵(−𝑧)|²𝐴𝑧=1,₋₁𝐵(𝑧)+𝐴−𝑧⁻¹𝐵(−𝑧)=0. (40)(41)

Proof. Note that 𝑃₀𝑧₁,𝑧₂²=𝑚₀𝑧₂²×𝐴𝑧₁²+𝐵𝑧₁²+𝑧₂⁻²𝐴𝑧₁𝐵𝑧₁⁻¹+𝑧²₂𝐴𝑧₁⁻¹𝐵𝑧₁,𝑚₀(𝑧)²+𝑚₀(−𝑧)²=1. (42) If (40) and (41) hold, we have _{𝑘1,𝑘2∈{}×{}0,10,1}𝑃₀(−1)^𝑘1𝑧₁,(−1)^𝑘2𝑧₂²=𝑚₀𝑧₂²𝐴𝑧₁²+𝐵𝑧₁²+𝐴−𝑧₁²+𝐵−𝑧₁²+𝑧₂⁻²×𝐴𝑧₁𝐵𝑧₁⁻¹+𝐴−𝑧₁𝐵−𝑧₁⁻¹+𝑧²₂𝐴𝑧₁⁻¹𝐵𝑧₁+𝐴−𝑧₁⁻¹𝐵−𝑧₁+𝑚₀−𝑧₂²𝐴𝑧₁²+𝐵𝑧₁²𝐴−𝑧₁²+𝐵−𝑧₁²+𝑧₂⁻²×𝐴𝑧₁𝐵𝑧₁⁻¹+𝐴−𝑧₁𝐵−𝑧₁⁻¹+𝑧²₂𝐴𝑧₁⁻¹𝐵𝑧₁+𝐴−𝑧₁⁻¹𝐵−𝑧₁=𝑚₀𝑧₂²+𝑚₀−𝑧₂²=1. (43) Hence (16) holds.
Conversely, if (16) holds, we have _{𝑘1,𝑘2∈{}×{}0,10,1}𝑃₀(−1)^𝑘1𝑧₁,(−1)^𝑘2𝑧₂²=𝑚₀𝑧₂²+𝑚₀−𝑧₂²×𝐴𝑧₁²+𝐵𝑧₁²+𝑧₂⁻²𝐴𝑧₁𝐵𝑧₁⁻¹+𝑧²₂𝐴𝑧₁⁻¹𝐵𝑧₁⁻¹+𝑚₀𝑧₂²+𝑚₀−𝑧₂²×𝐴−𝑧₁²+𝐵−𝑧₁²+𝑧₂⁻²𝐴−𝑧₁𝐵−𝑧₁⁻¹+𝑧²₂𝐴−𝑧₁⁻¹𝐵−𝑧₁⁻¹=𝐴𝑧₁²+𝐵𝑧₁²+𝐴−𝑧₁²+𝐵−𝑧₁²+𝑧₂⁻²𝐴𝑧₁𝐵𝑧₁⁻¹+𝐴−𝑧₁𝐵−𝑧₁⁻¹×𝑧²₂𝐴𝑧₁⁻¹𝐵𝑧₁+𝐴−𝑧₁⁻¹𝐵−𝑧₁=1, (44) for all |𝑧₁|=1 and |𝑧₂|=1. Hence (40) and (41) hold.

Remark 4. Condition (41) means that 𝐴(𝑧⁻¹)𝐵(𝑧) contains no even powers of 𝑧. Therefore, deg𝐴(𝑧) must be odd if 𝐵(0)≠0.

Lemma 5. The following conditions are equivalent: (i)the univariate real polynomials 𝐴(𝑧) and 𝐵(𝑧) satisfy condition (41) and 𝐵(0)≠0;(ii)there exist an odd integer number 𝜈≥deg𝐴(𝑧) and real polynomials 𝑠(𝑧), 𝑞(𝑧), and 𝑙(𝑧), such that 𝑧^𝜈𝐴(𝑧⁻¹)=𝑠(𝑧²)𝑙(𝑧), 𝐵(𝑧)=𝑞(𝑧²)𝑙(−𝑧), and 𝑔𝑐𝑑(𝑙(𝑧),𝑙(−𝑧))=1.

Proof. See the proof of Lemmas 3.2 and 3.3 of [13].

Lemma 6. Let the polynomials 𝐴(𝑧) and 𝐵(𝑧) be given by (34)–(38) for 𝑧=𝑒^𝑖𝜉; then 𝐴𝑒^𝑖𝜉+𝑒^2𝑖𝜂𝐵𝑒^𝑖𝜉=0 (45) if and only if 𝜉≡𝜋(𝑚𝑜𝑑2𝜋).

Proof. See the proof of Lemma 4.1 and Lemma 4.7 of [13].

By collecting the previous results and definition we can easily show the main result of this paper that is as follows.

Theorem 7 (existence). Let 𝑟 and 𝜈 be, respectively, a nonnegative integer and an odd integer number with 𝜈≥deg𝐴(𝑧); let the polynomials 𝐴(𝑧) and 𝐵(𝑧) be given by (34)–(37); let the filter 𝑚₀(𝑧) be given by (28); let the polynomials 𝑃₀(𝑧₁,𝑧₂), 𝑃₁(𝑧₁,𝑧₂), 𝑃₂(𝑧₁,𝑧₂), and 𝑃₃(𝑧₁,𝑧₂) be defined by (30)–(33). Then {𝑃_𝑘(𝑧₁,𝑧₂)}_0≤𝑘≤3 is a bivariate nonseparable filter bank satisfying (16) and (21). If, in addition, the polynomials 𝐴(𝑧) and 𝐵(𝑧) satisfy condition (38), then 𝑃₀(𝑧₁,𝑧₂) satisfies Cohen’s criterion (18), and the filter bank {𝑃_𝑘(𝑧₁,𝑧₂)}_0≤𝑘≤3 can be used, according to (17), (22), to define bivariate nonseparable compactly supported orthogonal wavelets.

Proof. According to Lemmas 3 and 5 by letting 𝑠𝑧²=𝑆𝑧²,𝑞𝑧²=1−𝑧²2^𝑟𝑄𝑧²,𝑙(𝑧)=1+𝑧2^𝑟𝐿(𝑧), (46) in order to prove that 𝑃₀(𝑧₁,𝑧₂) satisfies (16), we need only to show that 𝐴(𝑧), 𝐵(𝑧) satisfy (40). Taking into account (34) and (35), we have |𝐴(𝑧)|²+|𝐵(−𝑧)|²=1+𝑧2^2𝑟|𝐿(𝑧)|²𝑆𝑧²²+1−𝑧2^2𝑟1+𝑧2^4𝑟|𝐿(𝑧)|²𝑄𝑧²²=1+𝑧2^2𝑟|𝐿(𝑧)|²𝑆𝑧²²+1−𝑧2^2𝑟1+𝑧2^2𝑟𝑄𝑧²²=1+𝑧2^2𝑟|𝐿(𝑧)|². (47) Hence, (40) follows from (39).
By a simple calculation, we can also get that {𝑃_𝑘(𝑧₁,𝑧₂)}_0≤𝑘≤3 satisfies (21). In fact, let us assume that 𝑃₀𝑧₁,𝑧₂=𝑢𝑧₁𝑣𝑧₂, (48) where 𝑢(𝑧) and 𝑣(𝑧) are univariate low-pass filters. Then we have 𝑃₀1,𝑧₂=𝑣𝑧₂=𝑚₀𝑧₂. (49) This leads to the contradiction that 𝑢𝑧₁=𝐴𝑧₁+𝑧²₂𝐵𝑧₁. (50) Hence 𝑃₀(𝑧₁,𝑧₂) is nonseparable.
Finally, suppose that 𝐴(𝑧) and 𝐵(𝑧) satisfy (34)–(38). Let Ω=[−𝜋,𝜋]×[−𝜋,𝜋], then Ω is a compact fundamental domain of the lattice 2𝜋ℤ². When 𝑗≥1 and (𝜉,𝜂)∈Ω, it is (𝜉2^𝑗,𝜂2^𝑗𝜋)∈[−2,𝜋2]², so that we have 𝐴𝑒^{𝑖𝜉/2𝑗}+𝑒^{2𝑖𝜂/2𝑗}𝐵𝑒^{𝑖𝜉/2𝑗}≠0,𝑚₀𝑒^{𝑖𝜂/2𝑗}≠0, (51) for all 𝑗≥1 and (𝜉,𝜂)∈Ω by using Lemma 6 and the result of [27]. Hence 𝑃₀(𝑧₁,𝑧₂) satisfies Cohen’s criterion (18).

3.2. Accuracy

We will use the same method as given in [13] to study the accuracy of our scaling function. Let recall its definition and how the accuracy of the scaling function implies some restrictions on the low-pass filter.

According to the result of [14], an orthogonal scaling function 𝜑(𝐱) with the dilation matrix 2002 has accuracy 𝑟+1 if and only if its low-pass filter 𝑃₀(𝑧₁,𝑧₂) satisfies the following accuracy condition: 𝜕^𝑖+𝑗𝜕𝑧^𝑖₁𝜕𝑧^𝑗₂𝑃₀𝑧₁,𝑧₂∣_{𝑧1=−1,𝑧2=1}𝜕=0,_𝑖+𝑗𝜕𝑧^𝑖₁𝜕𝑧^𝑗₂𝑃₀𝑧₁,𝑧₂∣_{𝑧1=1,𝑧2=−1}𝜕=0,_𝑖+𝑗𝜕𝑧^𝑖₁𝜕𝑧^𝑗₂𝑃₀𝑧₁,𝑧₂∣_{𝑧1=−1,𝑧2=−1}=0, (52) for all 𝑖, 𝑗≥0 with 𝑖+𝑗≤𝑟. Therefore, when 𝑃₀(𝑧₁,𝑧₂) satisfies (52), we say that it has accuracy 𝑟+1.

Lemma 8. Let the polynomials 𝑚₀(𝑧) and 𝑃₀(𝑧₁,𝑧₂) be given by (28) and (30), respectively; then 𝑃₀(𝑧₁,𝑧₂) satisfies (52) if and only if the polynomials 𝐴(𝑧),𝐵(𝑧) in (30) defined by (34) have the form 𝐴(𝑧)=1+𝑧2^𝑟𝐴₀(𝑧),𝐵(𝑧)=1+𝑧2^𝑟𝐵₀(𝑧), (53) where 𝐴₀(𝑧), 𝐵₀(𝑧) are polynomials with real coefficients satisfying 𝐴₀(−1)+𝐵₀(−1)=0. (54)

Proof. It is easy to check that 𝜕^𝑖𝑃₀𝑧₁,𝑧₂𝜕𝑧^𝑖₁=𝑚₀𝑧₂𝐴^(𝑖)𝑧₁+𝑧²₂𝐵^(𝑖)𝑧₁,𝜕^𝑖+1𝑃₀𝑧₁,𝑧₂𝜕𝑧^𝑖₁𝜕𝑧₂=𝑚^′₀𝑧₂𝐴^(𝑖)𝑧₁+𝑧²₂𝐵^(𝑖)𝑧₁+2𝑧₂𝑚₀𝑧₂𝐵^(𝑖)𝑧₁,⋮𝜕^𝑖+𝑗𝑃₀𝑧₁,𝑧₂𝜕𝑧^𝑖₁𝜕𝑧^𝑗₂=𝑚₀^(𝑗)𝑧₂𝐴^(𝑖)𝑧₁+𝑧²₂𝐵^(𝑖)𝑧₁+2𝑗𝑧₂𝑚₀^(𝑗−1)𝑧₂𝐵^(𝑖)𝑧₁+𝑗𝑗−1𝑚₀^(𝑗−2)𝑧₂𝐵^(𝑖)𝑧₁, (55) for all 𝑖≥0 and 𝑗≥2.
If 𝐴(𝑧), 𝐵(𝑧) satisfy (53) and (54), it is not difficult to show that 𝑃₀(𝑧₁,𝑧₂) satisfies also (52). On the other hand, if 𝑃₀(𝑧₁,𝑧₂) satisfies (52) and noting that 𝑚₀(1)=1, we have 𝐴^(𝑖)(−1)+𝐵^(𝑖)()𝐴−1=0,for0≤𝑖≤𝑟,_(𝑖)(−1)=0,𝐵^(𝑖)(−1)=0,for0≤𝑖<𝑟. (56) So that 𝑧=−1 is a root of multiplicity 𝑟 of 𝐴(𝑧),𝐵(𝑧). Hence 𝐴(𝑧),𝐵(𝑧) have the form of (53), and 𝐴^(𝑟)(−1)+𝐵^(𝑟)(−1)=0 yields (54).

Lemma 9. Let the polynomials 𝑚₀(𝑧) and 𝑃₀(𝑧₁,𝑧₂) be given by (28) and (30), respectively, and assume also that 𝐴(𝑧) and 𝐵(𝑧) in (30) are polynomials with 𝐴(1)≠0 and 𝐵(0)≠0. Then 𝑃₀(𝑧₁,𝑧₂) satisfies (41) with accuracy 𝑟+1 if and only if there exist an odd integer number 𝜈≥deg𝐴(𝑧) and real polynomials 𝑆(𝑧), 𝑄(𝑧), and 𝐿(𝑧), such that 𝐴(𝑧) and 𝐵(𝑧) satisfy (34) and (37).

Proof. If there exist an odd integer 𝜈≥deg𝐴(𝑧) and real polynomials 𝑆(𝑧),𝑄(𝑧), and 𝐿(𝑧) such that 𝐴(𝑧) and 𝐵(𝑧) satisfy (34), and (37), then 𝑃₀(𝑧₁,𝑧₂) satisfies (41) and has accuracy 𝑟+1 by Lemmas 5 and 8 by letting 𝑠𝑧²=𝑆𝑧²,𝑞𝑧²=1−𝑧²2^𝑟𝑄𝑧²,𝑙𝑧²=1+𝑧2^𝑟𝐿(𝑧),𝐴₀(𝑧)=𝑧^𝜈−𝑟𝐿𝑧⁻¹𝑆𝑧⁻²,𝐵₀(𝑧)=𝐿(−𝑧)1−𝑧2^2𝑟𝑄𝑧². (57) Now, suppose that 𝑃₀(𝑧₁,𝑧₂) satisfies (41) with accuracy 𝑟+1. By Lemma 8, 𝑧=−1 is a root of multiplicity 2𝑟 of 𝐺(𝑧)∶=𝑧^𝜈𝐴(𝑧⁻¹)𝐵(𝑧). Since 𝐺(−𝑧)=𝐺(𝑧) by (41), 𝑧=1 is another root of multiplicity 2𝑟 of 𝐺(𝑧). But 𝐴(1)≠0, so 𝑧=1 is a root of multiplicity 2𝑟 of 𝐵(𝑧). Hence, we have 𝐴𝑧⁻¹=1+𝑧2^𝑟𝐴₁𝑧⁻¹,𝐵(𝑧)=1+𝑧2^𝑟1−𝑧2^2𝑟𝐵₁(𝑧). (58) By substituting (58) in (41), we have 𝐴₁𝑧⁻¹𝐵₁(𝑧)+𝐴₁−𝑧⁻¹𝐵₁(−𝑧)=0. (59) Then, by applying Lemma 5 to 𝐴₁(𝑧) and 𝐵₁(𝑧), we get that there exist an odd integer 𝜈≥deg𝐴(𝑧) and real polynomials 𝑆(𝑧), 𝑄(𝑧), and 𝐿(𝑧), such that (34) hold.
Since 𝐵(1)=0, and therefore 𝐴(1)=1, we can normalize 𝑆(𝑧) and 𝐿(𝑧), so that 𝑆(1)=𝐿(1)=1. And 𝑄(1)=(−1)^𝑟𝐿(−1) follows from condition (54). Hence 𝐴(𝑧) and 𝐵(𝑧) satisfy (34), and (37).

Theorem 10 (accuracy). Let the polynomials 𝐴(𝑧), 𝐵(𝑧), and 𝑃₀(𝑧₁,𝑧₂) be given by (34)–(38) and (30) with 𝑚₀(𝑧₂) according to (28). Then 𝑃₀(𝑧₁,𝑧₂) is a low-pass filter which can generate bivariate nonseparable compactly supported orthogonal scaling function with accuracy 𝑟+1.

Proof. This theorem follows from Lemmas 8, 9, and Theorem 7.

Remark 11. (i) According to the proof of Theorem 2.1 of [13], if 𝐴(𝑧), 𝐵(𝑧) satisfy (40), the polynomials 𝑆(𝑧) and 𝐿(𝑧) in Lemma 9 must satisfy (35) and (36).
(ii) The existence of 𝑆(𝑧) and 𝐿(𝑧) satisfying (35)–(38) is ensured by Corollary 4.2 of [13].

3.3. Smoothness

In this section, we measure the smoothness of the scaling function 𝜑(𝐱) generated by 𝑃₀(𝑧₁,𝑧₂) by its Hölder exponent (see [13, 23]). Let 𝑠=𝑛+𝑡, where 𝑛 is a nonnegative integer and 0≤𝑡<1. For nonnegative integers 𝛼₁,𝛼₂, and 𝛼=𝛼₁+𝛼₂, 𝐷^𝛼𝜑(𝐱) denotes the partial derivative of 𝜑(𝐱); that is, 𝐷^𝛼𝜑(𝐱)𝜕∶=^𝛼𝜑(𝐱)𝜕^𝛼1𝑥₁𝜕^𝛼2𝑥₂,𝐱=𝑥₁,𝑥₂∈ℝ². (60) If there exists a constant 𝑐, such that 𝐷^𝛼𝜑(𝐱)−𝐷^𝛼𝜑(𝐲)≤𝑐|𝐱−𝐲|^𝑡,∀𝛼≤𝑛and𝐱,𝐲∈ℝ², (61) we say that 𝜑(𝐱)∈𝐶^𝑠(ℝ²). It is well known that if the Fourier transform 𝜑(𝜔) of 𝜑(𝐱) satisfies 𝜑(𝜔)≤𝑐1+𝜉^{−𝑠−1−𝜀}1+𝜂^{−𝑠−1−𝜀},𝜔=𝜉,𝜂∈ℝ², (62) for some constants 𝑐>0, 𝜀>0, then 𝜑(𝐱)∈𝐶^𝑠(ℝ²) (see [13, 23]).

Let us assume that the polynomials 𝐴(𝑧) and 𝐵(𝑧) satisfying (34)–(38) are 𝐴(𝑧)=1+𝑧2^𝑟𝐴₀(𝑧),𝐵(𝑧)=1+𝑧2^𝑟𝐵₀(𝑧), (63) and let 𝑙(𝑧)=𝐴₀(𝑧)+𝐵₀(𝑧). (64)

Define a polynomial 𝑇_𝑟(𝑦) as follows: 𝑇_𝑟𝑦=^𝑟−1_𝑘=0𝑘𝑦𝑟+𝑘−1_𝑘,for𝑦∈ℝ. (65)

Lemma 12. Let 𝑙(𝑧) be given by (64). Then there exist constants 𝑐>0 and 0<𝜀<1, such that ^∞_𝑘=1𝑙𝑧^2−𝑘≤𝑐𝜉^{𝜇𝑟+1−𝜀}, (66) where 𝑧=𝑒^𝑖𝜉 and 𝜇_𝑟∶=(1/2)log₂𝑇_𝑟(3/4). Furthermore,it is 𝑟−𝜇_𝑟>11−2log₂31𝑟+2log₂3,lim_𝑟→∞𝑟−𝜇_𝑟−2𝑟1=1−2log₂3≈0.2075. (67)

Proof. See the proof of Lemma 4.1, 4.3–4.5 of [13].

Theorem 13 (smoothness). Let the polynomials 𝑚₀(𝑧), 𝐴(𝑧), and 𝐵(𝑧) be given by (28) and (34)–(38), respectively. If 𝜑(𝐱) is the bivariate scaling function generated by the low-pass filter 𝑃₀(𝑧₁,𝑧₂) according to (30) and 𝑟≥5, then 𝜑(𝐱)∈𝐶^{𝑟−𝜇𝑟−2}(ℝ²).

Proof. We only need to prove that the Fourier transform 𝜑(𝜔) of 𝜑(𝐱) satisfies (62). In fact, by applying Lemmas 7.1.1–7.1.8 of [2], we have ^∞_𝑘=1𝑚₀𝑒^{𝑖𝜂/2𝑘}≤𝑐1+𝜂^{−𝑟+𝜇𝑟}≤𝑐1+𝜂^{−𝑟+𝜇𝑟+1−𝜀}, (68) for the same 𝜀 as Lemma 12, and ^∞_𝑘=11+𝑒^{𝑖𝜉/2𝑘}2^𝑟≤𝑐1+𝜉^−𝑟. (69) Thus, according to (63), (64), (68), and (69) and Lemma 12, we get ^∞_𝑘=1𝑃𝑒^{𝑖𝜉/2𝑘},𝑒^{𝑖𝜂/2𝑘}≤^∞_𝑘=1𝑚₀𝑒𝑖𝜂2^𝑘∞_𝑘=11+𝑒𝑖𝜉2^𝑘2^𝑟∞_𝑘=1𝑙𝑒^{𝑖𝜉/2𝑘}≤𝑐1+𝜉^{−𝑟+𝜇𝑟+1−𝜀}1+𝜂^{−𝑟+𝜇𝑟+1−𝜀}. (70) Hence by using (17) it is shown that (62) holds for 𝑠=𝑟−𝜇_𝑟−2.

Remark 14. (i) The wavelets have the same smoothness as the corresponding scaling function. So we only need to study the smoothness of the scaling function.
(ii) As the remark of paper [13], (67) mean that the smoothness of the bivariate orthogonal wavelets improves asymptotically by 0.2075 when 𝑟 is incremented by 1. Hence, (67) guarantees the existence of bivariate orthogonal wavelets of any desired smoothness. In particular, if 𝑟>4, then 𝜑(𝐱)∈𝐶⁰(ℝ²), if 𝑟>8, then 𝜑(𝐱)∈𝐶¹(ℝ²), and if 𝑟>12, then 𝜑(𝐱)∈𝐶²(ℝ²).
According to the results of this section, we can summarize the following procedure for constructing the filter banks.

Step 1. Choose a nonnegative integer 𝑟, and let 𝜈=4𝑟−1.

Step 2. Use some spectral factorization method to make a list of all polynomials 𝑆(𝑧) and 𝐿(𝑧) that satisfy (35)–(38). Substitute 𝑆(𝑧),𝐿(𝑧), and 𝑄 in (34) to get 𝐴(𝑧) and 𝐵(𝑧).

Step 3. Construct the low-pass filter 𝑃₀(𝑧₁,𝑧₂) by substituting 𝐴(𝑧),𝐵(𝑧), and univariate Daubechies low-pass filter 𝑚₀(𝑧₂) in (30).

Step 4. Construct the high-pass filters 𝑃_𝑖(𝑧₁,𝑧₂),𝑖=1,2,3 by substituting 𝐴(𝑧),𝐵(𝑧), and univariate Daubechies filters 𝑚₀(𝑧₂) and 𝑚₁(𝑧₂) in (31)–(33).

3.4. Example

Finally, we give an example of bivariate nonseparable compactly supported orthogonal low-pass filter by using our method.

Let 𝑟=2. The following polynomial 𝑚₀(𝑧) is the univariate Daubechies low-pass filter with accuracy 3 (see [27]): 𝑚₀(𝑧)=⁵_𝑘=0𝑝_𝑘𝑧^𝑘, (71) where 𝑝₀=0.332670552950,𝑝₁𝑝=0.806891509311,₂=0.459877502118,𝑝₃𝑝=−0.135011020010,₄=−0.085441273882,𝑝₅=0.035226291886. (72)

If (38) holds, then there are four polynomials 𝑆(𝑧) satisfying (35), (37) and two polynomials 𝐿(𝑧) satisfying (36), (37). Hence, one can get 8 kinds of 𝐴(𝑧) and 𝐵(𝑧) satisfying (34)–(38) (see [13]). One of them is given as follows: 𝐴(𝑧)=⁷_𝑘=0𝑎_𝑘𝑧^𝑘,𝐵(𝑧)=⁷_𝑘=0𝑏_𝑘𝑧^𝑘, (73) where the coefficients 𝑎_𝑘,𝑏_𝑘,𝑘=0,1,…,7 are displayed in Table 1.

Table 1 

The coefficients of the polynomials 𝐴(𝑧) and 𝐵(𝑧) of the example given in Section 3.4

Then, by using (30), we can get a bivariate nonseparable compactly supported orthogonal low-pass filter (Figure 1) given as follows: 𝑃₀𝑧₁,𝑧₂=⁷_𝑘,𝑙=0𝑝_𝑘,𝑙𝑧^𝑘₁𝑧^𝑙₂, (74) where the coefficients 𝑝_𝑘,𝑙 are displayed in Table 2.

Table 2 

The coefficients of the low-pass filter 𝑃₀(𝑧₁,𝑧₂) of the example given in Section 3.4

(a)

(b)

Figure 1 

Plot of the real (a) and imaginary (b) parts of the bivariate compactly supported low-pass filter (74).

For instance, by explicit computation and the approximation of coefficients the first few terms of the filter look like 0.002629−0.004553𝑒^(𝑖/2)𝜉−0.045936𝑒^𝑖𝜉+0.040313𝑒^{(3𝑖/2)𝜉}+0.182944𝑒^2𝑖𝜉+0.115161𝑒^{(5𝑖/2)𝜉}+0.026698𝑒^3𝑖𝜉+0.015413𝑒^{(7𝑖/2)𝜉}−0.011045𝑒^{(𝑖/2)𝜉+(𝑖/2)𝜂}−0.111418𝑒^{𝑖𝜉+(𝑖/2)𝜂}+0.09778𝑒^{(3𝑖/2)𝜉+(𝑖/2)𝜂}+0.44373𝑒^{2𝑖𝜉+(𝑖/2)𝜂}⋯. (75) By the inverse Fourier transform we finally obtain the scaling function, so that we have a sequence in the following form: ()𝛿()𝛿()𝛿(𝑥)𝛿()𝛿()𝛿()𝛿()𝛿−0.00921𝛿−3+𝑥−3+𝑦−0.001424𝛿−2+𝑥−3+𝑦+0.030481𝛿−1+𝑥−3+𝑦−0.019846𝛿−3+𝑦−0.010634𝛿−7+2𝑥−3+𝑦+0.090022𝛿−5+2𝑥−3+𝑦−0.148136𝛿−3+2𝑥−3+𝑦+0.06875𝛿−1+2𝑥−3+𝑦+⋯. (76)

The scaling function generated by the low-pass filter (Figure 2) has accuracy 3. By (31)–(33), one can easily obtain the corresponding high-pass filters. We omit it here.

Figure 2 

The scaling function generated by the low-pass filter (74).

4. Conclusion

We gave a method for constructing bivariate nonseparable compactly supported orthogonal wavelets by combining three methods using the dilation matrix 2002. We studied the accuracy of the scaling functions and the smoothness of the wavelets constructed by this method. We proved that the scaling functions have the same accuracy order as the univariate Daubechies low-pass filter used in this method and the wavelets can be made arbitrarily smooth by choosing the accuracy order.

Acknowledgments

The authors wish to thank the anonymous reviewers for their valuable comments and suggestions that have improved the presentation of this paper. This work was supported by the National Natural Science Foundation of China (11271001, 61170311), 973 Program (2013CB329404), and Sichuan Province Science and Technology Research Project (12ZC1802).