In this section we briefly review the theory of multiwavelets and multiwavelet filter banks, to a large extent based on the work in [9, 10]. We will restrict to orthogonal multiwavelets having compact support. Compact support translates into finite impulse response (FIR) filters, whereas orthogonality is captured conveniently when switching from a filter bank description to polyphase filtering: it corresponds to the polyphase FIR filter being para-unitary, i.e. lossless. For lossless filters parameterizations are available in the literature, which offer opportunities to build in additional properties. We will exploit these when addressing vanishing moments and balancing conditions.

A multiresolution structure for the space L 2 ( R ) is defined, as usual in wavelet theory, to consist of nested approximation spaces V _m , which are linear subspaces of L 2 ( R ) such that …, ⊂ V ₋₁ ⊂ V ₀ ⊂ V ₁ ⊂ …, with intersection ∩ m ∈ Z V m = { 0 } and with completeness ∪ m ∈ Z V m ̄ = L 2 ( R ) . For each m ∈ Z we define the detail space W _m to be such that V _m+1 = V _m ⊕ W _m and V _m ⊥ W _m . Thus, V _m and W _m are orthogonal complements within the enveloping space V _m+1. Furthermore, we assume that this structure is both shift-invariant ∀ f ( t ) ∈ L 2 ( R ) , k , m ∈ Z : f ( t ) ∈ V m ⇔ f ( t − k ) ∈ V m and scale-invariant ∀ f ( t ) ∈ L 2 ( R ) , m ∈ Z : f ( t ) ∈ V m ⇔ f ( 2 t ) ∈ V m + 1 .

In orthogonal multiwavelet theory, it is assumed next that there exists a multiscaling function. This is a vector of r scaling functions Φ ( t ) = ( ϕ [ 0 ] ( t ) , ϕ [ 1 ] ( t ) , … , ϕ [ r − 1 ] ( t ) ) T of which the entries, together with their integer translates, generate an orthonormal basis of the approximation space V ₀. Thus, for each such multiscaling function we have that ϕ [ j ] ( t − k ) , ϕ [ p ] ( t − ℓ ) = δ j , p δ k , ℓ for all integers k , l ∈ Z and all indices j, p ∈ {0, 1, …, r − 1}. Here δ _j,k is the Kronecker delta, equal to 1 if j = k and equal to 0 otherwise. In this set-up, the orthonormal basis of V ₀ is generated by the r scaling functions jointly rather than by any single one of them. Because of shift and scale-invariance of the multiresolution structure, the orthonormal basis of V ₀, for each m, induces an orthonormal basis of V _m which is constituted by the entries of the vectors in { 2 m / 2 Φ ( 2 m t − k ) | k ∈ Z } .

Likewise, it is assumed that there exists an associated multiwavelet function Ψ(t) which together with its integer translates generates an orthonormal basis of the detail space W ₀. This multiwavelet is a vector of r wavelet functions Ψ ( t ) = ( ψ [ 0 ] ( t ) , ψ [ 1 ] ( t ) , … , ψ [ r − 1 ] ( t ) ) T , with the property ψ [ j ] ( t − k ) , ψ [ p ] ( t − ℓ ) = δ j , p δ k , ℓ for all k , l ∈ Z and all j, p ∈ {0, 1, …, r − 1}. Since V ₀ ⊥ W ₀ it holds that ϕ [ j ] ( t − k ) , ψ [ p ] ( t − ℓ ) = 0 , for all j, p, k, and ℓ. As before for the spaces V ₀ and V _m , the orthonormal basis of W ₀ induces an orthonormal basis of W _m which consists of the entries of the vectors in { 2 m / 2 Ψ ( 2 m t − k ) | k ∈ Z } .

From the fact that Φ ( t ) ∈ V 0 r and V ₀ ⊂ V ₁, it follows that Φ(t) can be represented in terms of the induced basis 2 Φ ( 2 t − k ) | k ∈ Z of V ₁. So, there exist unique r × r matrix coefficients C _k , k ∈ Z , for which the multiscaling function Φ(t) satisfies a two-scale vector equation (a refinement equation) called the dilation equation: Φ t = 2 ∑ k = − ∞ ∞ C k Φ 2 t + k . (1)

Likewise, we have that Ψ ( t ) ∈ W 0 r and W ₀ ⊂ V ₁, so it can also be represented uniquely in terms of the same induced basis of V ₁. This implies that there exist unique r × r matrix coefficients D _k , k ∈ Z , which express the multiwavelet Ψ(t) in terms of the multiscaling function on V ₁ by means of a two-scale vector equation called the wavelet equation: Ψ t = 2 ∑ k = − ∞ ∞ D k Φ 2 t + k . (2)

If we also impose that the multiscaling and multiwavelet functions have compact support (i.e., they are nonzero only on a finite interval, which makes them localized in space—a desirable property of wavelets), it follows that only a finite number of basis functions in the induced basis of V ₁ can contribute to representing Φ(t) and Ψ(t). As a consequence, the infinite sums in the dilation and wavelet equation become finite sums, and by shifting Φ(t) and Ψ(t) appropriately, if necessary, it can be assumed without loss of generality that the index k in each sum runs from k = 0 to k = 2n − 1 for some integer n ≥ 1. ²

The finite matrix coefficient sequences {C _k |k = 0, 1, …2n − 1} and {D _k |k = 0, 1, …2n − 1} are used to define two r × r polynomial matrices in z ⁻¹ according to C ( z ) = ∑ k = 0 2 n − 1 C k z − k and D ( z ) = ∑ k = 0 2 n − 1 D k z − k . The function C(z) is the transfer function of a FIR low-pass multiwavelet filter, whereas D(z) is the transfer function of a FIR high-pass multiwavelet filter. Together they make up an orthogonal FIR filter bank, satisfying the Smith-Barnwell orthogonality conditions [29]: C z C z − 1 T + C − z C − z − 1 T = 2 I r , (3) D z D z − 1 T + D − z D − z − 1 T = 2 I r , (4) C z D z − 1 T + C − z D − z − 1 T = 0 . (5)

The converse of what has been presented so far, is also largely true: if one starts from an orthogonal FIR filter bank with r × r FIR filters C(z) and D(z) satisfying Eqs 3–5, then this induces a multiresolution structure with orthonormal translation invariant bases generated by the multiscaling and multiwavelet functions satisfying the dilation equation and wavelet equation under relatively mild conditions. What is needed, is what we have assumed earlier, namely that the dilation equation admits a solution and that the spaces eventually span all of L 2 ( R ) .

