Matrix Spectral Factorization for SA4 Multiwavelet

10/15/2019 ∙ by Vasil Kolev, et al. ∙ 0

In this paper, we investigate Bauer's method for the matrix spectral factorization of an r-channel matrix product filter which is a halfband autocorrelation matrix. We regularize the resulting matrix spectral factors by an averaging approach and by multiplication by a unitary matrix. This leads to the approximate and exact orthogonal SA4 multiscaling functions. We also find the corresponding orthogonal multiwavelet functions, based on the QR decomposition.



There are no comments yet.


page 23

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

If is a matrix polynomial which represents a causal multifilter bank, the product filter bank is given by , where . The coefficients are matrices. Matrix Spectral Factorization (MSF) describes the problem of finding , given .

MSF is still quite unknown in the signal processing community. Because of its numerous possible applications, in particular in multiple-input and multiple-output (MIMO) communications, image processing, multidimensional control theory, and others, it deserves to receive greater attention. We are interested in MSF as a tool for designing orthogonal multiwavelets.

It is known that MSF is possible as long as is positive semidefinite for in the unit circle in the complex plane Ephremidze et al. (2009); Hardin et al. (2004). We call degenerate if it is not strictly positive definite, that is, at one or more points.

A number of numerical approaches to MSF have been proposed, but they usually cannot handle the degenerate case. Unfortunately, those are exactly the cases of greatest interest for applications in the construction of multiwavelets.

We use Bauer’s method, as described by Youla and Kazanjian Youla and Kazanjian (1978), which is based on the Cholesky factorization of a block Toeplitz matrix. This algorithm can handle the degenerate case, but convergence is quite slow.

We study the SA4 multiwavelet as a test case, and as a benchmark for future applications. Starting from a known , we compute and factor it again.

It is easy to see that matrix spectral factorization is not unique. If is any factor, then

is also a factor for any orthogonal matrix

, or more generally for paraunitary . For the factor produced by Bauer’s method, the matrix is always lower triangular, so we are forced to use further modifications to recover the original SA4.

In this paper, we only concentrate on modifications which speed up convergence and produce factors close to original SA4. This leads to two factorizations: the approximate SA4 multifilter, which is close to SA4, and the exact SA4 multifilter, which is equal to SA4 except for insignificant roundoff. In future papers, we plan to impose desirable conditions directly on the spectral factor, without working towards a known result.

The main contributions of this paper are

  • Considering a novel method for computing the multichannel spectral factorization of degenerate matrix polynomial matrices;

  • Speeding up the slow convergence of the algorithm;

  • Minimizing the numerical errors which appear in the algorithm;

  • Demonstrating by example that this approach can be used to construct an orthogonal symmetric multiwavelet filter;

  • Comparing the frequency responses and experimental numerical performance of the approximate and exact SA4 multiwavelets.

The rest of the paper is organized as follows. In chapter 2 we give some background material, state the problem, and describe related research. We also introduce matrix product filter theory, and present the product filter of the SA4 multiscaling function as a benchmark test case. In chapter 3 we consider Bauer’s method, and discuss its numerical behavior. In chapter 4 we describe the approach for obtaining the approximate and exact SA4 multiscaling and multiwavelet functions. The performance analysis of these multiwavelets is shown in chapter 5. Chapter 6 gives conclusions.

2 Background and Problem Statement

We will use the following notation conventions, illustrated with the letter ‘a’:

– lowercase letters refer to scalars or scalar-valued functions;

– lowercase bold letters are vectors or vector-valued functions;

– uppercase letters are matrices or matrix-valued functions.

The symbols and denote the identity and zero matrices of appropriate size, respectively.

We will be using polynomials in a complex variable on the unit circle, but all coefficients will be real-valued. Thus, the complex conjugate of a matrix polynomial is given by

2.1 Matrix Product Filters

A two-channel multiple-input multiple-output (MIMO) filter bank of multiplicity is shown in fig. 1. Here is the input signal vector of length , and the analysis multifilters are

where , are matrices of size . Likewise, and are the synthesis multifilters, and is the output vector signal. For simplicity, we assume that all coefficients are real.

The input-output relation of this filter bank is

It is convenient to write this equation in matrix form as

Figure 1: Two-channel multifilter bank.

In general, the design of the multifilter bank requires four multifilters: two on the analysis and two on the synthesis side. In perfect-reconstruction orthogonal MIMO filter banks, the analysis modulation matrix

is paraunitary, that is


