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
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.
. 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
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.
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.
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.
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.
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 .
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.
Step 2: The symbol of the product lowpass multifilter has the form
Numerical values of are given in table 2.
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.
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.
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.
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.
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.