The cascade algorithm is a tool to compute Φ(t) for a given choice of filter coefficients C _k , k = 0, 1, …, 2n − 1 by iteration: from a current estimate of Φ(t) a new estimate is computed by substituting the current estimate into the right-hand side of the dilation equation and reading off the new estimate from the resulting left-hand side. This requires a suitable initialisation, which is achieved by the Haar multiwavelet, for which (consistent with our convention) the multiscaling function Φ(t) is given by ϕ [ i ] ( t ) = r on the interval ( − i + 1 r , − i r ] and ϕ _[i](t) = 0 elsewhere. Clearly, this Haar multiwavelet is an example of an orthogonal multiwavelet (albeit not a spectacular one, as it simply mimics the scalar Haar wavelet). It is easily verified that orthonormality of the induced basis is preserved in each iteration step and that compact support holds too due to the finite number of terms in the summation; so if the cascade algorithm converges it will produce a feasible solution for Φ(t). If the algorithm happens to not converge, then a multiwavelet interpretation is lacking, but the filter bank may still prove valuable from a signal processing perspective. The cascade algorithm is useful to visualize the multiwavelets which correspond to a chosen filter bank. When discussing balancing conditions, it is instrumental for determining the vector v ₀ in Eq. 39.

Assume a function (a scalar signal) s(t) ∈ V _m+1 to be represented by a sequence of r × 1 vector coefficients {s _ℓ } in terms of the induced basis of V _m+1: s t = ∑ ℓ ∈ Z 2 m + 1 / 2 s ℓ T Φ 2 m + 1 t − ℓ . (6)

Since V _m+1 = V _m ⊕ W _m we can decompose s(t) into s t = a t + b t , (7) in which a(t) is an approximation signal with coefficient vectors {a _ℓ } with respect to the induced multiscaling basis in V _m , and b(t) is a detail signal (orthogonal to a(t)) with coefficient vectors {b _ℓ } with respect to the induced multiwavelet basis in W _m : a t = ∑ ℓ ∈ Z 2 m / 2 a ℓ T Φ 2 m t − ℓ , (8) b t = ∑ ℓ ∈ Z 2 m / 2 b ℓ T Ψ 2 m t − ℓ . (9)

It now holds that these coefficient vectors can be computed from those of s(t) as follows: a ℓ = ∑ k = 0 2 n − 1 C k s 2 ℓ − k , (10) b ℓ = ∑ k = 0 2 n − 1 D k s 2 ℓ − k . (11)

In digital signal processing terms, the sequences {a _ℓ } and {b _ℓ } are obtained from {s _ℓ } by first filtering {s _ℓ } with the filters C(z) and D(z) (in parallel), and then dyadically down-sampling the resulting sequences: only the even-indexed coefficient vectors are kept (and relabeled), all the odd-indexed coefficient vectors are discarded. ³

This process can be reorganized in a more efficient way, by first splitting the coefficient vector sequence {s _ℓ } into two phases (using down-sampling on the sequence and on a delayed duplicate) and then passing both phases jointly through a suitable polyphase filter H _p (z) to directly produce {a _ℓ } and {b _ℓ }. More generally, if we start from a discrete-time representation of the scalar signal s(t) as a sequence of scalar values {f _ℓ } (for instance obtained by regular sampling) and corresponding z-transform S ( z ) = ∑ ℓ ∈ Z f ℓ z − ℓ , then one can split {f _ℓ } into 2r phases and submit these jointly to the 2r × 2r polyphase filter H _p (z). These 2r phases are defined by rewriting S(z) as: S z = S 0 z 2 r + z 1 S 1 z 2 r + ⋯ + z 2 r − 1 S 2 r - 1 z 2 r , (12) where S _k (z) denotes the z-transform of the k-th phase of {f _ℓ } (k = 0, 1, …, 2r − 1). This matches the situation with a coefficient vector sequence {s _ℓ } if each s _ℓ contains the phases 0, 1, …, r − 1 of the scalar coefficient sequence {f _ℓ } for even indices ℓ and the phases r, r + 1, …, 2r − 1 for odd indices ℓ.

The 2r × 2r FIR polyphase filter H _p (z) is constructed as: where C _even(z), C _odd(z), D _even(z), and D _odd(z) are r × r polynomial matrices (in z ⁻¹) that split C(z) and D(z) into two phases, such that C z = C even z 2 + z − 1 C odd z 2 (14) D z = D even z 2 + z − 1 D odd z 2 . (15)

This architecture is illustrated in Figure 1 and leads to the same approximation coefficients a _ℓ and detail coefficients b _ℓ as before, where we used double subscripts to anticipate a multiresolution structure with repeated down-sampling and filtering steps later on.

For the 2r × 2r polyphase filter H _p (z) defined in Eq. 13 the orthogonality conditions of Eqs 3–5 translate into the equivalent condition: H p z H p z − 1 T = I 2 r . (16)

This implies that H _p (z) is a lossless FIR filter, with coefficient matrices H k = C 2 k C 2 k + 1 D 2 k D 2 k + 1 , for k = 0, 1, …, n − 1.

Complementary to Eq. 13, which expresses H _p (z) in terms of C(z) and D(z), the low-pass and high-pass FIR filters C(z) and D(z) are reobtained from the polyphase filter H _p (z) by: C z D z = H p z 2 I r z − 1 I r . (17)

Lossless polyphase filters have a rich structure and have been well studied in the literature. See for instance [30, 31] for an accessible review of lossless filters and their key properties, including application areas and examples of the various uses they have. The basic property of a discrete-time lossless system, from which its name derives, is that the total energy of the output signals equals the total energy of the input signals, regardless of the signals being used. Here, energy is measured by the sum of squares of the values in the (discrete-time) signals. For H _p (z) this precise property is captured by Eq. 16. The point is that for |z| = 1 it holds that z ⁻¹ = z ^∗ (i.e., the complex conjugate) and since H _p (z) has real coefficients it follows that H p ( z − 1 ) T = H p ( z ) ∗ . This shows that H _p (z) is unitary for all z on the unit circle, which causes the conservation of energy property.

For the construction of arbitrary lossless FIR polyphase filters (for given choices of r and n), the condition of Eq. 16 is not convenient to work with and impose directly. Equating coefficients for all the entries of the corresponding expressions on the left and on the right, can of course be done but does not provide orthogonality conditions in a form suitable for analytic or numerical computation. Instead, there is a body of literature which describes how the class of all 2r × 2r lossless systems of a fixed given order can be parameterized, by carrying out a recursion with respect to the system order, known as the tangential Schur algorithm; see [8]. This approach is flexible, as it allows the user to build in certain properties: for instance, by making specific choices for some parameters in the general procedure it allows one to parameterize the subclass of lossless FIR filters.