In this case we can define the synthesis filters in terms of the analysis filters:

Eq. (1) can also be expressed as

or in terms of the coefficients

for any integer .

Given , the matrix polynomial

is called the matrix lowpass product filter (or scalar product filter in the case ). The coefficients satisfy . is a half-band polynomial filter Cooklev et al. (1996), that is, it satisfies the equation


This implies that , and for all nonzero integers .

2.2 Multiwavelet Theory

2.2.1 Multiscaling and Multiwavelet Functions

Consider the iteration of a MIMO filter bank along the channel containing . After iterations, the equivalent filters will be

or in the time domain

where , .

The filterbank coefficients are associated with function vectors , called the multiscaling function, and , called the multiwavelet function. These functions satisfy the recursion equations

The support of , is contained in the interval , but it could be strictly smaller than the interval. See Massopust et al. (1996) for more details.

2.2.2 Multiwavelet Properties

Pre- and Postfiltering

Since multifilters are MIMO, they operate on several streams of input data rather than one. Therefore, input data needs to be vectorized. There are many prefilters available for various applications, but three kinds are preferred:

  • Oversampling - Repeated Row Preprocessing

  • Critical Sampling - Matrix Preprocessing

  • Embedded Orthogonal Symmetric Prefilter Bank

In Strela and Walden (2000), the first and second approach is presented for the GHM and CL multiwavelets. In practical applications, a popular choice is based on the Haar transform Cotronei et al. (1998); Tan et al. (1999). Haar pre- and postfilters have the advantage of simultaneously possessing symmetry and the orthogonality, and no multiplication is needed if one ignores the scaling factor. The problem of constructing symmetric prefilters is considered in detail in Hsung et al. (2003), with proposed solutions for the DGHM and CL multifilters.

The third technique is applied in Gan and Ma (2005); Hsung et al. (2007)

. They find a three-tap prefilter for the DGHM multifilter by searching among the parameters which minimize the first-order absolute moment of the filter coefficients.

Similarly, at the filter bank output, a postfilter is needed to reconstruct the signal. Pre- and postprocessing is not needed for scalar wavelets.

Balancing Order

The orthogonal multiscaling function is said to be balanced of order if the signals are preserved by the operator for . That is,


Multiwavelets that do not satisfy this property are said to be unbalanced. See Lebrun and Vetterli (2001, 1998); Plonka and Strela (1998); Li and Peng (2011) for more details. Balanced multiwavelets do not require preprocessing.

Approximation Order

In the scalar case, a certain approximation order refers to the ability of the low-pass filter to reproduce discrete-time polynomials, while the wavelet annihilates discrete-time polynomials up to the same degree. Since many real-world signals can be modeled as polynomials, this property typically leads to higher coding gain (CG).

In the case of non-balanced multiwavelets, appropriate preprocessing is required to take advantage of approximation orders.

The approximation order of multiscaling and multiwavelet functions can be determined using the following result established in Wu et al. (2010). A multiscaling function provides approximation order iff there exist vectors , , , which satisfy

where is the derivative of order , and the masks of multiscaling and multiwavelet functions are


To fully characterize the multifilter bank with respect to approximation order we can add the condition

Good Multifilter Properties (GMPs)

A multiwavelet system can be represented by an equivalent scalar filter bank system Tham et al. (2000). The equivalent scalar filters are, in fact, the polyphases of the corresponding multifilter. This relationship motivates a new measure called good multifilter properties (GMPs) Tham et al. (1998), (Tham et al., 2000, (Def. 1))

, which characterizes the magnitude response of the equivalent scalar filter bank associated with a multifilter. GMPs provide a set of design criteria imposed on the scalar filters, which can be translated directly to eigenvector properties of the designed multiwavelet filters.

An orthogonal multiwavelet system has GMPs of order at least if the following conditions are satisfied Tham et al. (2000)

where and are masks of the lowpass and highpass filters (see eq. (3)), and

A class of symmetric-antisymmetric biorthogonal multiwavelet filters which possess GMPs was introduced in Tan et al. (1999).

A GMP order of at least is critical for ensuring no frequency leakage across bands, hence improving compression performance. Multifilters possessing GMPs have better performance than those which do not possess GMPs Li and Yang (2010), due to better frequency responses for energy compaction, greater regularity and greater approximation order of the corresponding wavelet/scaling functions. GMPs help prevent both DC and high-frequency leakage across bands; this contributes to reduced smearing, blocking, ringing artifacts, and also helps to prevent checkerboard artifacts in reconstructed images for image coding.

