Matrix Spectral Factorization for SA4 Multiwavelet

Multidimensional Systems and Signal Processing manuscript No. (will be inserted by the editor) Matrix Spectral F actorization - SA4 Multiwa velet V asil Kole v · T odor Cooklev · Fritz Keinert Receiv ed: date / Accepted: date Abstract In this paper, we in vestigate Bauer’ s method for the matrix spectral factorization of an r -channel matrix product filter which is a half-band autocorrelation matrix. W e regu- larize 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. W e also find the corresponding orthogonal multiwav elet functions, based on the QR decomposition. Keyw ords SA4 orthogonal multiwav elet · matrix spectral factorization 1 Introduction If H ( z ) = C 0 + C 1 z − 1 + · · · C n z − n is a matrix polynomial which represents a causal multi- filter bank, the product filter bank is gi ven by P ( z ) = H ( z ) H ∗ ( z ) , where H ∗ ( z ) = H T ( z − 1 ) . The coefficients C k are r × r matrices. Matrix Spectral F actorization (MSF) describes the problem of finding H ( z ) , given P ( z ) . V asil K olev Institute of Information and Communication T echnologies Bulgarian Academy of Sciences Bl. 2 Acad. G. Bonchev St. 1113 Sofia, Bulgaria E-mail: kole v acad@abv .bg T odor Cooklev , Senior Member IEEE W ireless T echnology Center Indiana Univ ersity – Purdue Univ ersity Fort W ayne, IN 46805, USA E-mail: cooklevt@ipfw .edu Fritz Keinert Dept. of Mathematics Iow a State Univ ersity Ames, IA 50011, USA E-mail: keinert@iastate.edu 2 V asil K olev et al. MSF is still quite unkno wn in the signal processing community . Because of its numerous possible applications, in particular in multiple-input and multiple-output (MIMO) communi- cations, image processing, multidimensional control theory , and others, it deserves to receive greater attention. W e are interested in MSF as a tool for designing orthogonal multiwa velets. It is known that MSF is possible as long as P ( z ) is positive semidefinite for z in the unit circle in the complex plane [13, 20]. W e call P ( z ) deg enerate if it is not strictly positi v e definite, that is, det ( P ( z )) = 0 at one or more points. A number of numerical approaches to MSF ha ve 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 multiwa velets. W e use Bauer’ s method, as described by Y oula and Kazanjian [50], which is based on the Cholesky f actorization of a block T oeplitz matrix. This algorithm can handle the de generate case, but con v ergence is quite slo w . W e study the SA4 multiwavelet as a test case, and as a benchmark for future applications. Starting from a known H ( z ) , we compute P ( z ) and factor it again. It is easy to see that matrix spectral factorization is not unique. If H 1 ( z ) is any factor , then H 2 ( z ) = H 1 ( z ) U is also a f actor for any orthogonal matrix U , or more generally H 2 ( z ) = H 1 ( z ) U ( z ) for paraunitary U ( z ) . For the f actor produced by Bauer’ s method, the matrix C 0 is always lower triangular , so we are forced to use further modifications to reco ver the original SA4. In this paper , we only concentrate on modifications which speed up con vergence 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 tow ards a known result. The main contributions of this paper are – Considering a nov el method for computing the multichannel spectral factorization of degenerate matrix polynomial matrices; – Speeding up the slow con vergence 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 multiwa velet filter; – Comparing the frequency responses and experimental numerical performance of the approximate and exact SA4 multiwa velets. The rest of the paper is organized as follows. In chapter 2 we give some background material, state the problem, and describe related research. W e 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 e xact SA4 multiscal- ing and multiwa velet functions. The performance analysis of these multiwavelets is shown in chapter 5. Chapter 6 giv es conclusions. 2 Background and Pr oblem Statement W e will use the following notation con ventions, illustrated with the letter ‘a’: a – lowercase letters refer to scalars or scalar -valued functions; a – lowercase bold letters are v ectors or vector-v alued functions; Matrix Spectral Factorization - SA4 Multiwa velet 3 Fig. 1 T wo-channel multifilter bank. A – uppercase letters are matrices or matrix-valued functions. The symbols I and 0 denote the identity and zero matrices of appropriate size, respectively . W e will be using polynomials in a complex variable z on the unit circle, but all co- efficients will be real-valued. Thus, the complex conjugate A ∗ ( z ) of a matrix polynomial A ( z ) = ∑ k A k z k is giv en by A ∗ ( z ) = A T ( z − 1 ) = ∑ k A T k z − k . 2.1 Matrix Product Filters A two-channel multiple-input multiple-output (MIMO) filter bank of multiplicity r is shown in fig. 1. Here x ( z ) is the input signal vector of length r , and the analysis multifilters are H ( z ) = C 0 + C 1 z − 1 + · · · C n z − n , G ( z ) = D 0 + D 1 z − 1 + · · · D n z − n , where C k , D k are matrices of size r × r . Likewise, E ( z ) and F ( z ) are the synthesis multifilters, and ˆ x ( z ) is the output vector signal. For simplicity , we assume that all coefficients are real. The input-output relation of this filter bank is ˆ x ( z ) = 1 2 [ E ( z ) H ( z ) + F ( z ) G ( z )] x ( z ) + 1 2 [ E ( z ) H ( − z ) + F ( z ) G ( − z )] x ( − z ) . It is con venient to write this equation in matrix form as  ˆ x ( z ) ˆ x ( − z )  = 1 2  E ( z ) F ( z ) E ( − z ) F ( − z )   H ( z ) H ( − z ) G ( z ) G ( − z )   x ( z ) x ( − z )  . In general, the design of the multifilter bank requires four multifilters: two on the anal- ysis and two on the synthesis side. In perfect-reconstruction orthogonal MIMO filter banks, the analysis modulation matrix M ( z ) =  H ( z ) H ( − z ) G ( z ) G ( − z )  . 4 V asil K olev et al. is paraunitary , that is M ∗ ( z ) M ( z ) = M ( z ) M ∗ ( z ) = I . (1) In this case we can define the synthesis filters in terms of the analysis filters: E ( z ) = H ∗ ( z ) , F ( z ) = G ∗ ( z ) . Eq. (1) can also be expressed as H ( z ) H ( z ) ∗ + H ( − z ) H ( − z ) ∗ = I , G ( z ) G ( z ) ∗ + G ( − z ) G ( − z ) ∗ = I , H ( z ) G ( z ) ∗ + H ( − z ) G ( − z ) ∗ = 0 , or in terms of the coefficients ∑ k C k C T k + 2 ` = ∑ k D k D T k + 2 ` = δ 0 ` I , ∑ k C k D T k + 2 ` = 0 , for any inte ger ` . Giv en H ( z ) , the matrix polynomial P ( z ) = H ( z ) H ∗ ( z ) = n ∑ k = 0 C k z − k ! n ∑ k = 0 C T k z k ! = n ∑ k = − n P k z k , is called the matrix lowpass pr oduct filter (or scalar product filter in the case r = 1). The coefficients P k satisfy P − k = P T k . P ( z ) is a half-band polynomial filter [7], that is, it satisfies the equation P ( z ) + P ( − z ) = 2 I . (2) This implies that P 0 = I , and P 2 ` = 0 for all nonzero inte gers ` . 2.2 Multiwa velet Theory 2.2.1 Multiscaling and Multiwavelet Functions Consider the iteration of a MIMO filter bank along the channel containing H ( z ) . After k iterations, the equiv alent filters will be H ( k ) ( z ) = k ∏ j = 0 H ( z 2 j ) = H ( z 2 k ) H ( z 2 k − 1 ) · · · H ( z ) , G ( k ) ( z ) = G ( z 2 k ) k − 1 ∏ j = 0 H ( z 2 j ) = G ( z 2 k ) H ( z 2 k − 1 ) · · · H ( z ) , or in the time domain C ( k ) j = ∑ j C ( k − 1 ) j C ( k − 1 ) j − 2 k − 1 , D ( k ) j = ∑ j D ( k − 1 ) j D ( k − 1 ) j − 2 k − 1 , Matrix Spectral Factorization - SA4 Multiwa velet 5 where C ( 0 ) j = C j , D ( 0 ) j = D j . The filterbank coef ficients are associated with function vectors ϕ ϕ ϕ = [ φ 0 , φ 1 , . . . , φ r − 1 ] T , called the multiscaling function , and ψ ψ ψ = [ ψ 0 , ψ 1 , . . . , ψ r − 1 ] T , called the multiwavelet func- tion . These functions satisfy the recursion equations ϕ ϕ ϕ ( t ) = √ 2 n ∑ k = 0 C k ϕ ϕ ϕ ( 2 t − k ) , ψ ψ ψ ( t ) = √ 2 n ∑ k = 0 D k ϕ ϕ ϕ ( 2 t − k ) . The support of ϕ ϕ ϕ , ψ ψ ψ is contained in the interval [ 0 , n ] , but it could be strictly smaller than the interval. See [32] for more details. 2.2.2 Multiwavelet Pr operties Pre- and P ostfiltering 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 av ailable for various applications, b ut three kinds are preferred: – Oversampling - Repeated Row Preprocessing – Critical Sampling - Matrix Preprocessing – Embedded Orthogonal Symmetric Prefilter Bank In [41], 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 [8, 42]. Haar pre- and postfilters have the adv antage 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 [22], with proposed solutions for the DGHM and CL multifilters. The third technique is applied in [16, 21]. They find a three-tap prefilter for the DGHM multifilter by searching among the parameters which minimize the first-order absolute mo- ment 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 wa velets. Balancing Order The orthogonal multiscaling function is said to be balanced of or der q if the signals u k = [ · · · , ( − 2 ) k , ( − 1 ) k , 0 k , 1 k , 2 k , · · · ] T are preserved by the operator L T for k = 0 , 1 , . . . , q − 1. That is, L T u k = 2 − k u k , where L =     · · · C 0 C 1 · · · C n C 0 C 1 · · · C n · · ·     . Multiwa velets that do not satisfy this property are said to be unbalanced . See [29, 28, 34, 30] for more details. Balanced multiwa velets do not require preprocessing. 6 V asil K olev et al. Appr oximation 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 wa velet annihilates discrete-time polyno- mials 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 multiwa velets, appropriate preprocessing is required to tak e advantage of approximation orders. The approximation order of multiscaling and multiwa velet functions can be determined using the following result established in [49]. A multiscaling function provides approxima- tion order p iff there e xist vectors u k ∈ R r , 0 ≤ k < p , u 0 6 = 0 , which satisfy k ∑ ` = 0 1 ( 2 i ) `  k `  u T k − ` D ` H ( 0 ) = 2 − k u T k , k ∑ ` = 0 1 ( 2 i ) `  k `  u T k − ` D ` H ( π ) = 0 , where D ` is the deri vati ve of order ` , and the masks of multiscaling and multiwa velet func- tions are H ( ω ) = 1 √ 2 n ∑ k = 0 C k e ik ω , G ( ω ) = 1 √ 2 n ∑ k = 0 D k e ik ω . (3) T o fully characterize the multifilter bank with respect to approximation order we can add the condition k ∑ ` = 0 1 ( 2 i ) `  k `  u T k − ` D ` G ( 0 ) = 0 . Good Multifilter Properties (GMPs) A multiwav elet system can be represented by an equiv alent scalar filter bank system [44]. The r equiv alent scalar filters are, in fact, the r polyphases of the corresponding multi- filter . This relationship motiv ates a new measure called good multifilter pr operties (GMPs) [43], [44, (Def. 1)], which characterizes the magnitude response of the equiv alent 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 eigen vector properties of the designed multiwa velet filters. An orthogonal multiwav elet system has GMPs of order at least ( 1 , 1 , 1 ) if the following conditions are satisfied [44] H ( 0 ) e ( 0 ) = e ( 0 ) H ( r π ) e ( π ) = 0 G ( 0 ) e ( 0 ) = 0 where H ( ω ) and G ( ω ) are masks of the lowpass and highpass filters (see eq. (3)), and e ( ω ) = [ 1 , e − i ω , e − 2 i ω , . . . , e − ( r − 1 ) i ω ] T . A class of symmetric-antisymmetric biorthogonal multiwave let filters which possess GMPs was introduced in [42]. Matrix Spectral Factorization - SA4 Multiwa velet 7 A GMP order of at least ( 1 , 1 , 1 ) is critical for ensuring no frequency leakage across bands, hence improving compression performance. Multifilters possessing GMPs hav e bet- ter performance than those which do not possess GMPs [31], due to better frequency re- sponses for energy compaction, greater regularity and greater approximation order of the corresponding wav elet/scaling functions. GMPs help pre vent both DC and high-frequency leakage across bands; this contributes to reduced smearing, blocking, ringing artifacts, and also helps to prev ent checkerboard artifacts in reconstructed images for image coding. The SA4 multiwav elet, introduced belo w in section 3.3, is a member of a one-parameter family of orthogonal multiwav elets with four coefficient matrices [44]. Different members of the family hav e different approximation order , but all of them hav e a GMP order of at least ( 1 , 1 , 1 ) . This manifests itself in the smooth decay to zero of the magnitude responses near ω = π , compared to the following 4-tap orthogonal multiw av elet filters: – The GHM multiwav elet [10] has symmetric orthogonal scaling functions and an approx- imation order of 2 for its filter length 4. – Chui and Lian’ s CL multiwav elet [5] has the highest possible approximation order of 3 for its filter length 3. – Jiang’ s multiwav elet JOPT4 [26] has optimal time-frequency localization for its filter length. Although the GHM and CL systems are the most commonly used orthogonal multi- wa velet systems, and hav e higher approximation order than the SA4 system, they do not satisfy GMPs. Smoothness The smoothness of ϕ ϕ ϕ , ψ ψ ψ can be characterized by Sobolev re gularity S ; this is discussed in [33, 6, 25]. The regularity estimate is related to both the eigenv alues and the correspond- ing eigen vectors of the transition operator and to spectral radius of the transition operator . Symmetry and Antisymmetry A function g ( t ) is symmetric about a point c if g ( c − t ) = g ( c + t ) for all t . It is antisym- metric if g ( c − t ) = − g ( c + t ) . For scaling functions and wa velets of compact support, the only possible symmetry is about the center of their support. F or multiscaling and multiwav elet functions, this could be a different point for dif ferent components. In the scalar case, the scaling function cannot be antisymmetric, since it must have a nonzero integral. In the multiwav elet case, some components ev en of ϕ ϕ ϕ can have antisym- metry . 2.3 Problem Statement The design of perfect reconstruction multifilter banks and multiwav elets remains a signif- icant problem. In the scalar case, spectral factorization of a half-band filter that is positi ve definite on the unit circle was, in fact, the first design technique, suggested by Smith and Barnwell [38]. This provides a motiv ation to extend the technique to the multiwav elet case. Spectral factorization of a Hermitian product filter in the multiwavelet case is more chal- lenging than in the scalar case, and has not been in vestigated previously . 8 V asil K olev et al. The present paper considers the application of MSF to the problem of finding an orthog- onal lowpass multifilter H ( z ) , given a product filter P ( z ) . Since ev ery orthogonal multiscal- ing filter H ( z ) is a spectral factor of some matrix product filter , this can be used as a tool for designing new filters. W ith an eye towards future applications, we are interested in constructing multifilters with desirable properties, especially GMPs. This implies that the determinant of P ( z ) will hav e at least one zero of higher multiplicity on the complex unit circle. That is, P ( z ) will be degenerate. As a test case, we want to use a product filter P ( z ) which satisfies the half-band con- dition (2), and is derived from an H ( z ) which has more than 2 taps (i.e. n ≥ 3) and good regularity properties and GMPs. Our benchmark test case will be the 2-channel orthogonal SA4 multiwa velet. While there are many matrix spectral factorization algorithms [7, 24, 27], most cannot handle the degenerate case. W e use Bauer’ s method, based on the T oeplitz method of spec- tral factorization of the Y oula and Kazanjian algorithm [2, 50], which can handle this case. Howe ver , conv ergence becomes very slo w . It is easy to see that matrix spectral factorization is not unique. If H 1 ( z ) is a factor of a giv en P ( z ) , then H 2 ( z ) = H 1 ( z ) U is also a factor for any orthogonal matrix U , or more generally H 2 ( z ) = H 1 ( z ) U ( z ) for paraunitary U ( z ) . Other operations may be permissible in some cases. In our experiments with Bauer’ s algorithm, the initial factors do not correspond to smooth functions. W e apply regularization techniques, both to speed up conv ergence and to work to wards recov ering the original filter in this test case. In this manner , we deri ve tw o filters which we call the appr oximate and exact SA4 multiwavelets . Let us briefly summarize some previous w ork. Matrix Spectral Factorization (MSF) plays a crucial role in the solution of various ap- plied problems for MIMO systems in communications and control engineering [46]. It has been applied to designing minimum phase FIR filters and the associated all-pass filter [19], quadrature-mirror filter banks [3], MIMO systems for optimum transmission and reception filter matrices for precoding and equalization [15], precoders [12], and many other applica- tions. The analysis of linear systems corresponding to a given spectral density function was first established by W iener, who used linear prediction theory of multi-dimensional stochas- tic processes [47]. The method of Wilson for MSF was developed for applications in the prediction theory of multiv ariate discrete time series, with known (or estimated) spectral density [48]. Many numerical approaches to MSF have been proposed, for example [7, 24, 27]. A surve y of such algorithms is giv en in [36]; ho wev er, 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 ˇ cera’ 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 [35] 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 ac- curacy around 10 − 17 after 10,000 recursions on a 64-bit floating point processor . Matrix Spectral Factorization - SA4 Multiwa velet 9 3 Matrix Spectral Factorization The following theorem answers the e xistence question for MSF . Theorem 1 (Matrix form of the Riesz-F ej ´ er lemma [20, 13]) A matrix polynomial P ( z ) = n ∑ k = − n P k z k can be factor ed as P ( z ) = H ( z ) H ∗ ( z ) , wher e H ( z ) = n ∑ k = 0 C k z − k is an r-channel causal polynomial matrix, if and only if P ( z ) is symmetric positive semidefi- nite for z on the complex unit cir cle. 3.1 Bauer’ s Method If H ( z ) and P ( z ) are kno wn, we can construct doubly infinite block T oeplitz matrices L and T from their coefficients, by setting L i j = C i − j , T i j = T j − i . T is symmetric and block banded with bandwidth n ; L is block lower triangular , also with bandwidth n . The relation P ( z ) = H ( z ) H ∗ ( z ) corresponds to T = L L T , that is, L is a Cholesky factor of T . In Bauer’ s method, we pick a large enough integer f , and truncate T to (block) size f × f : T ( f ) =                P 0 P 1 · · · P n P − 1 P 0 P 1 · · · P n . . . . . . . . . P − n . . . P n . . . . . . . . . . . . . . . P 1 P − n · · · P − 1 P 0                If T ( f ) is positiv e definite, we can compute the Cholesky factorization T = L ( f ) L ( f ) T , where L ( f ) =             C ( 1 ) 0 C ( 2 ) 1 C ( 2 ) 0 . . . . . . C ( n + 1 ) n C ( n + 1 ) 0 . . . . . . C ( f ) n · · · C ( f ) 1 C ( f ) 0             It is shown in [50] that this f actorization is possible, and that C ( f ) k → C k as f → ∞ . Bauer’ s method works e ven in the case of highly degenerate P ( z ) . Unfortunately , con- ver gence in this case is very slo w . 10 V asil K olev et al. 0 10 20 30 40 10 -6 10 -4 10 -2 10 0 f σ (T) σ (L SF ) Fig. 2 A sharp drop of the singular v alues of the T oeplitz matrix T and of its spectral factor L SF , for f = 20; (filled circles) singular values of T ; (hollow circles) singular v alues of L SF . 3.2 Numerical Behavior For chosen size f , the matrix T ( f ) is of size 2 f × 2 f . Its singular values are all in the range [ 0 , 2 ] . The first ( f − 1 ) singular v alues are close to 2, then σ f = σ f + 1 = 1, and the remaining ( f − 1 ) singular v alues 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 slo w con vergence. For example, in the scalar case p ( z ) = z − 1 + 2 + z , where the factor is h ( z ) = 1 + z − 1 , the value on the diagonal in ro w f is s 1 + 1 f ≈ 1 + 1 2 f , which con verges 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, con vergence speed becomes e ven worse. 3.3 T esting Bauer’ s Method A flowchart for testing Bauer’ s method for MSF is shown in fig. 3. The steps are Step 1: Choice of the benchmark multiscaling function H ( z ) ; Step 2: Construction of the matrix lowpass product filter P ( z ) ; Step 3: Find a spectral factor L SF = L ( f ) , and assess its quality; Step 4: Do different kinds of postprocessing to obtain spectral factors H ( i ) , i = 1 , 2 , 3; Step 5: Compare the obtained factors with the benchmark H ( z ) . Step 1: W e consider the well-known orthogonal SA4 multiwav elet [5, 23, 44] as a bench- mark testing case. The recursion coefficients C k and D k of SA4 depend on a parameter t . They are C 0 = α  1 t 1 − t  ; C 1 = α  t 2 t − t 2 t  ; C 2 = α  t 2 − t t 2 t  ; C 3 = α  1 − t − 1 − t  ; D 0 = α  t − 1 t 1  ; D 1 = α  − t t 2 t t 2  ; D 2 = α  − t − t 2 − t t 2  ; D 0 = α  t 1 − t 1  ; (4) Matrix Spectral Factorization - SA4 Multiwa velet 11 Fig. 3 Steps in testing matrix spectral factorization using Bauer’ s method. T able 1 Coefficients of the original SA4 multiscaling and multiw avelet functions. n C n D n 0 0 . 011226792152545 0 . 088388347648318 − 0 . 088388347648318 0 . 011226792152545 0 . 011226792152545 − 0 . 088388347648318 − 0 . 088388347648318 − 0 . 011226792152545 1 0 . 695879989034003 0 . 088388347648318 0 . 088388347648318 − 0 . 695879989034003 − 0 . 695879989034003 0 . 088388347648318 − 0 . 088388347648318 − 0 . 695879989034003 2 0 . 695879989034003 − 0 . 088388347648318 0 . 088388347648318 0 . 695879989034003 0 . 695879989034003 0 . 088388347648318 0 . 088388347648318 − 0 . 695879989034003 3 0 . 011226792152545 − 0 . 088388347648318 − 0 . 088388347648318 − 0 . 011226792152545 − 0 . 011226792152545 − 0 . 088388347648318 0 . 088388347648318 − 0 . 011226792152545 where α = 1 / ( √ 2 ( 1 + t 2 )) . In this paper , we use the filter with t = 4 + √ 15, which leads to α = ( 4 − √ 15 ) / ( 8 √ 2 ) . Numerical values of C k , D k are listed in table 1. These coefficients generate the 2-band compactly supported orthogonal multiscaling function ϕ ϕ ϕ = [ φ 0 , φ 1 ] T and multiwavelet function ψ ψ ψ = [ ψ 0 , ψ 1 ] T , all with support [ 0 , 3 ] . Graphs of these functions can be found in fig. 8 in a later chapter . Step 2: The symbol of the product lowpass multifilter P has the form P ( z ) = P T 3 z − 3 + P T 1 z − 1 + P 0 + P 1 z + P 3 z 3 , where P 0 = I , P 1 = P T − 1 = 1 64  4 √ 15 + 17 4 √ 15 + 16 − 4 √ 15 − 16 − 4 √ 15 − 17  , P 2 = 0 , P 3 = P T − 3 = 1 64  15 − 4 √ 15 4 √ 15 − 16 16 − 4 √ 15 4 √ 15 − 15  . Numerical values of P k are giv en in table 2. The impulse responses of the product filter (see fig. 4) hav e the character of nearly Haar- type scaling and wa velet functions. The frequency responses ha ve good selectivity . 12 V asil K olev et al. T able 2 The matrix coefficients P k of the half-band product filter P k P k 0 1 0 0 1 1 0 . 507686459137964 0 . 492061459137964 − 0 . 492061459137964 − 0 . 507686459137964 3 − 0 . 007686459137964 − 0 . 007938540862036 0 . 007938540862036 0 . 007686459137964 4 5 6 -5 0 5 0 0.2 0.4 0.6 0.8 1 -150 -100 -50 Normalized Frequency ( ×π rad/sample) Magnitude (dB) P 0 P 1 Fig. 4 The impulse response (left) and frequency response (right) of the two components of the product filter P = [ P 0 ; P 1 ] T . Step 3: The matrix T ( f ) in this case looks like this: T ( f ) =                P 0 P T 1 0 P T 3 P 1 P 0 P T 1 0 P T 3 0 . . . . . . P 3 . . . P T 3 . . . . . . 0 . . . . . . P T 1 P 3 0 P 1 P 0                For any fixed f , the coefficients from the last row of L ( f ) are approximate factors of P ( z ) . Suppressing the dependence on f in the notation, we define the filter L SF ( z ) = n ∑ k = 0 C ( f ) k z − k and the matrix L SF = h C ( f ) 0 C ( f ) 1 C ( f ) 2 C ( f ) 3 i (5) =  L 00 L 01 L 02 L 03 L 04 L 05 L 06 L 07 L 10 L 11 L 12 L 13 L 14 L 15 L 16 L 17  . (6) Matrix Spectral Factorization - SA4 Multiwa velet 13 T able 3 Coefficients of the lo wpass spectral factor L ( f ) ( z ) for f = 81, with ∆ P ( 81 ) = 3 . 169 · 10 − 9 . k C ( 81 ) k 0 0.094428373297668 0 -0.091754813647953 0.022310801334572 1 0.175943193428843 0.679731045265413 0.010360133828403 -0.702056243347588 2 0.000129414745433 0.700756678587017 0.165030422671601 0.681046913905671 3 -0.081190837390599 0.021002135513698 -0.083853536287580 -0.001275235049147 This notation will be used again later . T o measure the quality of the factorization, we compute ∆ P ( z ) = P ( z ) − L SF ( z ) L SF ∗ ( z ) = n ∑ k = − n ∆ P k z k . The residual is then ∆ P = max k = − n ,..., n max i j | ( ∆ P k ) i j | . In numerical experiments, the minimal error ∆ P is achie ved for f = 81; the coef ficients C ( 81 ) k are giv en in table 3. The numerical error ∆ P ( 81 ) = 3 . 169 · 10 − 9 for Bauer’ s method is much better than that of W ilson’ s method for MSF , which cannot be lower then 10 − 5 [14]. The impulse response of the obtained spectral factor is non-regular , but the frequency response has good selectivity . This approximate f actor does not quite provide perfect recon- struction. This is addressed in more detail in subsection 5.2. Steps 4 and 5 will be cov ered in chapter 4. 4 Construction of Appr oximate and Exact SA4 Multiwavelets Abov e, we showed how to obtain the approximate spectral factor L SF for a giv en size f . Giv en an y f actor of P , other possible f actors 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 rev ersal to produce a lowpass filter that is similar to, but not identical to, the original SA4 multiscaling function. W e call this the appr oximate SA4 multiwavelet . In subsection 4.2, we add a rotation postfactor and an averaging process to produce another multiw av elet which is even closer to the original SA4 multiwa velet. W e call this the exact SA4 multiwavelet . The newly deri ved multiscaling filters are denoted by H ( 1 ) = h C ( 1 ) 0 , C ( 1 ) 1 , C ( 1 ) 2 , C ( 1 ) 3 i for the approximate SA4 multiwa velet, and H ( 3 ) = h C ( 3 ) 0 , C ( 3 ) 1 , C ( 3 ) 2 , C ( 3 ) 3 i 14 V asil K olev et al. for the exact SA4 multiwa velet. ( H ( 2 ) is an intermediate step in the calculation of H ( 3 ) ). In each case, we find the corresponding multiwavelet filter by using a QR decomposition [18]. The multiscaling filter is factorized as     C 0 C 1 C 2 C 3     = QR =   | | q 1 · · · q 8 | |   ·     I 0 0 0     The coefficients D k are then obtained from the third and fourth columns of Q . W e define the absolute error by ∆ H ( k ) = H − H ( k ) = h ∆ C ( k ) 0 , ∆ C ( k ) 1 , ∆ C ( k ) 2 , ∆ C ( k ) 3 i , where ∆ C ( k ) n = C n − C ( k ) n for k = 1 , 2 , 3, and like wise for ∆ G ( k ) . T o measure the deviation of the new filters from the original, we introduce the mean squar e error (MSE) of the multiscaling and multiw av elet functions, MSE-MF ( k ) = ∑ i j | ∆ H ( k ) i j | 2 , MSE-MwF ( k ) = ∑ i j | ∆ G ( k ) i j | 2 , and the maximal absolute err ors (MAE) of the matrix coefficients, the multiscaling function and the multiwa velet function as MAE-MC ( k ) ` = max i j | [ ∆ C ( k ) ` ] i j | , MAE-MF = max i j | ∆ H ( k ) i j | = max i j | H i j − H ( k ) i j | , MAE-MwF = max i j | ∆ G ( k ) i j | = max i j | G i j − G ( k ) i j | . 4.1 Approximate SA4 Multiwa velet The spectral factor L SF 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 f actor: L SF 0 , 6 = 1 4  | L SF 00 | + | L SF 10 | + | L SF 06 | + | L SF 16 |  , L SF 1 , 7 = 1 4  | L SF 01 | + | L SF 11 | + | L SF 07 | + | L SF 17 |  , L SF 2 , 4 = 1 4  | L SF 02 | + | L SF 12 | + | L SF 04 | + | L SF 14 |  , L SF 3 , 5 = 1 4  | L SF 05 | + | L SF 15 | + | L SF 05 | + | L SF 15 |  . and construct av erage matrices  L SF 2 , 4 L SF 3 , 5 L SF 2 , 4 L SF 3 , 5  ,  L SF 0 , 6 L SF 1 , 7 L SF 0 , 6 L SF 1 , 7  . Matrix Spectral Factorization - SA4 Multiwa velet 15 20 40 60 80 100 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 f MAE-MF (1) MSE-MC 0 (1) MSE-MC 1 (1) Fig. 5 Dependence of maximum absolute error MAE-MF ( 1 ) and mean squared errors MSE-MC ( 1 ) 0 = MSE- MC ( 1 ) 3 , MSE-MC ( 1 ) 1 = MSE-MC ( 1 ) 2 on matrix size f (log. scale). Second, the matrix coefficients C ( 1 ) 1 and C ( 1 ) 3 are obtained by multiplying the a veraged ma- trices from the right by J = diag ( 1 , − 1 ) , but keep the signs, while the matrix coefficients C ( 1 ) 0 and C ( 1 ) 2 are obtained by multiplying C ( 1 ) 1 and C ( 1 ) 3 on the left and right by U = antidiag ( 1 , 1 ) : C ( 1 ) 3 = sign ( L SF 3 )  L SF 0 , 6 L SF 1 , 7 L SF 0 , 6 L SF 1 , 7  · J , C ( 1 ) 0 = U · C ( 1 ) 3 · U , C ( 1 ) 1 = sign ( L SF 1 )  L SF 2 , 4 L SF 3 , 5 L SF 2 , 4 L SF 3 , 5  · J , C ( 1 ) 2 = U · C ( 1 ) 1 · U . This produces coefficients quite close to the original coef ficients in (4). In order to determine the best approximate solution, we inv estigate the influence of the size f of the T oeplitz matrix T ( f ) for the numerical errors MAE-MF ( 1 ) , MSE-MC ( 1 ) 0 , and MSE-MC ( 1 ) 1 . These errors are shown in fig. 5 in logarithmic scale, for f up to 100. Note that MSE-MC ( 1 ) 3 = MSE-MC ( 1 ) 0 , and MSE-MC ( 1 ) 2 = MSE-MC ( 1 ) 1 , by construction. The important minimal errors MAE-MF ( 1 ) and MSE-MwF ( 1 ) for f = 21 are tabulated in table 4; the matrix coefficients of are listed in table 5. Obviously , the de velopment of the errors shows the con vergence of the algorithm, as well as the dependency of the errors on the size f . It also shows that the main error comes from the matrix coef ficients C ( 1 ) 0 and C ( 1 ) 3 , rather than C ( 1 ) 1 and C ( 1 ) 2 . The error MSE-MF ( 1 ) behav es essentially the same way as MAE-MC ( 1 ) 0 . Unfortunately , the con vergence of MAE-MC ( 1 ) 0 and MAE-MC ( 1 ) 3 stalls at a relatively large number . Thus, the approximate SA4 multiwav elet is not very accurate. Better matrix coefficients are achiev ed in the next subchapter . 16 V asil K olev et al. T able 4 The minimum MAEs and MSEs for the approximate and exact SA4 multiscaling and multiwavelet functions, the size f and the corresponding angle θ . f MAE-MF ( 1 ) MAE-MwF ( 1 ) MSE-MF ( 1 ) MSE-MwF ( 1 ) 21 4 . 2797 · 10 − 3 4 . 6139 · 10 − 6 f θ MAE-MF ( 2 ) MAE-MwF ( 2 ) MSE-MF ( 2 ) MSE-MwF ( 2 ) 15,400 1 . 444628609880879 1 . 220 · 10 − 4 7 . 376 · 10 − 9 f θ MAE-MF ( 3 ) MAE-MwF ( 3 ) MSE-MF ( 3 ) MSE-MwF ( 3 ) 12,042 1 . 444630025395427 8 . 285 · 10 − 8 8 . 557 · 10 − 8 3 . 835 · 10 − 15 3 . 717 · 10 − 15 T able 5 Coefficients of the approximate SA4 multiscaling and multiw avelet functions for f = 21. n C ( 3 ) n D ( 3 ) n 0 0 . 011165766264837 0 . 088552225597447 − 0 . 088552225597447 0 . 011165766264837 0 . 011165766264837 − 0 . 088552225597447 − 0 . 088552225597447 − 0 . 011165766264837 1 0 . 691600252880066 0 . 088718768987217 0 . 088718768987217 − 0 . 691600252880066 − 0 . 691600252880066 0 . 088718768987217 − 0 . 088718768987217 − 0 . 691600252880066 2 0 . 691600252880066 − 0 . 088718768987217 0 . 088718768987217 0 . 691600252880066 0 . 691600252880066 0 . 088718768987217 0 . 088718768987217 − 0 . 691600252880066 3 0 . 011165766264837 − 0 . 088552225597447 − 0 . 088552225597447 − 0 . 011165766264837 − 0 . 011165766264837 − 0 . 088552225597447 0 . 088552225597447 − 0 . 011165766264837 4.2 Exact SA4 Multiwa velet The second algorithm leads to more regular scaling functions. W e multiply the coefficients of L SF from the right by the unitary matrix U ( θ ) C ( 2 ) k = C SF k · U ( θ ) = C SF k ·  cos θ sin θ − sin θ cos θ  , k = 0 , . . . , 3 . This produces the approximate factor H ( 2 ) . The angle θ (in radians) is calculated as the average v alue θ = 1 2  cos − 1 ( ev en ) + sin − 1 ( odd )  , where ev en = 1 4 ( | L 00 + L 04 + L 10 + L 14 | + | L 02 + L 06 + L 12 + L 16 | ) and odd = 1 4 ( | L 01 + L 05 + L 11 + L 15 | + | L 03 + L 07 + L 13 + L 17 | ) . Again, we are using the notation from eq. (5). Figure 6 shows the dependence of the angle θ on the size f (log. scale). It can be seen that in the angle there is one ov ershoot, after that it is decreasing very slo wly . The minimum errors are obtained at the maximal possible size of T oeplitz matrix f = 15 , 400, with the values shown in table 4. This shows the influence of the very slow conv ergence; despite right multiplication with the unitary matrix, desirable precison cannot be an achieved. The errors MSE-MF ( 2 ) and MAE-MF ( 2 ) for H ( 2 ) are shown in fig. 7. Matrix Spectral Factorization - SA4 Multiwa velet 17 10 0 10 1 10 2 10 3 10 4 1.36 1.38 1.4 1.42 1.44 1.46 f θ Fig. 6 Dependence of the angle θ (radians) on size f (log scale). T able 6 Coefficients of the e xact SA4 multiscaling and multiwavelet functions for f = 12 , 042. n C ( 3 ) n D ( 3 ) n 0 0 . 01122679201336641 0 . 08838843031594616 − 0 . 08838843285624247 0 . 01122679262734716 0 . 01122679228175923 − 0 . 08838843013543155 − 0 . 08838843302944295 − 0 . 01122679235802515 1 0 . 6958799459541121 0 . 08838843050191729 0 . 08838843089794699 − 0 . 6958799676294173 − 0 . 6958799462085388 0 . 08838842817713602 − 0 . 08838843321541406 − 0 . 6958799673174058 2 0 . 6958799459541121 − 0 . 08838843050191729 0 . 08838843089794699 0 . 6958799676294172 0 . 6958799462085388 0 . 08838842817713602 0 . 08838843321541406 − 0 . 6958799673174058 3 0 . 01122679201336641 − 0 . 08838843031594616 − 0 . 08838843285624255 − 0 . 01122679262734722 − 0 . 01122679228175923 − 0 . 08838843013543155 0 . 08838843302944295 − 0 . 01122679235802536 T o increase the numerical precision and obtain the exact SA4 multiscaling function, we again apply an av eraging approach. H ( 3 ) = [ C ( 3 ) 0 , C ( 3 ) 1 , C ( 3 ) 2 , C ( 3 ) 3 ] , ( C ( 3 ) 0 ) i j = sign (( C ( 2 ) 0 ) i j ) · 1 2     ( C ( 2 ) 0 ) i j    +    ( C ( 2 ) 3 ) i j     , ( C ( 3 ) 1 ) i j = sign (( C ( 2 ) 1 ) i j ) · 1 2     ( C ( 2 ) 1 ) i j    +    ( C ( 2 ) 2 ) i j     , ( C ( 3 ) 2 ) i j = sign (( C ( 2 ) 2 ) i j ) · 1 2     ( C ( 2 ) 1 ) i j    +    ( C ( 2 ) 2 ) i j     , ( C ( 3 ) 3 ) i j = sign (( C ( 2 ) 3 ) i j ) · 1 2     ( C ( 2 ) 0 ) i j    +    ( C ( 2 ) 3 ) i j     . The errors for H ( 3 ) 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 ( 3 ) . There are 3 minima, as can be seen in the zoomed part of fig. 7. The first local minimum is at the value f = 4 , 269, the second at f = 9 , 139, the third at f = 12 , 042. This means that increasing the size of the T oeplitz matrix leads to slow improv ement in the precision of the exact multiscaling function. Nev ertheless, due to applying the averaging method, the minimal errors in H ( 3 ) and G ( 3 ) for f = 12 , 042 are smaller by 2 or 3 orders of magnitude than the errors in H ( 2 ) and G ( 2 ) , and smaller by 5 orders of magnitude than the errors in H ( 1 ) and G ( 1 ) . W e choose H ( 3 ) for f = 12 , 042 for the coefficients of the exact SA4 multiwa velet, giv en in table 6. They are almost the same as the original coef ficients in (4). 18 V asil K olev et al. 0 5000 10000 15000 10 -15 10 -10 10 -5 f 4200 4250 4300 10 -7 X: 4269 Y: 1.004e-07 9120 9140 9160 10 -7 X: 9139 Y: 8.624e-08 1.202 1.204 1.206 1.208 x 10 4 10 -7 X: 1.204e+04 Y: 8.285e-08 MAE-MF (2) MSE-MF (2) MAE-MF (3) MSE-MF (3) Fig. 7 Dependence of the errors MSE-MF ( i ) and MAE-MF ( i ) i = 2 , 3, on the size of f (log scale). Matrix Spectral Factorization - SA4 Multiwa velet 19 T able 7 Comparison of coding gain (CG), Sobole v regularity (S), symmetry/antisymmetry (S/A), and length, for various multiw avelets. Multifilter CG S S/A Length Biorthogonal Hermitian cubic spline (bih34n, [39]) 1.51 2.5 yes 3/5 Biorthogonal (dual to Hermitian cubic spline) (bih32s, [1]) 1.01 2.5 yes 3/5 Integer Haar [4] 1.83 0.5 yes 2 Chui-Lian (CL, [5]) 2.06 1.06 yes 3 Biorthogonal (dual to Hermitian cubic spline) (bih54n, [1]) 2.42 0.61 yes 5/3 Biorthogonal (from GHM) (bighm2, [1]) 2.43 0.5 yes 2/6 Biorthogonal (from GHM, dual to bighm2) (bighm6, [1]) 3.53 2.5 yes 6/2 Biorthogonal (dual to Hermitian cubic spline) (bih52s, [45]) 3.69 0.83 yes 5/3 Approximate SA4 3.73 0.99 yes 4 Exact SA4 ([37]) 3.73 0.99 yes 4 Geronimo-Hardin-Massopust (GHM, [17]) 4.41 1.5 yes 4 5 Perf ormance Analysis 5.1 Comparison W ith Other Multiwavelets T able 7 giv es the results of our experiments with the approximate and exact SA4 multi- wa velets, compared with other well-known orthogonal and biorthogonal filter banks. Com- parisons include coding gain (CG), Sobolev smoothness (S) [8], symmetry/antisymmetry , and length. All of the multiwa velets considered ha ve symmetry/antisymmetry . The CG for orthogonal transforms is a good indication of the performance in signal processing. It is the ratio of arithmetic and geometric means of channel variances σ 2 i : CG = 1 r ∑ r i = 1 σ 2 i  ∏ r i = 1 σ 2 i  1 / r Coding gain is one of the most important factors to be considered in many applications. It is always greater than or equal to 1; greater values are better . CG is equal to 1 if all the v ariances are equal, which means that it is not possible to clearly distinguish between the smooth and the detailed components of the multiwav elet transformation coefficients. T o estimate CG, the variance is computed using a first order Markov model AR(1) with intersample autocorrelation coefficient ρ = 0 . 95 [11]. The Sobolev exponent S of a filter bank measures the L 2 -differentiability of the cor- responding multiscaling function ϕ ϕ ϕ = [ φ 0 , φ 1 ] T , and thus also the multiwav elet function ψ ψ ψ = [ ψ 0 , ψ 1 ] T . It is completely determined by the multiscaling symbol H ( z ) . The obtained Sobolev regularity and CGs of the approximate and exact SA4 multi- wa velets are equal: S ( 1 ) = S ( 3 ) = 0 . 9919, and C G ( 1 ) = C G ( 3 ) = 3 . 7323 dB. This means that for some applications we can use the SA4 multiwa velet. According to table 7, the approximate and exact SA4 multiwavelets are better than most commonly used filter banks. They also allow an economical lifting scheme for future imple- mentation. Fig. 8 shows a comparison of the approximate and exact SA4 multiwav elets in time domain and impulse response. The influence of the lo wer precision of the approximate SA4 multiwa velet can be seen in the time domain image. The magnitude of the approximate SA4 in the time domain is observably smaller than for the exact SA4, while the frequency responses are essentially identical (see fig. 8). There- 20 V asil K olev et al. 0 1 2 3 0 0.5 1 φ 0 φ 0 A 0 1 2 3 -1 0 1 φ 1 φ 1 A 0 1 2 3 -1 -0.5 0 0.5 1 ψ 0 ψ 0 A 0 1 2 3 -1 -0.5 0 0.5 1 ψ 1 ψ 1 A 0 0.5 1 -80 -60 -40 -20 Normalized Frequency ( ×π rad/sample) Magnitude (dB) φ A 0 φ 0 φ A 1 φ 1 0 0.2 0.4 0.6 0.8 1 -80 -60 -40 -20 Normalized Frequency ( ×π rad/sample) Magnitude (dB) ψ A 0 ψ 0 ψ A 1 ψ 1 Fig. 8 T op: The impulse responses of the components of the exact SA4 multiwa velet ϕ ϕ ϕ = ( φ 0 , φ 1 ) T , ψ ψ ψ = ( ψ 0 , ψ 1 ) T (in blue) , compared with the approximate multiwa velet ϕ ϕ ϕ A = ( φ A 0 , φ A 1 ) T , ψ ψ ψ A = ( ψ A 0 , ψ A 1 ) T (in red) ; left to right: φ 0 , φ 1 , ψ 0 , ψ 1 ; Bottom: The frequency response of the components of the exact SA4 multiwav elet (in blue) compared with the approximate multiwav elet (in red) ; left: φ 0 , φ 1 ; right: ψ 0 , ψ 1 . fore, in applications where the frequency response is important, we can use the approximate multiscaling function. 5.2 Influence of Inaccuracy of Filter Coef ficients with Approximate SA4 Multiwavelet The approximate SA4 multiwa velet H ( 1 ) is not quite an orthogonal perfect reconstruction filter . In this section, we explore whether it can still be useful in applications. 5.2.1 1D Signal (No Additional Pr ocessing) The influence of the inaccuracy of the matrix filter coefficients in the case of the approxi- mate SA4 multiwavelet is shown by 1D applications with no additional processing. By “no additional processing” we mean that we simply decompose and reconstruct a signal through sev eral levels. No compression, denoising, or other processing is done. W e use Haar balancing pre- and postfilters [8, 40] Q = 1 √ 2  1 1 − 1 1  , whose small length preserves the time localization of the multiw avelet decomposition, sim- plicity , orthogonality , and symmetry . Results for the balanced version are sho wn in fig. 9(b). Matrix Spectral Factorization - SA4 Multiwa velet 21 The error in orthogonality of the approximate balanced SA4 multiwa velet is ∆ H = I − Q · H ( 1 ) · ( H ( 1 ) ) T Q T = 0 . 0117 I . For the quality measure, we decompose and reconstruct fi ve normed test signals of length N = 2 7 , s (without noise) and ˆ s = s + ε (with noise). The number of decomposition lev els is 6. The noise components ε i are independent identically distributed random v ariables with mean 0 and standard deviation σ . W e use maximum absolute errors (MAE) ∆ s = max i | reconstructed s i − original s i | , where s i refers to the i th test signal. The test signals are ’Cusp’ , ’HiSine’ , ’LoSine’ , ’Piece-r e gular’ , and ’Piece-polynomial’ , implemented in the Matlab en vironment. See fig. 9(a). For the test signal ’Cusp’ , increasing the number of decomposition and reconstruction lev els j leads to a linear increase of MAE. From the obtained small differences (about 3 . 5 · 10 − 4 ) in error between noisy and noiseless signal, it follows that the presence of noise does not make much dif ference in signals of this type. For the test signal ’HiSine’ we observe nearly constant MAEs for the noisy and noiseless signals with increasing lev el j ; the difference between the noisy and noiseless signal is in the interval [ 2 − 5 . 3 ] · 10 − 3 . Again, the presence of noise has only a minimal influence. For the test signal ’LoSine’ , both noisy and noiseless MAEs sho w a linear increase up to the second lev el, with only a small increase at the third and fourth le vel of − 5 . 7 · 10 − 3 and 3 . 6 · 10 − 3 , as shown in fig. 9(b). Therefore, after the second lev el the influence of noise is weak. For the test signal ’Piece-r e gular’ , the MAE grows quadratically with increasing j . The difference between noisy and noiseless signals is smallest at a pre-/post-processing step (6 · 10 − 5 and largest at the second le vel (1 . 2 · 10 − 3 ). For the test signal ’Piece-polynomial’ , the MAE grows non-uniformly and non-linearly with increasing j . The dif ference between noisy and noiseless signals is smallest at the pre-/post-processing step (3 . 8 · 10 − 4 ) and largest at the third lev el (2 · 10 − 3 ). 5.2.2 2D Signal (No Additional Pr ocessing) The performance of the approximate SA4 multiwav elet was tested by decomposing and reconstructing se veral gray le vel images, through 6 or 7 le vels. Fig. 10 shows the details for four of the images, with the quality measure PSN R = log 10 ( 255 2 / M SE ) dB. The images used are ’Lena’ , ’P eppers’ , ’Girlface’ , and ’Barbara’ . The approximate balanced SA4 multiwav elet applied to the four images leads to an exponential decrease of the PSNRs. According to these results, it is preferable to use no more than 3 or 4 le vels of decomposition. After that, the reconstructed image has very lo w PSNRs and visibly worse quality . In some applications, such as big data archiv es, higher lev els of decomposition may be useful. 5.3 Image Denoising In this section, we compare the balanced and non-balanced version of the exact SA4 multi- wa velet with the GHM multiwavelet (Do wnie and Silverman 1998) and the Chui-Lian mul- tiwa velet CL (Downie and Silverman 1998), by considering image denoising with vector 22 V asil K olev et al. (a) 0 50 100 -1 0 1 Piece-regular 0 50 100 -1 0 1 Piece-polynomial 0 50 100 0 0.5 1 Cusp 0 50 100 -1 0 1 LoSine 0 10 20 30 40 50 60 70 -1 0 1 HiSine (b) ·�                                                                      &    &   Fig. 9 Decomposition and reconstruction with the balanced approximate SA4 multiwa velet; (a) The normed test signals; (b) The MAEs obtained at lev el j . hard thresholding (Donoho and Johnstone 1994), using 1-5 lev els. The images are ’Lena’ , ’Zelda’ , and ’House’ , of size 512 × 512 pixels, with white additiv e Gaussian noise with variance σ = 10. See fig. 10. The multiw av elet coef ficients of the white noise are reduced at each le vel, b ut uniformly distributed within each lev el. Therefore, the best approach to denoising is to find an appro- priate threshold value at each le vel. The exact SA4 multiwavelet, both balanced and non-balanced, achiev ed the maximal PSNRs. The PSNR differences for the both versions of the test images are 3 . 57 dB at the Matrix Spectral Factorization - SA4 Multiwa velet 23 (a) (b) (c) (d) (e) (f) Fig. 10 T est images of size 256 × 256 pixels) (a) Lena, (b) Peppers, and (512 × 512 pixels) (c) House, (d) Girlface, (e) Barbara, and (f) Zelda. Fig. 11 The PSNRs of four images for 2D decomposition and reconstruction without additional processing of 1 through 6 or 7 lev els with the approximate SA4 multiwav elet. pre-/post-processing step and 4 . 06 dB at the fifth lev el for “Lena” ; 4 . 17 dB (pre-/post- processing step) and 5 . 13 dB (fifth level) for “Zelda” , 5 . 09 dB (pre-/post-processing step) and 6 . 53 dB (fifth level) for House” . Although balancing of multiwa velets destroys the symmetry , it leads to increasing PSNRs and better image denoising with the exact SA4 multiwa velet for the three test images for the both v ersion multiwa velets, while for the non- balanced GHM and CL multiwa velets PSNRs decrease (see fig. 12(a)). According to PSNRs, image decomposition and reconstruction through two lev els with the approximate SA4 multiwa velet is comparable to image denoising with the non-balanced SA4 multiwavelet through three lev els, while one lev el is comparable with the balanced multiwa velet pre-/post-processing step (see fig. 11 and fig. 12). 24 V asil K olev et al. Fig. 12 Comparative analysis of PSNRs for image denoising with the exact SA4 multiwavelet and vector hard threshold through 2-6 lev els of the images ’House’ , ’Zelda’ , and ’Lena’ with 512 × 512 pixels and white additi ve Gaussian noise with variance σ = 10. Left: Non-balanced multiwav elets; Right: Balanced multiwav elets. Obviously , the PSNRs of the exact balanced SA4 multiwav elet are better than the well- known orthogonal multiw avelets GHM and CL. They are useful in man y applications. 6 Conclusion In this paper , we consider for the first time the problem of obtaining an orthogonal multi- scaling function by matrix spectral factorization from a degenerate polynomial matrix. W e show benchmark testing of MSF , and apply Bauer’ s method to factoring the product filter of the SA4 multiwa velet. In addition, we sho w how to remov e numerical errors and improve the properties of the factors obtained from Cholesky factorization, leading to fast con vergent algorithms. A very important part is obtaining the ke y angle θ in explicit form. Based on the pro- posed av eraging approach, we de velop two filter banks, the approximate and the exact SA4 orthogonal multiwa velets. Experimental results hav e shown that the performances of the resulting multiwav elets are better than those of the Chui-Lian multiwavelet and biorthogonal multifilters, and are highly comparable to that of longer multiwa velets. Theoretical analyses for the influence of the size f of the T oeplitz matrix are considered, as well as simple 1D and 2D applications. After comparing both types of multifilters, we concluded that the proposed av eraging approach is a better way to remove numerical errors and find the exact SA4 multiwav elet filter bank. It is important to note that the performance of the balanced exact SA4 multi- wa velet for image denoising is better than the well-known orthogonal multiwav elets GHM and CL, which are longer . Acknowledgements The authors would like to thank the three anonymous referees for their critical revie w and helpful suggestions that allowed impro ving the exposition of the manuscript. Matrix Spectral Factorization - SA4 Multiwa velet 25 References 1. A verbuch, A.Z., Zheludev , V .A., Cohen, T .: Multiwavelet frames in signal space originated from Hermite splines. IEEE Trans. Signal Process. 55 (3), 797–808 (2007) 2. Bauer, F .L.: Beitr ¨ age zur Entwicklung numerischer Verfahren f ¨ ur programmgesteuerte Rechenanlagen. II. Direkte Faktorisierung eines Polynoms. Bayer . Akad. Wiss. Math.-Nat. Kl. S.-B. 1956 , 163–203 (1957) (1956) 3. Charoenlarpnopparut, C.: One-dimensional and multidimensional spectral factorization using gr ¨ obner basis approach. In: 2007 Asia-Pacific Conference on Communications, pp. 201–204 (2007) 4. Cheung, K.W ., Po, L.M.: Integer multiwa velet transform for lossless image coding. In: Intelligent Mul- timedia, V ideo and Speech Processing, 2001. Proceedings of 2001 International Symposium on, pp. 117–120 (2001) 5. Chui, C.K., Lian, J.A.: A study of orthonormal multi-wavelets. Appl. Numer . Math. 20 (3), 273–298 (1996). Selected keynote papers presented at 14th IMA CS W orld Congress (Atlanta, GA, 1994) 6. Cohen, A., Daubechies, I., Plonka, G.: Regularity of refinable function vectors. J. Fourier Anal. Appl. 3 (3), 295–324 (1997) 7. Cooklev , T ., Nishihara, A., Kato, M., Sablatash, M.: T wo-channel multifilter banks and multiwavelets. In: Conference Proceedings, IEEE International Conference on Acoustics, Speech, and Signal Processing ICASSP-96, vol. 5, pp. 2769–2772 (1996) 8. Cotronei, M., Montefusco, L.B., Puccio, L.: Multiwa velet analysis and signal processing. IEEE T rans- actions on Circuits and Systems II: Analog and Digital Signal Processing 45 (8), 970–987 (1998) 9. Donoho, D.L., Johnstone, I.M.: Ideal spatial adaptation by wav elet shrinkage. Biometrika 81 (3), 425– 455 (1994) 10. Donovan, G.C., Geronimo, J.S., Hardin, D.P ., Massopust, P .R.: Construction of orthogonal wavelets using fractal interpolation functions. SIAM J. Math. Anal. 27 (4), 1158–1192 (1996) 11. Downie, T .R., Silverman, B.W .: The discrete multiple wa velet transform and thresholding methods. IEEE T ransactions on Signal Processing 46 (9), 2558–2561 (1998) 12. Du, B., Xu, X., Dai, X.: Minimum-phase FIR precoder design for multicasting over MIMO frequency- selectiv e channels. Journal of Electronics (China) 30 (4), 319–327 (2013) 13. Ephremidze, L., Janashia, G., Lagvilava, E.: A simple proof of the matrix-v alued Fej ´ er-Riesz theorem. Journal of Fourier Analysis and Applications 15 (1), 124–127 (2009) 14. Ephremidze, L., Saied, F ., Spitkovsk y , I.: On the algorithmization of Janashia-Lagvilava matrix spectral factorization method. Submitted to IEEE Trans. Inform. Theory 15. Fischer, R.F .: Sorted spectral factorization of matrix polynomials in MIMO communications. IEEE T ransactions on Communications 53 (6), 945–951 (2005) 16. Gan, L., Ma, K.K.: On minimal lattice factorizations of symmetric-antisymmetric multifilterbanks. IEEE T rans. Signal Process. 53 (2, part 1), 606–621 (2005) 17. Geronimo, J.S., Hardin, D.P ., Massopust, P .R.: Fractal functions and wa velet expansions based on several scaling functions. J. Approx. Theory 78 (3), 373–401 (1994) 18. Golub, G.H., V an Loan, C.F .: Matrix computations, fourth edn. Johns Hopkins Studies in the Mathemat- ical Sciences. Johns Hopkins Univ ersity Press, Baltimore, MD (2013) 19. Hansen, M., Christensen, L.P .B., Winther , O.: Computing the minimum-phase filter using the QL- factorization. IEEE Trans. Signal Process. 58 (6), 3195–3205 (2010) 20. Hardin, D.P ., Hogan, T .A., Sun, Q.: The matrix-valued Riesz lemma and local orthonormal bases in shift-in variant spaces. Advances in Computational Mathematics 20 (4), 367–384 (2004) 21. Hsung, T .C., Lun, D.P .K., Shum, Y .H., Ho, K.C.: Generalized discrete multiwavelet transform with em- bedded orthogonal symmetric prefilter bank. IEEE Trans. Signal Process. 55 (12), 5619–5629 (2007) 22. Hsung, T .C., Sun, M.C., Lun, D.K., Siu, W .C.: Symmetric prefilters for multiwa velets. IEE Proceedings- V ision, Image and Signal Processing 150 (1), 59–68 (2003) 23. Huo, G., Miao, L.: Cycle-slip detection of GPS carrier phase with methodology of SA4 multi-wa velet transform. Chinese Journal of Aeronautics 25 (2), 227–235 (2012) 24. Je ˇ zek, J., Ku ˇ cera, V .: Efficient algorithm for matrix spectral factorization. Automatica 21 (6), 663–669 (1985) 25. Jiang, Q.: On the regularity of matrix refinable functions. SIAM J. Math. Anal. 29 (5), 1157–1176 (1998) 26. Jiang, Q.: Orthogonal multiwav elets with optimum time-frequency resolution. IEEE T rans. Signal Pro- cess. 46 (4), 830–844 (1998) 27. Lawton, W .: Applications of complex valued wa velet transforms to subband decomposition. IEEE T rans- actions on Signal Processing 41 (12), 3566–3568 (1993) 28. Lebrun, J., V etterli, M.: Balanced multiwav elets theory and design. IEEE T rans. Signal Process. 46 (4), 1119–1125 (1998) 26 V asil K olev et al. 29. Lebrun, J., V etterli, M.: High-order balanced multiwa velets: theory , factorization, and design. IEEE T rans. Signal Process. 49 (9), 1918–1930 (2001) 30. Li, B., Peng, L.: Balanced multiwav elets with interpolatory property . IEEE T rans. Image Process. 20 (5), 1450–1457 (2011) 31. Li, Y .F ., Y ang, S.Z.: Construction of paraunitary symmetric matrices and parametrization of symmetric orthogonal multiwav elet filter banks. Acta Math. Sinica (Chin. Ser .) 53 (2), 279–290 (2010) 32. Massopust, P .R., Ruch, D.K., V an Fleet, P .J.: On the support properties of scaling v ectors. Appl. Comput. Harmon. Anal. 3 (3), 229–238 (1996) 33. Micchelli, C.A., Sauer, T .: Regularity of multiwav elets. Adv . Comput. Math. 7 (4), 455–545 (1997) 34. Plonka, G., Strela, V .: Construction of multiscaling functions with approximation and symmetry . SIAM J. Math. Anal. 29 (2), 481–510 (1998) 35. Roux, J.L.: 2D spectral factorization and stability test for 2D matrix polynomials based on the radon projection. In: Acoustics, Speech, and Signal Processing, IEEE International Conference on ICASSP ’86., vol. 11, pp. 1041–1044 (1986) 36. Sayed, A.H., Kailath, T .: A survey of spectral factorization methods. Numer . Linear Algebra Appl. 8 (6-7), 467–496 (2001). Numerical linear algebra techniques for control and signal processing 37. Shen, L., T an, H.H., Tham, J.Y .: Symmetric-antisymmetric orthonormal multiw avelets and related scalar wav elets. Appl. Comput. Harmon. Anal. 8 (3), 258–279 (2000) 38. Smith, M., Barnwell, T .: Exact reconstruction techniques for tree-structured subband coders. IEEE T ransactions on Acoustics, Speech, and Signal Processing 34 (3), 434–441 (1986) 39. Strela, V .: A note on construction of biorthogonal multi-scaling functions. In: W avelets, multiwav elets, and their applications (San Diego, CA, 1997), Contemp. Math. , v ol. 216, pp. 149–157. Amer . Math. Soc., Providence, RI (1998) 40. Strela, V ., Heller , P .N., Strang, G., T opiwala, P ., Heil, C.: The application of multiwavelet filterbanks to image processing. IEEE Transactions on Image Processing 8 (4), 548–563 (1999) 41. Strela, V ., W alden, A.: Signal and image denoising via wa velet thresholding: orthogonal and biorthog- onal, scalar and multiple wavelet transforms. In: Nonlinear and nonstationary signal processing (Cam- bridge, 1998), pp. 395–441. Cambridge Univ . Press, Cambridge (2000) 42. T an, H.H., Shen, L.X., Tham, J.Y .: New biorthogonal multiwav elets for image compression. Signal Processing 79 (1), 45–65 (1999) 43. Tham, J.Y ., Shen, L., Lee, S.L., T an, H.H.: Good multifilter properties: a new tool for understanding multiwav elets. In: Proc. Intern. Conf. on Imaging, Science, Systems and T echnology CISST -98, Las V egas, USA, pp. 52–59 (1998) 44. Tham, J.Y ., Shen, L., Lee, S.L., T an, H.H.: A general approach for analysis and application of discrete multiwav elet transforms. IEEE Trans. Signal Process. 48 (2), 457–464 (2000) 45. Turcajo v ´ a, R.: Hermite spline multiwavelets for image modeling. In: Proc. SPIE 3391, W avelet Appli- cations V , 46, Orlando, FL, USA, pp. 46–56 (1998) 46. W ang, Z., McWhirter, J.G., W eiss, S.: Multichannel spectral factorization algorithm using polynomial matrix eigenv alue decomposition. In: 2015 49th Asilomar Conference on Signals, Systems and Com- puters, pp. 1714–1718 (2015) 47. Wiener , N., Masani, P .: The prediction theory of multiv ariate stochastic processes. I. The regularity condition. Acta Math. 98 , 111–150 (1957) 48. Wilson, G.T .: The factorization of matricial spectral densities. SIAM J. Appl. Math. 23 , 420–426 (1972) 49. Wu, G., Li, D., Xiao, H., Liu, Z.: The M -band cardinal orthogonal scaling function. Appl. Math. Comput. 215 (9), 3271–3279 (2010) 50. Y oula, D.C., Kazanjian, N.N.: Bauer-type factorization of positive matrices and the theory of matrix polynomials orthogonal on the unit circle. IEEE Trans. Circuits and Systems CAS-25 (2), 57–69 (1978)