The following theorem was used in [27, 28] to parameterize (real) lossless FIR polyphase filters of a given order; see also [10, 31] for a similar yet slightly different construction. Here, the order k of a lossless FIR polyphase filter G ^(k)(z) is the McMillan degree of this rational matrix. For lossless functions this equals the degree of the denominator of the rational function det(G ^(k)(z)). For the lossless FIR polyphase filter H _p (z), the order n − 1 traditionally is the maximum lag of the filter, which relates to the length of the interval on which the multiwavelet lives. In the scalar case r = 1 these definitions coincide. In the multiwavelet case r ≥ 2 they do not. Generically, the McMillan degree and the maximum lag are the same for the lossless FIR polyphase filters in our class, but for special choices of the unit vectors u _k they may be different: the maximum lag of H _p (z) is less than its McMillan degree n − 1 if and only if u k + 1 T u k = 0 for some k ∈ {1, …, n − 2}. We prefer to adopt the McMillan degree definition of the filter order when addressing parameterizations, to avoid having to consider special cases resulting from this mismatch. ⁴

Theorem 3.1. Let G ^(k)(z) be a real 2r × 2r lossless FIR polyphase filter of order k ≥ 1. Then there exists a real 2r × 1 vector u _k of norm ‖u _k ‖ = 1 such that G ^(k)(z) can be factored as: G k z = I 2 r + z − 1 − 1 u k u k T G k − 1 z , (18) where G^(k−1)(z) is a real 2r × 2r lossless FIR polyphase filter of order k − 1.

This theorem shows that a lossless FIR polyphase filter of order k can be reduced in k iteration steps to a lossless FIR polyphase filter G ⁽⁰⁾(z) of order 0, which is nothing else than a constant orthogonal matrix. We therefore drop the variable z for this last filter and simply write G ⁽⁰⁾.

For the polyphase matrix H _p (z) of the filter bank, the recursion makes clear that G ⁽⁰⁾ is in fact equal to the value of H _p (z) at z = 1: G ⁽⁰⁾ = H _p (1). This will make it more convenient to build in balanced vanishing moments than for the construction described in [10]. Generically, there will exist unit vectors u ₁, …, u _n−1 which allow H _p (z) to be factored as: H p z = I 2 r + z − 1 − 1 u n − 1 u n − 1 T … , I 2 r + z − 1 − 1 u 1 u 1 T G 0 . (19)

To make the parameterization explicit, all the vectors u _k of norm 1 as well as the orthogonal matrix G ⁽⁰⁾ still can be parameterized in terms of scalar parameters. This is easily achieved. A 2r × 1 unit vector u _k can be recursively parameterized as u k = cos ( θ k , 1 ) g k ⁡ sin ( θ k , 1 ) with θ _k,1 ∈ [0, π] and g _k again a unit vector but of smaller size (2r − 1) × 1. For orthogonal matrices it is well-known how they can be parameterized, for instance, with Givens rotations and Householder matrices. ⁵ When addressing balanced vanishing moments of several orders, this parameterization will be refined further.

When processing a signal with a multiwavelet filter bank, the differences in spectral behavior of the r components of the multiscaling function may result in unbalanced channels of the low-pass filter. The consequences will further escalate in a multiresolution structure, when the filter bank is repeatedly applied to part of the output (the approximation signals) of a previous filtering step. The imbalance will adversely affect performance in compression and detection applications in signal and image processing, due to the fact that high and low frequencies start to mix. This problem is known as the balancing problem and was pointed out and characterized in detail in [9, 12, 24, 32].

The balancing problem is first encountered when processing a single scale of a multiwavelet filter bank with a sampled univariate signal as in Figure 1. Assume that the sampled input signal with z-transform S(z) is constant, then all down-sampled signals with z-transforms S _k (z) will be constant too and in fact be equal. A constant signal is supposed to pass unaltered through the polyphase filter, as the multiwavelet function Ψ(t) should pick up details and the multiscaling function Φ(t) the trend (an approximation to the signal); indeed, a constant function is all about trends and has no details. This can be achieved by imposing a zero-order vanishing moment on the filter bank, as is commonly done. This is the multiwavelet counterpart to the admissibility condition for scalar wavelets: ∫ R Ψ ( t ) d t = 0 .

However, though each of the low-pass outputs a _1,ℓ (z) will then be constant (due to the vanishing moment), their values may differ between components, since the vanishing moment by itself does not enforce coherence between the low-pass output channels. The signals a _1,ℓ can therefore no longer be considered as phases from the constant approximation signal a(t). Though in itself this is already undesired as it hampers interpretation of the decomposition by the filter, this is further complicated when the signals are reused unaltered in a multiresolution structure. Spurious frequencies will then be introduced and severely compromise the usefulness and interpretability of the multiwavelet analysis. If a reconstruction of the approximation signal a(t) is done, possibly involving a number of upsampling and inverse filter steps, then this imbalance could be undone after which the signal can be split in 2r phases again as in Figure 1 for each wavelet scale. But to ensure that the filters will work properly on dyadic halfbands, this will involve extra effort whereas the interpretability of the approximation signals remains unfixed.

From a practical perspective, the approach as illustrated in Figure 2 is more appealing. Here, intermediate reconstruction and phase splitting between the scales is omitted. The approximation coefficient sequences are directly split into two phases and propagated to the next scale, as is common for scalar wavelets and consistent with the signal processing formulas (10)–(11). Without balancing, the consequence of this setup is that for a constant input signal at scale 2 the low-pass output a _2,ℓ will no longer even be constant. This happens despite the fact that the zero-order vanishing moment ensures that constant signals are retained in the approximation low-pass output. With balancing, a constant input signal produces equal output signals a _1,ℓ which remain well-behaved when further propagated to the next scale.

This issue generalizes to higher order vanishing moments which aim to have all polynomials up to a given degree in the approximation spaces. Including balanced vanishing moments of each order ≤ p ensures that every polynomial of degree ≤ p is annihilated in the high-pass output and retained as a polynomial of degree ≤ p in the low-pass output.

The p-th order vanishing moment condition on the multiwavelet function Ψ(t) is given by: ∫ R Ψ t t p d t = 0 . (20)

It makes clear that polynomials of degree ≤ p are fully suppressed in the detail spaces spanned by the multiwavelet function. The general form of the balancing conditions on the multiscaling function Φ(t) up to order p are given by [11]. In view of our convention for the dilation and wavelet equation they take the form: ∫ R ϕ 0 t t + 0 r p d t = ∫ R ϕ 1 t t + 1 r p d t = … = ∫ R ϕ r − 1 t t + r − 1 r p d t . (21)

Constructing specific orthogonal multiscaling functions which obey both (20) and (21) up to a given order p can be straightforward, e.g. by using a concatenation of r shifted versions of a Daubechies-p scaling function, where each ϕ _[j](t) has been shifted by − j/r. Then, due to orthogonality of the shifted Daubechies scaling functions, clearly an orthogonal multiscaling function is obtained, which mimics the scalar setup entirely by having the wavelet forms ψ _[j](t) now capture what is normally achieved by integer translations of a single wavelet. Using such Daubechies-p wavelets ensures vanishing moments up to and including order p − 1. Since each component ϕ _[j](t) of the multiscaling function is just a time-shifted version of the others, it is also balanced.