The SA4 multiwavelet, introduced below in section 3.3, is a member of a one-parameter family of orthogonal multiwavelets with four coefficient matrices Tham et al. (2000). Different members of the family have different approximation order, but all of them have a GMP order of at least . This manifests itself in the smooth decay to zero of the magnitude responses near , compared to the following 4-tap orthogonal multiwavelet filters:

  • The GHM multiwavelet Donovan et al. (1996) has symmetric orthogonal scaling functions and an approximation order of 2 for its filter length 4.

  • Chui and Lian’s CL multiwavelet Chui and Lian (1996) has the highest possible approximation order of 3 for its filter length 3.

  • Jiang’s multiwavelet JOPT4 Jiang (1998b) has optimal time-frequency localization for its filter length.

Although the GHM and CL systems are the most commonly used orthogonal multiwavelet systems, and have higher approximation order than the SA4 system, they do not satisfy GMPs.


The smoothness of , can be characterized by Sobolev regularity ; this is discussed in Micchelli and Sauer (1997); Cohen et al. (1997); Jiang (1998a)

. The regularity estimate is related to both the eigenvalues and the corresponding eigenvectors of the transition operator and to spectral radius of the transition operator.

Symmetry and Antisymmetry

A function is symmetric about a point if for all . It is antisymmetric if .

For scaling functions and wavelets of compact support, the only possible symmetry is about the center of their support. For multiscaling and multiwavelet functions, this could be a different point for different components.

In the scalar case, the scaling function cannot be antisymmetric, since it must have a nonzero integral. In the multiwavelet case, some components even of can have antisymmetry.

2.3 Problem Statement

The design of perfect reconstruction multifilter banks and multiwavelets remains a significant problem. In the scalar case, spectral factorization of a half-band filter that is positive definite on the unit circle was, in fact, the first design technique, suggested by Smith and Barnwell Smith and Barnwell (1986). This provides a motivation to extend the technique to the multiwavelet case. Spectral factorization of a Hermitian product filter in the multiwavelet case is more challenging than in the scalar case, and has not been investigated previously.

The present paper considers the application of MSF to the problem of finding an orthogonal lowpass multifilter , given a product filter . Since every orthogonal multiscaling filter is a spectral factor of some matrix product filter, this can be used as a tool for designing new filters.

With an eye towards future applications, we are interested in constructing multifilters with desirable properties, especially GMPs. This implies that the determinant of will have at least one zero of higher multiplicity on the complex unit circle. That is, will be degenerate.

As a test case, we want to use a product filter which satisfies the half-band condition (2), and is derived from an which has more than 2 taps (i.e. ) and good regularity properties and GMPs. Our benchmark test case will be the 2-channel orthogonal SA4 multiwavelet.

While there are many matrix spectral factorization algorithms Cooklev et al. (1996); Ježek and Kučera (1985); Lawton (1993), most cannot handle the degenerate case. We use Bauer’s method, based on the Toeplitz method of spectral factorization of the Youla and Kazanjian algorithm Bauer (1956); Youla and Kazanjian (1978), which can handle this case. However, convergence becomes very slow.

It is easy to see that matrix spectral factorization is not unique. If is a factor of a given , then is also a factor for any orthogonal matrix , or more generally for paraunitary . Other operations may be permissible in some cases.

In our experiments with Bauer’s algorithm, the initial factors do not correspond to smooth functions. We apply regularization techniques, both to speed up convergence and to work towards recovering the original filter in this test case. In this manner, we derive two filters which we call the approximate and exact SA4 multiwavelets.

Let us briefly summarize some previous work.

Matrix Spectral Factorization (MSF) plays a crucial role in the solution of various applied problems for MIMO systems in communications and control engineering Wang et al. (2015). It has been applied to designing minimum phase FIR filters and the associated all-pass filter Hansen et al. (2010), quadrature-mirror filter banks Charoenlarpnopparut (2007), MIMO systems for optimum transmission and reception filter matrices for precoding and equalization Fischer (2005), precoders Du et al. (2013), and many other applications.

The analysis of linear systems corresponding to a given spectral density function was first established by Wiener, who used linear prediction theory of multi-dimensional stochastic processes Wiener and Masani (1957). The method of Wilson for MSF was developed for applications in the prediction theory of multivariate discrete time series, with known (or estimated) spectral density Wilson (1972). Many numerical approaches to MSF have been proposed, for example Cooklev et al. (1996); Ježek and Kučera (1985); Lawton (1993). A survey of such algorithms is given in Sayed and Kailath (2001); however, this does not include algorithms for the degenerate case.