However, constructing non-trivial orthogonal multiwavelet functions with different (non-shifted) wave forms and balanced vanishing moments up to a given order p is not straightforward. To address this question, we proceed as follows. (1) First, we derive vanishing moment conditions on the polyphase filter H _p (z). (2) Next, we derive additional balancing conditions. (3) We investigate how to satisfy all these conditions, up to a chosen order p, by building them into the parameterization of lossless FIR polyphase filters of order n. This is achieved for orders 0 and 1, whereas for order 2 the conditions can be brought into a form which can be used for numerical search.

Let us introduce the r × 1 vectors v _k and w _k (for k = 0, 1, 2, …) by defining: v k = ∫ R Φ t t k d t , w k = ∫ R Ψ t t k d t . (22)

These are the k-th moments of the multiscaling and multiwavelet functions. From the dilation and wavelet Eqs 1–2, we derive the following relationships. First, integration gives v 0 = 1 2 2 ∑ k = 0 2 n − 1 C k v 0 , (23) w 0 = 1 2 2 ∑ k = 0 2 n − 1 D k v 0 . (24)

Here it is used that ∫ R Φ ( 2 t + k ) d t = 1 2 ∫ R Φ ( τ ) d τ by substituting τ = 2t + k.

Next, if we first multiply left and right hand sides of the dilation and wavelet equations by t and then integrate, we obtain in a similar fashion: v 1 = 1 4 2 ∑ k = 0 2 n − 1 C k v 1 − ∑ k = 0 2 n − 1 k C k v 0 , (25) w 1 = 1 4 2 ∑ k = 0 2 n − 1 D k v 1 − ∑ k = 0 2 n − 1 k D k v 0 . (26)

Here it is used that ∫ R Φ ( 2 t + k ) t d t = 1 4 ∫ R Φ ( τ ) τ d τ − k ∫ R Φ ( τ ) d τ by writing t = 1 2 ( ( 2 t + k ) − k ) and substituting τ = 2t + k.

More generally, if we first multiply left and right hand sides of the dilation and wavelet equations by t ² and then integrate, we obtain: v 2 = 1 8 2 ∑ k = 0 2 n − 1 C k v 2 − 2 ∑ k = 0 2 n − 1 k C k v 1 + ∑ k = 0 2 n − 1 k 2 C k v 0 , (27) w 2 = 1 8 2 ∑ k = 0 2 n − 1 D k v 2 − 2 ∑ k = 0 2 n − 1 k D k v 1 + ∑ k = 0 2 n − 1 k 2 D k v 0 . (28)

Here we wrote t 2 = 1 4 ( ( 2 t + k ) − k ) 2 and substituted τ = 2t + k to get ∫ R Φ ( 2 t + k ) t 2 d t = 1 8 ∫ R Φ ( τ ) τ 2 d τ − 2 k ∫ R Φ ( τ ) τ d τ + k 2 ∫ R Φ ( τ ) d τ . It is clear how this generalizes to arbitrary order p, using t p = 1 2 p ( τ − k ) p .

The vanishing moment conditions up to and including order p are: w 0 = 0 , w 1 = 0 , … , w p = 0 , (29) because w _k represents the (vector) wavelet coefficient of the polynomial signal t ^k for the untranslated multiwavelet Ψ(t), and these wavelet coefficients should all vanish. The vanishing moment conditions state that the filter coefficients C _k and D _k (k = 0, 1, …, 2n − 1) should be such as to allow these equations to be satisfied for some vectors v ₀, v ₁, …, v _p , representing the moments of the multiscaling function.

For orders p = 0, p = 1, and p = 2, the vanishing moment conditions above can be rewritten in terms of the filters C(z) and D(z). Focusing on the equations involving C(z) (the other equations with D(z) are handled entirely analogously) we see from differentiation, that C z = ∑ k = 0 2 n − 1 C k z − k , C ′ z = − ∑ k = 0 2 n − 1 k C k z − k + 1 , C ″ z = ∑ k = 0 2 n − 1 k k + 1 C k z − k + 2 , (30) whence ∑ k = 0 2 n − 1 C k = C 1 , ∑ k = 0 2 n − 1 k C k = − C ′ 1 , ∑ k = 0 2 n − 1 k 2 C k = C ″ 1 + C ′ 1 . (31)

The relationship (17) links the filters C(z) and D(z) to the polyphase filter H _p (z). Differentiation after multiplication by any constant vector v gives: C ′ z D ′ z v = 2 z H p ′ z 2 v z − 1 v − z − 2 H p z 2 0 v , (32) which produces at z = 1: C ′ 1 D ′ 1 v = 2 H p ′ 1 v v − H p 1 0 v . (33)

An additional differentiation step yields: C ″ z D ″ z v = 4 z 2 H p ″ z 2 v z − 1 v + 2 H p ′ z 2 v z − 1 v − 4 z − 1 H p ′ z 2 0 v + 2 z − 3 H p z 2 0 v , (34) which gives at z = 1: C ″ 1 D ″ 1 v = 4 H p ″ 1 v v + 2 H p ′ 1 v − v + 2 H p 1 0 v . (35)

These results are combined and summarized in the following theorem.

Theorem 3.2. Let the vectors v_k (k = 0, 1, 2) be defined as in Eq. 22. The polyphase filter H_p(z) imposes a vanishing moment of order 0 on the multiwavelet if it holds that: v 0 0 = 1 2 2 H p 1 v 0 v 0 . (36)

A vanishing moment of order 1 occurs if: v 1 0 = 1 4 2 2 H p ′ 1 v 0 v 0 + H p 1 v 1 v 1 − v 0 . (37)

A vanishing moment of order 2 occurs if: v 2 0 = 1 8 2 ( 4 H p ″ 1 v 0 v 0 + 4 H p ′ 1 v 0 + v 1 v 1 + H p 1 v 2 v 2 − 2 v 1 + v 0 ) . (38)

To have balanced vanishing moments of order up to and including p, requires extra conditions on the vectors v_k, k = 0, 1, 2, …, p. This is addressed next.

The balancing conditions can be computed from Eq. 21 by integration.

For order p = 0 this gives that all entries of v ₀ are equal, say with value α ₀. To find α ₀, the cascade algorithm is useful. For the Haar multiwavelet in the initialisation step, the scaling function ϕ _[j](t) is positive and constant on the interval ( − j + 1 r , j r ] of width 1 r ; hence, as the basis is orthonormal, this constant value equals r . The value ( v 0 ) [ j ] therefore initially equals 1 r . At the initialisation step, the vector v ₀ already has the required form and therefore it is preserved during the cascade algorithm iterations. It follows that α 0 = 1 r , and so: v 0 = 1 r 1 ⋮ 1 . (39)

For order p = 1 it is obtained that ( v 1 ) [ j ] + j r ( v 0 ) [ j ] should have the same value, say λ r for all j = 0, 1, …, r − 1. Imposing balancing of order 1 only makes sense when balancing of order 0 is also imposed. Then, with the earlier result ( v 0 ) [ j ] = 1 r , it is obtained that: v 1 = − 1 r r 0 1 ⋮ r − 1 + λ r 1 1 ⋮ 1 . (40)

By inspecting what happens for balanced orthogonal multiwavelets generated from scalar orthogonal wavelets with vanishing moments through shifting, just as explained earlier for the Haar multiwavelet, we find that the value of λ varies between multiwavelets. It therefore is left as a free parameter.

Next, for order p = 2 it is obtained that ( v 2 ) [ j ] + 2 j r ( v 1 ) [ j ] + j r 2 ( v 0 ) [ j ] should have the same value, say μ r for all j = 0, 1, …, r − 1. With the earlier results ( v 0 ) [ j ] = 1 r and ( v 1 ) [ j ] = − j r r + λ r we obtain: v 2 = 1 r 2 r 0 1 2 ⋮ r − 1 2 − 2 λ r r 0 1 ⋮ r − 1 + μ r 1 1 ⋮ 1 (41)

The parameter μ is again left free.

From a signal processing point of view, if a signal s(t) is a polynomial in t of degree ≤ p and passed through the multiwavelet filter bank starting from its sampled sequence, then the output channels retain the interpretation of phases of sampled signals.

These additional balancing conditions on the vectors v ₀, v ₁ and v ₂ can now be combined with the previous vanishing moment conditions. This gives the following result.

Theorem 3.3. (a) The polyphase filter H_p(z) imposes a balanced vanishing moment of order 0 on the multiwavelet structure if it holds that:

(b) Balanced vanishing moments of orders 0 and 1 occur if in addition to the condition of (a) it holds that there exists a constant λ such that:

(c) Balanced vanishing moments of orders 0, 1 and 2 occur if in addition to the conditions of (a) and (b) it holds that there exists a constant μ such that:

The proof of this result is by direct computation, using Eqs 39–41 together with Thm. 3.2. Part (a) of this theorem was shown earlier in the work of [9] and used before in [27]. Parts (b) and (c) are novel characterizations, using the lossless polyphase filter and its derivatives at z = 1.

A parametrization of lossless FIR polyphase filters with a balanced vanishing moment of order 0 has previously been described in [27, 28]. We present it here for completeness and also because some of the techniques and notation will be reused when we address balancing of order 1. In fact, the parameterization in Eq. 19 is slightly different from the parameterization used in the literature cited above, but the basic ideas are similar.

From part (a) of Thm. 3.3 we have that for a zero-order balanced vanishing moment, the orthogonal matrix H _p (1) must map the vector 1 … 1 1 … 1 T into the vector 2 1 … 1 0 … 0 T . This can be achieved with Householder transformation matrices. Let R ₁ be the Householder matrix which maps 1 … 1 0 … 0 T to r 1 0 … 0 0 … 0 T and let R ₂ be the Householder matrix which maps 1 … 1 1 … 1 T to 2 r 1 0 … 0 0 … 0 T . A Householder matrix R _v,w which maps a vector v into another vector w with ‖w‖ = ‖v‖ is orthogonal and symmetric and of the form R v , w = I − 2 w − v w − v T w − v T w − v . (45)

Clearly R v , w = R w , v = R v , w T = R v , w − 1 . For the matrices R ₁ and R ₂ we have the explicit forms

Note that the zero-order condition can be rewritten as and in this way turns into

The matrix R ₁ H _p (1)R ₂ is again orthogonal, and it therefore necessarily is of the form with Q an orthogonal matrix of size (2r − 1) × (2r − 1). The converse clearly also holds true, which motivates the following result.

Theorem 3.4. All lossless FIR polyphase filters H_p(z) of order n − 1 with a balanced vanishing moment of order 0 are obtained as: where R ₁ and R ₂ are the fixed Householder matrices given in Eqs 46, 47, where Q ranges over the set of (2r − 1) × (2r − 1) orthogonal matrices, and where u ₁, u ₂, …, u _n−1 ranges over the set of n − 1 unit vectors of size 2r × 1.

Explicit parameterization of the orthogonal matrix Q can be done with Givens rotations and Householder matrices as indicated earlier, see Section 2.2. There we also showed how a parameterization of the unit vectors u ₁, …, u _n−1 can be obtained. We see that manipulating the vectors u _k (or increasing their number) does not affect the zero order balancing property as it is fully implied by the structure of Eq. 51. This is conveniently exploited when considering an additional balanced vanishing moment of order 1 below.

To build an extra balanced vanishing moment of order 1 into the lossless FIR polyphase filters H _p (z), retaining the balanced vanishing moment of order 0, the idea is to start from the parameterization in Thm. 3.4 and to refine the structure of the orthogonal matrix Q and the unit vectors u ₁, …, u _n−1 to meet the condition of Eq. 43. The latter condition involves an additional scalar parameter λ, which can be freed up by premultiplication of all the terms in the equation by the Householder matrix R ₁ of Eq. 46 as we will now show. In view of the structure of H _p (1) = G ⁽⁰⁾ given by Eq. 50, it is useful to first work out the following two matrix-vector products (the horizontal lines are just for clarity):

The matrix H p ′ ( 1 ) results from differentiation of Eq. 19 followed by substitution of z = 1, which gives:

Using all of this, condition (43) takes the form:

The parameter λ only appears in the top row of this vector equation. Selecting this row admits computation of λ as: λ = − 2 ∑ k = 1 n − 1 R 1 u k 1 2 − 1 2 r . (55)

Here, ( R 1 u k ) 1 denotes the first entry in the unit vector R ₁ u _k . It is convenient to parameterize these vectors R ₁ u _k (rather than the unit vectors u _k ) as R 1 u k = cos θ k u ̃ k ⁡ sin θ k (56) for some scalar θ _k ∈ [0, π] and with u ̃ k a unit vector of size (2r − 1) × 1. Then λ = − 2 ∑ k = 1 n − 1 cos 2 θ k − 1 2 r , (57) whereas, upon division by 2 r , the remaining part of the balancing condition attains the form

Introducing the following (r − 1) × 1 vectors h _r , for positive integers r: h r = 2 r r 1 ⋮ r − 1 − 1 + r r 1 ⋮ 1 (59) allows us to summarize these findings concisely in the following theorem.

Theorem 3.5. All lossless FIR polyphase filters H _p (z) of order n − 1 with two balanced vanishing moments, of orders 0 and 1, are obtained as: where R ₁ and R ₂ are the fixed Householder matrices given in Eqs 46, 47, where Q is (2r − 1) × (2r − 1) orthogonal, and where u _k (k = 1, …, n − 1) are unit vectors of size 2r × 1 parameterized as u k = R 1 cos θ k u ̃ k ⁡ sin θ k , (61) with scalar θ _k ∈ [0, π], unit vectors u ̃ k of size (2r − 1) × 1, and such that the following condition is satisfied: with the vectors h _r and h _2r defined as in Eq. 59.