It is well known that all MSF methods have difficulties in this situation, and some of them cannot handle zeros on the unit circle at all. For example, Kučera’s algorithm, an otherwise popular matrix spectral factorization algorithm, has this limitation, and is not suitable for our purpose.

Bauer’s method was successful applied in Roux (1986) to the Radon projection. In that paper, factorization of the autocorrelation matrix of the Radon projection of a minimum phase pseudo-polynomial restored the coefficients of the original pseudo-polynomial with an accuracy around after 10,000 recursions on a 64-bit floating point processor.

3 Matrix Spectral Factorization

The following theorem answers the existence question for MSF.

Theorem 3.1

(Matrix form of the Riesz-Fejér lemma Hardin et al. (2004); Ephremidze et al. (2009)) A matrix polynomial

can be factored as , where

is an -channel causal polynomial matrix, if and only if is symmetric positive semidefinite for on the complex unit circle.

3.1 Bauer’s Method

If and are known, we can construct doubly infinite block Toeplitz matrices and from their coefficients, by setting , . is symmetric and block banded with bandwidth ; is block lower triangular, also with bandwidth .

The relation corresponds to , that is, is a Cholesky factor of .

In Bauer’s method, we pick a large enough integer , and truncate to (block) size :

If is positive definite, we can compute the Cholesky factorization , where

It is shown in Youla and Kazanjian (1978) that this factorization is possible, and that as .

Bauer’s method works even in the case of highly degenerate . Unfortunately, convergence in this case is very slow.

3.2 Numerical Behavior

For chosen size , the matrix is of size

. Its singular values are all in the range

. The first singular values are close to 2, then , and the remaining singular values are close to 0. See fig. 2.

Figure 2: A sharp drop of the singular values of the Toeplitz matrix and of its spectral factor , for ; (filled circles) singular values of ; (hollow circles) singular values of .

The numerical stability of the factorization algorithm is good, until the higher singular values get too close to 0.

A bigger problem is the slow convergence. For example, in the scalar case , where the factor is , the value on the diagonal in row is

which converges only very slowly to the limit 1. This simple problem has a double zero of the determinant on the unit circle. For higher order zeros, convergence speed becomes even worse.

3.3 Testing Bauer’s Method

A flowchart for testing Bauer’s method for MSF is shown in fig. 3. The steps are

  • Choice of the benchmark multiscaling function ;

  • Construction of the matrix lowpass product filter ;

  • Find a spectral factor , and assess its quality;

  • Do different kinds of postprocessing to obtain spectral factors , ;

  • Compare the obtained factors with the benchmark .

Figure 3: Steps in testing matrix spectral factorization using Bauer’s method.

Step 1: We consider the well-known orthogonal SA4 multiwavelet Chui and Lian (1996); Huo and Miao (2012); Tham et al. (2000) as a benchmark testing case.

The recursion coefficients and of SA4 depend on a parameter . They are


where . In this paper, we use the filter with , which leads to . Numerical values of , are listed in table 1.

These coefficients generate the 2-band compactly supported orthogonal multiscaling function and multiwavelet function , all with support . Graphs of these functions can be found in fig. 8 in a later chapter.

Table 1: Coefficients of the original SA4 multiscaling and multiwavelet functions.

Step 2: The symbol of the product lowpass multifilter has the form


Numerical values of are given in table 2.

0 1 0
0 1
Table 2: The matrix coefficients of the half-band product filter

The impulse responses of the product filter (see fig. 4) have the character of nearly Haar-type scaling and wavelet functions. The frequency responses have good selectivity.

Figure 4: The impulse response (left) and frequency response (right) of the two components of the product filter .

Step 3: The matrix in this case looks like this:

For any fixed , the coefficients from the last row of are approximate factors of . Suppressing the dependence on in the notation, we define the filter

and the matrix


This notation will be used again later.

To measure the quality of the factorization, we compute

The residual is then

In numerical experiments, the minimal error is achieved for ; the coefficients are given in table 3. The numerical error for Bauer’s method is much better than that of Wilson’s method for MSF, which cannot be lower then Ephremidze et al. . The impulse response of the obtained spectral factor is non-regular, but the frequency response has good selectivity. This approximate factor does not quite provide perfect reconstruction. This is addressed in more detail in subsection 5.2.