To complete the parameterization of Theorem 3.5, we now show how all tuples of orthogonal Q, unit vectors u ̃ k and scalars θ _k can be obtained which make up all the solutions of Eq. 62.

For h_r it is straightforward to compute its norm ‖h_r‖ as ‖ h r ‖ = 1 3 − 1 3 r 2 . (63)

For r = 1, we encounter the scalar orthogonal wavelet case, and the vector h _r is in fact empty (of size 0, ×, 1). Balancing is not meaningful then. Indeed the condition reduces to the scalar equation − 1 2 Q = − ∑ k = 1 n − 1 u ̃ k ⁡ sin ( 2 θ k ) in which Q as well as u ̃ 1 , … , u ̃ n − 1 are all ±1. Note that the scalars u ̃ k can all be fixed to 1, because u _k and therefore also u ̃ k need only be parameterized up to a sign. For Q it is remarked in general that the space of real orthogonal matrices has two connected components, characterized by the determinant which is either 1 or − 1. For r = 1, restricting to just one of these components corresponds to using a sign convention for ψ(t), which can be performed by fixing a sign for the last row of H _p (z) and therefore the sign of Q. Choosing Q = 1 and taking slight differences in parameterization into account, the resulting condition ∑ k = 1 n − 1 ⁡ sin ( 2 θ k ) = 1 2 is entirely consistent, unsurprisingly, with the vanishing moment condition of order 1 reported in [6, 10, 28] or in Proposition 2 in [7].

Focusing on the multiwavelet case with r ≥ 2, it is noted that ‖h _r ‖ increases monotonically from 1 2 = 0.50 (at r = 2) to 1 3 ≈ 0.58 (for r → ∞). It is also noted that ‖ Q h 2 r ‖ = ‖ h 2 r ‖ = 1 3 − 1 12 r 2 because Q is orthogonal.

Writing g k = − u ̃ k ⁡ sin ( 2 θ k ) , we see that this is an arbitrary vector of norm ≤ 1 , with its direction determined by the unit vector u ̃ k and its norm determined by   sin(2θ _k ). With a sum of m such vectors g _k one can precisely cover all points in a (2r − 1)-dimensional hyperball of radius m centered at the origin.

The idea of the parameterization, is to build the vectors g _k one by one, such that the condition remains feasible, i.e., the remaining vectors-to-be-constructed can still be assigned values to meet Eq. 62. Suppose g ₁, …, g _ℓ−1 have been chosen and denote (for ℓ = 1, 2, …, n):

Then, for ℓ = 1, 2, …, n − 2, the vector g _ℓ must be chosen such that the vector q _ℓ+1 = q _ℓ + g _ℓ is in (or on) the hyperball (centered at the origin) of radius r ℓ ≔ ‖ h 2 r ‖ + n − 1 − ℓ . (65)

This still allows the right-hand side of Eq. 62, viz. the expression q _ℓ+1 + g _ℓ+1 + ⋯ + g _n−1 = q _n , to eventually land on the hypersphere of vectors with norm ‖h _2r ‖. For the left-hand side expression Qh _2r , note that the norm equals ‖h _2r ‖ for every orthogonal Q, and that a suitable orthogonal matrix Q can always be constructed to make Qh _2r equal to that landing point (a Householder matrix would do). In fact, all such orthogonal Q can be constructed along the same lines as before when constructing G ⁽⁰⁾ for balancing of order 0, as we will demonstrate below in Eq. 74. Note that if one cannot land on that hypersphere the construction is infeasible, since the norm of Qh _2r is fixed; so this precisely characterizes feasibility. (Note, also, that points inside the hypersphere of radius ‖h _2r ‖ can always be brought to a point on the hypersphere in one step, by adding a single vector g _k , since ‖h _2r ‖ < 1; therefore one need not consider a constraint on the inside of the hypersphere.)

For ℓ = 1, 2, …, n − 2, we proceed as follows. If ‖q _ℓ ‖ ≤ r _ℓ − 1, then g _ℓ can be any vector of norm ≤ 1 , so there are no particular constraints on u ̃ ℓ (of norm 1) or θ _ℓ . Note that r ℓ − 1 = ‖ h 2 r ‖ + n − 2 − ℓ ≥ ‖ h 2 r ‖ ≥ 1 2 , so this is a non-empty set of vectors including the case q _ℓ = 0.

If ‖q _ℓ ‖ > r _ℓ − 1, then let q ℓ ⊥ denote an arbitrary vector of unit norm and orthogonal to q _ℓ . The space of such vectors is easily parameterized. ⁶ Next, note that each vector g _ℓ can be written uniquely as: g ℓ = α ℓ q ℓ + β ℓ q ℓ ⊥ , (66) with scalar coefficients α _ℓ (still unrestricted) and β _ℓ ≥ 0. For g _ℓ to be feasible it must hold that: (1) ‖g _ℓ ‖ ≤ 1, (2) ‖q _ℓ + g _ℓ ‖ ≤ r _ℓ . Because q _ℓ and q ℓ ⊥ are orthogonal to each other, setting β _ℓ = 0 gives the maximum range for α _ℓ . Requirement (1) implies: α ℓ ∈ [ − 1 ‖ q ℓ ‖ , 1 ‖ q ℓ ‖ ] . Requirement (2) implies: α ℓ ∈ [ − 1 − r ℓ ‖ q ℓ ‖ , − 1 + r ℓ ‖ q ℓ ‖ ] . Combining these two intervals, using ‖q _ℓ ‖ > r _ℓ − 1, it follows that both requirements are met if and only if: − 1 ‖ q ℓ ‖ ≤ α ℓ ≤ − 1 + r ℓ ‖ q ℓ ‖ . (67)

Once α _ℓ is chosen, β _ℓ can be chosen non-negative but such that the two conditions remain satisfied. This gives: 0 ≤ β ℓ ≤ min 1 − α ℓ 2 ‖ q ℓ ‖ 2 , r ℓ 2 − 1 + α ℓ 2 ‖ q ℓ ‖ 2 . (68)

For ℓ = n − 1, choosing the final vector g _n−1 (of norm ≤ 1 ) is special, in the sense that it must satisfy the equality ‖q _n−1 + g _n−1‖ = r _n−1 = ‖h _2r ‖ rather than an inequality as before.

If ‖q _n−1‖ ≤ 1 − r _n−1, then q _n can land on every point of the hypersphere with radius r _n−1 by choosing g _n−1 as the difference between that point and q _n−1.

If ‖q _n−1‖ > 1 − r _n−1, then the hypersphere centered at the origin with radius r _n−1 and the hypersphere centered at q _n−1 with radius 1 intersect at points q n = q n − 1 + g n − 1 = ( 1 + α n − 1 ) q n − 1 + β n − 1 q n − 1 ⊥ if the following two conditions hold: α n − 1 2 ‖ q n − 1 ‖ 2 + β n − 1 2 = 1 , (69) 1 + α n − 1 2 ‖ q n − 1 ‖ 2 + β n − 1 2 = r n − 1 2 . (70)

Subtracting the first equation from the second leaves a linear equation for α _n−1, which gives: α n − 1 = − 1 2 1 + 1 − r n − 1 2 ‖ q n − 1 ‖ 2 . (71)

This determines the minimal value for α _n−1 to give solutions on the hypersphere for q _n . The maximal value for α _n−1 is obtained with β _n−1 = 0 and equals α n − 1 = − 1 + r n − 1 ‖ q n − 1 ‖ . This implies that all feasible g n − 1 = α n − 1 q n − 1 + β n − 1 q n − 1 ⊥ have − 1 2 1 + 1 − r n − 1 2 ‖ q n − 1 ‖ 2 ≤ α n − 1 ≤ − 1 + r n − 1 ‖ q n − 1 ‖ (72) and β n − 1 = r n − 1 2 − 1 + α n − 1 2 ‖ q n − 1 ‖ 2 . (73)

Finally, to find all feasible orthogonal Q once the right-hand side expression q _n is constructed to have ‖q _n ‖ = r _n−1 = ‖h _2r ‖, let R ₃ be the Householder matrix which maps q _n to ‖ h 2 r ‖ 1 0 … 0 T and let R ₄ be the Householder matrix which maps h _2r to ‖ h 2 r ‖ 1 0 … 0 T . Then all feasible orthogonal matrices Q are obtained as with Q ̃ an arbitrary orthogonal matrix of size (2r − 2) × (2r − 2).

Incorporating an additional balanced vanishing moment of order p = 2, requires further refinement of the current parameterization such that condition (44) is also satisfied. As we have not developed a way to incorporate this into a parameterization, we will only show how the condition can be reworked into a form which makes it suitable to be added as a constraint to the parametrization in terms of the parameters used before.

A first step to achieve this, is to free up and remove the parameter μ by premultiplication of all terms in the condition by the Householder matrix R ₁. The effect of this is that μ only appears in the term ( 3 μ − 2 λ 2 ) r 2 2 R 1 1 1 … 1 0 … 0 T = ( 3 μ − 2 λ 2 ) r 2 2 r 1 0 … 0 0 … 0 T , of which just the first entry is nonzero. Therefore, as before for λ for balancing of order 1, the parameter μ can be left free and computed afterwards from the equation which corresponds to the first entries on the left and right-hand sides.

The actual conditions for balancing of order 2, are then captured by the remaining 2r − 1 equations in which μ no longer appears. The parameter λ which does show up in them, should be replaced by its value in terms of the chosen parameters as given in Eq. 57. This gives the constraints in a form which can readily be used for numerical optimization, using a routine for constrained optimization which admits nonlinear constraints.

For the scalar orthogonal wavelet design case, sparsity of a prototype signal was already advocated in [7]. This is appealing for detection and compression purposes. There it was also shown that maximization of sparsity of the vector of all approximation and detail coefficients w = {w _k } in the setting with a critically sampled orthogonal wavelet transform, boils down to either: 1) maximization of the variance of the sequence of absolute values |w _k | of the coefficients, or 2) maximization of the variance of the squared wavelet coefficients. This is due to the fact that Parseval’s identity holds [4], which makes that the sum of squares of all the wavelet and detail coefficients equals the sum of squares of all the values of the digital signal being processed. This result does obviously not change for orthogonal multiwavelet design. Hence [[7], Theorem 1], holds in the current setting too:

Theorem 3.6. Let w be the vector of the approximation coefficients at the coarsest scale and the detail coefficients at all scales, resulting from the processing of a signal s by means of an orthogonal multiwavelet filter bank across multiple scales. Then:

1) Maximization of the variance of the vector of absolute values {|w _k |} is equivalent to minimization of the L ₁-norm of the vector w.

For practical applications the latter criterion, L ₄-maximization, is more appealing than the former since it can be combined with weighing both in scale and location [27]. This allows for designing a matched multiwavelet, where each of the components of Φ(t) is forced to focus on a single aspect in a signal (see Section 3.2 for an example).

If instead of a critically sampled multiwavelet transform an undecimated multiwavelet transform [33] is used, this has advantages for detection purposes, since due to the redundant representation time-invariance is achieved. The lack of down-sampling, however, causes an abundance of coefficients at coarser scales. To preserve conservation of energy properties for orthogonal multiwavelets, this abundance can be counteracted by dyadic discounting of energies towards coarser scales as in [7]. This makes it possible to use adapted versions of the L ₁-minimization and L ₄-maximization criteria for undecimated multiwavelet transforms promoting sparsity.

Another criterion one might consider is entropy. One can quantify the effectiveness of compression algorithms in terms of entropy [34]. As suggested in [35], entropy can be used to select an orthogonal wavelet basis, and a similar approach can serve as a design criterion. This is certainly appealing for compression purposes.

The effects of having unbalanced vanishing moments are easily illustrated by considering what happens if a constant signal is fed into an unbalanced orthogonal multiwavelet filter bank with multiplicity r = 2 and of order n − 1 = 3, using the construction scheme in Theorem 3.1, but without balancing. To this end we took Q = I _2r−1 = I ₃ and we changed the Householder matrix R ₁ to produce a non-constant vector v ₀.

When a constant input signal is fed into an orthogonal multiwavelet filter bank with an order 0 unbalanced vanishing moment, as described in Section 3.1, the results (depending on the exact choice of coefficients) are as shown in Figure 3. What can be observed in the top figure, is that the detail coefficients are all zero for both channels (consistent with the imposed vanishing moments), and the approximation coefficients are constant per channel, but the values are different between the channels. In terms of the notation in Eqs 8, 9, the detail vector coefficients b _k are all zero, the approximation vector coefficients a _k are all constant, but their entries ( a k ) 1 and ( a k ) 2 are different.

In the bottom figure of Figure 3 it is displayed what the resulting signal looks like if the outputs are nevertheless interpreted as the two phases of a signal.

If the outputs are now processed further with the multiwavelet filter bank for a few more scales, this has the effect that at later scales the constant nature of the original signal is lost and the detail coefficients become nonzero, as demonstrated in Figure 4.

In this section we introduce an example for balanced multiwavelet design. We will employ L ₄-maximization with weighted masking to match the multiwavelet filter bank to a prototype signal. The example evolves around creating a representation for a prototype ECG signal as is commonly encountered in the field of cardiology. The goal of this design is for feature detection to distinguish and detect different complexes in the ECG signal.