0.094428373297668 0
-0.091754813647953 0.022310801334572
0.175943193428843 0.679731045265413
0.010360133828403 -0.702056243347588
0.000129414745433 0.700756678587017
0.165030422671601 0.681046913905671
-0.081190837390599 0.021002135513698
-0.083853536287580 -0.001275235049147
Table 3: Coefficients of the lowpass spectral factor for , with .

Steps 4 and 5 will be covered in chapter 4.

4 Construction of Approximate and Exact SA4 Multiwavelets

Above, we showed how to obtain the approximate spectral factor for a given size . Given any factor of , other possible factors can be found be postmultiplication with suitable orthogonal matrices, or in some cases by other manipulations.

In subsection 4.1 we use a simple averaging process and column reversal to produce a lowpass filter that is similar to, but not identical to, the original SA4 multiscaling function. We call this the approximate SA4 multiwavelet.

In subsection 4.2, we add a rotation postfactor and an averaging process to produce another multiwavelet which is even closer to the original SA4 multiwavelet. We call this the exact SA4 multiwavelet.

The newly derived multiscaling filters are denoted by

for the approximate SA4 multiwavelet, and

for the exact SA4 multiwavelet. ( is an intermediate step in the calculation of ). In each case, we find the corresponding multiwavelet filter by using a QR decomposition Golub and Van Loan (2013). The multiscaling filter is factorized as

The coefficients are then obtained from the third and fourth columns of .

We define the absolute error by

where for , and likewise for .

To measure the deviation of the new filters from the original, we introduce the mean square error (MSE) of the multiscaling and multiwavelet functions,

and themaximal absolute errors (MAE) of the matrix coefficients, the multiscaling function and the multiwavelet function as


4.1 Approximate SA4 Multiwavelet

The spectral factor obtained from Bauer’s method leads to non-symmetrical and non-regular functions. A simple algorithm can be used to symmetrize and regularize the scaling functions. Below, we are using the notation from eq. (5).

First, we average the absolute values of the first and fourth, as well as the second and third, matrix coefficients of the spectral factor:

and construct average matrices

Second, the matrix coefficients and are obtained by multiplying the averaged matrices from the right by , but keep the signs, while the matrix coefficients and are obtained by multiplying and on the left and right by :

This produces coefficients quite close to the original coefficients in (4).

In order to determine the best approximate solution, we investigate the influence of the size of the Toeplitz matrix for the numerical errors MAE-MF, MSE-MC, and MSE-MC. These errors are shown in fig. 5 in logarithmic scale, for up to 100. Note that MSE-MC = MSE-MC, and MSE-MC = MSE-MC, by construction.

Figure 5: Dependence of maximum absolute error MAE-MF and mean squared errors MSE-MC = MSE-MC, MSE-MC = MSE-MC on matrix size (log. scale).

The important minimal errors MAE-MF and MSE-MwF for are tabulated in table 4; the matrix coefficients of are listed in table 5. Obviously, the development of the errors shows the convergence of the algorithm, as well as the dependency of the errors on the size . It also shows that the main error comes from the matrix coefficients and , rather than and .

The error MSE-MF behaves essentially the same way as MAE-MC. Unfortunately, the convergence of MAE-MC and MAE-MC stalls at a relatively large number. Thus, the approximate SA4 multiwavelet is not very accurate. Better matrix coefficients are achieved in the next subchapter.

Table 4: The minimum MAEs and MSEs for the approximate and exact SA4 multiscaling and multiwavelet functions, the size and the corresponding angle .
Table 5: Coefficients of the approximate SA4 multiscaling and multiwavelet functions for .

4.2 Exact SA4 Multiwavelet

The second algorithm leads to more regular scaling functions. We multiply the coefficients of from the right by the unitary matrix

This produces the approximate factor .

The angle (in radians) is calculated as the average value



Again, we are using the notation from eq. (5).

Figure 6 shows the dependence of the angle on the size (log. scale). It can be seen that in the angle there is one overshoot, after that it is decreasing very slowly. The minimum errors are obtained at the maximal possible size of Toeplitz matrix , with the values shown in table 4. This shows the influence of the very slow convergence; despite right multiplication with the unitary matrix, desirable precison cannot be an achieved. The errors MSE-MF and MAE-MF for are shown in fig. 7.

To increase the numerical precision and obtain the exact SA4 multiscaling function, we again apply an averaging approach.

The errors for are shown in fig. 7, with minima shown in table 4. Let us consider the accuracy of the multiscaling functions obtained at the leading local minima of the error MAE-MF