Feature detection is an application that makes efficient use of the different components of the multiwavelet. The idea of feature detection via multiwavelets is to design each component of the multiwavelet to detect a specific feature in the signal. Orthogonality between the components of the scaling- and wavelet function allows the multiwavelet to detect up to r orthogonal features, which cannot all be accurately detected by a scalar wavelet at the same time. Since features in a signal need not be orthogonal but can be overlapping, making them hard to separate, the multiwavelet approach is fundamentally different from template based approaches. When training orthogonal multiwavelets on different features, the goal is to have them pick up orthogonal aspects from the features that help to distinguish between them.

For this application the variance maximizing L ₄-objective function is no longer measured across all the wavelet coefficients. Instead, a time-scale mask is employed for each component of the multiwavelet [28], and the value of the objective function of a component is only measured on the wavelet coefficients that are in the time-scale mask. In this way, maximization of the criterion forces energy into the masked areas, helping the multiwavelets to focus on the user-selected signal features covered by the time-scale masks. (This same idea is less conveniently pursued with L ₁-minimization, because minimization will force energy out of the masked areas, but while the energy leaks away this does not mean that useful multiwavelets are promoted.)

In this example, the design procedure employs a first-order balanced orthogonal multiwavelet, with focus on the approximation coefficients rather than the detail coefficients, and aiming to detect features by looking at the peaks of the approximation coefficients at masked areas at the coarsest scale l = 3.

It can be seen from Figure 5 that the first scaling function mainly represents the QRS-complex on which it was trained, whereas the second scaling function manages to capture the T-wave on which it was trained. By thresholding the approximation coefficients of the first channel the location of the QRS-complex can be obtained. By thresholding the second channel the location of the T-peak is obtained. Both channels can detect the feature that they were designed for independently from one another. This is remarkable, as the T-wave overlaps in spectrum and waveform with the QRS-complex, which carries far more energy. Though the second wavelet also picks up some energy from the QRS-complex, it still allows to detect the T-peak location with simple thresholding. This demonstrates that a matched multiwavelet is a powerful signal processing tool to detect and discriminate between different features in a signal.

In this example we show an orthogonal multiwavelet filter with multiplicity r = 3, involving a lossless FIR polyphase filter H _p (z) of order n − 1 = 3, with balanced vanishing moments of orders p = 0, p = 1, and p = 2 as a concrete example of multiwavelet design with three balanced vanishing moments. To obtain such a filter, we utilized the parameterization that was introduced in Section 2.7 for multiwavelets with balanced vanishing moments of orders 0 and 1, and we followed the approach of Section 2.8. For r = 3 and n − 1 = 3, we first chose random parameters to obtain an initial orthogonal multiwavelet only lacking a second order vanishing moment. From there, constrained optimization was used to adapt the parameters to satisfy condition (44), producing an orthogonal multiwavelet with three balanced vanishing moments.

For this example, the FIR polyphase filter H _p (z) = H ₀ + H ₁ z ⁻¹ + H ₂ z ⁻² + H ₃ z ⁻³ has the following coefficient matrices H _k : H 0 = 0.5698 0.8165 0.0101 0.0248 0.0105 − 0.0342 − 0.0317 − 0.0070 0.7024 0.6980 0.1093 − 0.0555 0.0167 0.0051 − 0.0687 0.0488 0.5488 0.8164 − 0.0351 0.1002 0.1538 − 0.2415 0.0724 − 0.0026 0.0577 − 0.0455 − 0.4025 0.2809 0.6175 − 0.3845 0.0170 0.0280 − 0.0517 0.0043 0.1349 − 0.0645 (75) H 1 = 0.0551 − 0.0351 0.0089 − 0.0026 − 0.0170 0.0190 − 0.0205 0.0109 0.0101 − 0.0165 0.0270 − 0.0225 0.1108 − 0.0792 0.0464 − 0.0487 0.0339 − 0.0205 − 0.4528 0.2714 0.2526 − 0.3512 0.4382 − 0.3106 − 0.3402 0.2485 − 0.0335 0.0410 − 0.0668 0.0795 − 0.0069 − 0.0331 − 0.4161 0.3877 − 0.0384 − 0.1749 (76) H 2 = − 0.0283 0.0203 − 0.0104 0.0109 − 0.0080 0.0053 0.0266 − 0.0187 0.0077 − 0.0076 0.0042 − 0.0024 0.0136 − 0.0092 0.0020 − 0.0012 − 0.0011 0.0013 0.3095 − 0.2110 0.0525 − 0.0384 − 0.0127 0.0188 − 0.1124 0.0831 − 0.0559 0.0630 − 0.0558 0.0394 0.4275 − 0.3234 0.2539 − 0.2942 0.2794 − 0.2014 (77) H 3 = − 0.0042 0.0029 − 0.0011 0.0010 − 0.0004 0.0002 0.0012 − 0.0008 0.0003 − 0.0003 0.0001 − 0.0001 − 0.0020 0.0014 − 0.0005 0.0005 − 0.0002 0.0001 − 0.0360 0.0251 − 0.0093 0.0086 − 0.0037 0.0017 − 0.0362 0.0252 − 0.0093 0.0087 − 0.0037 0.0018 0.1926 − 0.1343 0.0495 − 0.0462 0.0196 − 0.0093 (78)

Associated with this polyphase filter are the wavelet and scaling functions as shown in Figure 6.

By definition, the Householder matrices R ₁ and R ₂ are given by: R 1 = 1 3 1 3 1 3 0 0 0 1 3 1 2 1 − 1 3 − 1 2 1 + 1 3 0 0 0 1 3 − 1 2 1 + 1 3 1 2 1 − 1 3 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 = 0.5774 0.5774 0.5774 0 0 0 0.5774 0.2113 − 0.7887 0 0 0 0.5774 − 0.7887 0.2113 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 (79) R 2 = 1 6 1 6 1 6 1 6 1 6 1 6 1 6 1 5 4 − 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 1 6 − 1 5 1 + 1 6 1 5 4 − 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 1 5 4 − 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 1 5 4 − 1 6 − 1 5 1 + 1 6 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 − 1 5 1 + 1 6 1 5 4 − 1 6 = = 0.4082 0.4082 0.4082 0.4082 0.4082 0.4082 0.4082 0.7184 − 0.2816 − 0.2816 − 0.2816 − 0.2816 0.4082 − 0.2816 0.7184 − 0.2816 − 0.2816 − 0.2816 0.4082 − 0.2816 − 0.2816 0.7184 − 0.2816 − 0.2816 0.4082 − 0.2816 − 0.2816 − 0.2816 0.7184 − 0.2816 0.4082 − 0.2816 − 0.2816 − 0.2816 − 0.2816 0.7184 (80)