%% %% This is a LaTeX document generated by automatic conversion of a %% Mathematica notebook using Mathematica 4.0. %% %% This document uses special macros that are defined in a LaTeX %% package file with a name of the form notebook*.sty. Appropriate %% commands for loading the package have been included in this file. %% The LaTeX package files can be found in the Mathematica %% installation under %% %% Tigger:Applications:Mathematica 4.0:SystemFiles:IncludeFiles:TeX: %% %% The LaTeX package used in in this document may require that your %% LaTeX implementation support fonts other than the standard Computer %% Modern fonts that are used by LaTeX. See the Additional Information %% section under the TeXSave entry in the Mathematica reference guide %% for more information. %% \documentclass{article} \usepackage{notebook2e, latexsym} \begin{document} \Title{Automation of the lifting factorisation of wavelet transforms} \Subtitle{M. Maslen, P. Abbott\inlineTFinmath{{{}^{a}}}} \noindent \inlineTFinmath{{{}^a}} Department of Physics, University of Western Australia, Nedlands WA 6907, Australia \\ E-mail: {\tt paul@physics.uwa.edu.au} \\ URL: {\verb$http://physics.uwa.edu.au/~paul$} \\ \Section*{Abstract} Wavelets are sets of basis functions used in the analysis of signals and images. In contrast to Fourier analysis, wavelets have both spatial and frequency localisation, making them useful for the analysis of sharply-varying or non-periodic signals. The {\itshape lifting scheme} for finding the discrete wavelet transform was demonstrated by Daubechies and Sweldens (1996). In particular, they showed that this method depends on the factorisation of {\itshape polyphase matrices}, whose entries are Laurent polynomials, using the Euclidean algorithm extended to Laurent polynomials. Such factorisation is not unique and hence there are multiple factorisations of the polyphase matrix. In this paper we outline a {\itshape Mathematica} program that finds all factorisations of such matrices by automating the Euclidean algorithm for Laurent polynomials. Polynomial reduction using Gr\ODoubleDot{}bner bases was also incorporated into the program so as to reduce the number of wavelet filter coefficients appearing in a given expression through use of the relations they satisfy, thus permitting exact symbolic factorisations for any polyphase matrix.\\ \noindent PACS: 02.70, 07.05.K.\\ \noindent Keywords: Wavelets; Lifting; Euclidean algorithm; Laurent polynomials, \\ Gr\ODoubleDot{}bner bases, Polynomial reduction \\ \Subtitle{PROGRAM SUMMARY} \noindent Title of program: {\verb$LiftingFactorisation.nb$} \noindent Version number: 1.0 \noindent Catalogue identifier: \noindent Program obtainable from: {\tt ftp://ftp.physics.uwa.edu.au/-\\ pub/Wavelets/LiftingFactorisation.nb} \noindent Programming language used: {\itshape Mathematica} 3.0 \noindent Platform: Any platform supporting {\itshape Mathematica} 3.0 \noindent Number of bytes in distributed program, including test data, etc.: 36,764 \noindent Distribution format: ASCII \noindent Keywords: wavelets, lifting, Euclidean algorithm, Laurent polynomials, Gr\ODoubleDot{}bner bases, polynomial reduction \noindent Nature of physical problem \\ Spectral analysis and compression of signals or images. \noindent Method of solution \\ Symbolic computer algebra is used to automate the Euclidean algorithm for Laurent polynomials [1] so as to factorise wavelet transforms yielding a sequence of simple arithmetic operations suitable for parallel, in-place, implementation [2]. \noindent Limitations \\ The program requires that the Laurent polynomial quotients used in the algorithm have Laurent degree at most 1. The polyphase matrix must have unit determinant (though any polyphase matrix may be adjusted to satisfy this criterion [2]). \noindent Unusual features \\ The program can find symbolic greatest common divisors (gcds) of two Laurent polynomials. Using Gr\ODoubleDot{}bner bases and polynomial reduction on the filter coefficients reduces the number of unknown coefficients appearing in expressions. There is also capability to convert the lifting steps from mathematical notation to computer pseudo-code, suitable for implementation in C or Fortran. \\ \noindent {\itshape References} \noindent [1] M.J. Maslen, Factoring Wavelet Transforms into Lifting Steps, (University of Western Australia, 1997). URL: \\ {\verb$http://www.physics.uwa.edu.au/~paul/publications.html$} \noindent [2] I. Daubechies, W. Sweldens, J. Fourier Anal. Appl. 4 (1998) 247. URL: \\ {\tt http://cm.bell-labs.com/who/wim/papers/papers.html\#{}factor} \\ \Subtitle{LONG WRITE-UP} \Section*{1. Introduction} Scientists and engineers are familiar with Fourier analysis, which is an ideal tool for describing information such as waves, optical signals, and vibrations. It is the periodic nature of the Fourier basis which makes it so useful for these applications, for which the information is frequently oscillatory in nature. However, for highly non-periodic signals, the Fourier expansion requires many terms to achieve acceptable representation. In such cases an alternative is to use {\itshape wavelets}. There are a variety of wavelet bases, each having particular characteristics such as symmetry, smoothness or orthogonality. As such they are far more versatile than Fourier analysis, and can be tailored to specific situations. Wavelets bases of {\itshape compact support} have excellent spatial localisation (in contrast with the Fourier basis). Wavelets have been applied to areas such as geophysical survey analysis, image compression (FBI fingerprint database) [3], denoising, and encryption. Further details can be found in, for example, [1,4,5,6,7,8]. For a survey of applications of the wavelet transform to a wide range of physical fields --- including astrophysics, turbulence, meteorology, plasma physics, atomic and solid state physics, and mathematical physics --- see [9]. Wavelet research has proliferated over the last decade. Much of the relevant background theory has been developed, and the recent focus has been on applications and implementation techniques. One particularly powerful implementation method is that of lifting. This has been shown to be equivalent to the more traditional methods such as Mallat's algorithm and Fourier methods [10,11,12,13]. Lifting has several advantages over other methods. It is a very general technique, applicable to {\itshape n}-dimensional transforms and wavelets on manifolds. In some cases the method can yield integer-to-integer transforms [14] allowing for lossless compression. Computationally it allows a parallel, in-place implementation, which is faster by a factor of 2 for large datasets [2]. It reduces the implementation to a sequence of simple steps involving elementary algebraic operations, making it trivial to find the inverse transform. A major recent development is the {\itshape factorisation} of the lifting scheme [2]. This provides a simple method for deriving the sequence of lifting steps. Importantly, the factorisation is fundamentally non-unique, so there exist several possible sequences of lifting steps for a given wavelet transform. This is an important feature as it allows a choice of which set of steps to use -- some sets may be better suited to hardware implementation, numerical computation, or symmetric implementation, for example. Two questions that arise from this factorisation are (i) how many factorisations exist, and (ii) how does one find all factorisations. This paper addresses both of these issues, in the form of a {\itshape Mathematica} program that finds all gcds of two Laurent polynomials. \Section*{2. Mathematical preliminaries} We outline some of the background mathematical ideas relating to wavelets and the lifting scheme. \Subsection*{2.1. Laurent polynomials, finite filters, and {\itshape z}--transforms} A{\itshape Laurent polynomial} is an expression of the form \inlineTFinmath {p(z) = \sum _{k=a}^{b}{c_k}{z^{-k}}.} \noindent We define the {\itshape degree} of a nonzero Laurent polynomial as \inlineTFinmath{\LeftBracketingBar p(z)\RightBracketingBar = b-a}, and the degree of the zero polynomial is defined to be -\inlineTFinmath{\infty }. A {\itshape finite filter} may be described as a set of coefficients \inlineTFinmath{h=\{{h_k} : k \in \DoubleStruckCapitalZ \}}, which are zero outside the finite range \inlineTFinmath{{k_b} \leq k \leq {k_e}}. The \inlineTFinmath{z}{\itshape -transform} of a finite filter is the Laurent polynomial whose coefficients are the filter coefficients, namely \dispTFNumberedEquationmath{h(z) = \sum _{k={k_b}}^{{k_e}}{h_k}{z^{-k}}.} \noindent The limiting points \inlineTFinmath{{k_e}} and \inlineTFinmath{{k_b}} are arbitrary subject to the restriction that \inlineTFinmath{{k_e}-{k_b}} is the length of the filter. Equivalently, the {\itshape z}-transform is defined only up to a monomial factor. \Subsection*{2.2. The Euclidean algorithm} The {\itshape Euclidean algorithm} is a method for finding the greatest common divisor (gcd) of two integers. The {\itshape division algorithm} states that for any integers {\itshape a} and {\itshape b}, where \inlineTFinmath{\LeftBracketingBar a\RightBracketingBar >\LeftBracketingBar b\RightBracketingBar } (here \(\LeftBracketingBar \)\(\cdot \)\(\RightBracketingBar \) denotes absolute value) and \inlineTFinmath{b\neq 0}, there exist integers \inlineTFinmath{{q_1}}and\inlineTFinmath{ {r_1}} (the reason for the subscripts will become clear) such that \inlineTFinmath{a = {q_1}b + {r_1}} with \inlineTFinmath{0 \leq {r_1} < b}. The Euclidean algorithm can be used to find the gcd as follows. Firstly note that if \inlineTFinmath{{r_1}=0}, then \inlineTFinmath{\gcd(a,b)=b}. If \inlineTFinmath{{r_1}\neq 0}, we see that any common divisor of {\itshape b} and \inlineTFinmath{{r_1}} is a divisor of {\itshape a}, and thus is a common divisor of {\itshape a} and {\itshape b}. Rewriting the above expression as \inlineTFinmath{{r_1} = a - {q_1}b}, shows that any common factor of {\itshape a} and {\itshape b} is a factor of \inlineTFinmath{{r_1}}, and is thus a common factor of \inlineTFinmath{{r_1}} and {\itshape b}. We have shown that the common divisors of {\itshape a} and {\itshape b} are the same as those of \inlineTFinmath{{r_1}} and {\itshape b}, and hence \inlineTFinmath{\gcd(a,b) = \gcd({r_1},b)}. The problem has thus been reduced to finding the gcd of the smaller pair of numbers \inlineTFinmath{({r_1},b)}. We may now apply the division algorithm to \inlineTFinmath{{r_1}} and {\itshape b} to get \inlineTFinmath{b = {q_2}{r_1} + {r_2}} with \inlineTFinmath{0 \leq {r_2} < {r_1}}, and repeating the argument above, we see that \inlineTFinmath{\gcd({r_1},b) = \gcd({r_1}, {r_2})}. We iterate this process, each time obtaining a smaller remainder until eventually \inlineTFinmath{{r_i}=0} for some {\itshape i}. We then have \inlineTFinmath{{r_{i-2 }}= {q_i}{r_{i-1}}} and thus we have \inlineTFinmath{\gcd(a,b)=\gcd({r_{i-2}},{r_{i-1}})={r_{i-1}}}. The algorithm has a simple computational implementation, in recursive form: \dispTFNumberedEquationmath{\MathBegin{MathArray}{l} {a_i}\leftarrow {b_{i-1}}, \\ \noalign{\vspace{0.5ex}}{b_i}\leftarrow {a_{i-1}}-{a_i} {q_i}.\\ \MathEnd{MathArray}} The Euclidean algorithm can be generalised to the Laurent polynomials. The procedure is similar to that shown above, with the generalisation being that we require that \inlineTFinmath{\LeftBracketingBar r\RightBracketingBar <\LeftBracketingBar b\RightBracketingBar }, where \inlineTFinmath{\LeftBracketingBar \cdot \RightBracketingBar } is the Laurent degree. \Subsection*{2.3. An example} An important point is that the quotients are {\itshape not} necessarily unique in the Euclidean algorithm. That is, although the gcd is unique up to an invertible (monomial) factor, the {\itshape quotients} can be chosen in several nontrivially distinct ways. To illustrate this we give an example of the Euclidean algorithm for Laurent polynomials. Let \inlineTFinmath{a=\frac{1}{z}+3-2 z}, and \inlineTFinmath{b=\frac{-6}{{z^2}}+\frac{3}{z}}. We wish to find \inlineTFinmath{\gcd(a,b)}. Note that \inlineTFinmath{\LeftBracketingBar a\RightBracketingBar =1-(-1)=2}, while \inlineTFinmath{\LeftBracketingBar b\RightBracketingBar =-1-(-2)=1}. We set up the first step of the Euclidean algorithm: \inlineTFinmath {\frac{1}{z}+3-2 z={q_1} \big(\frac{-6}{{z^2}}+\frac{3}{z}\big)+{r_1}\cdot } \noindent Since the remainder must have degree less than that of {\itshape b}, it must be monomial. Rearranging this equation for \inlineTFinmath{{r_1}}, we get \inlineTFinmath {{r_1}=\big(\frac{1}{z}+3-2 z\big)-{q_1} \big(\frac{-6}{{z^2}}+\frac{3}{z}\big).} \noindent We must choose \inlineTFinmath{{q_1}} so that \inlineTFinmath{{r_1}} will be monomial. One way to do this is to choose \inlineTFinmath{{q_1}} so that the highest monomials in each term in this expression are equal, and will thus cancel. This means that \inlineTFinmath{{q_1}} will need to be of the form \inlineTFinmath{{z^2} \big(c+\frac{d}{z}\big)} in order for the terms to have the same highest exponents. The parameters {\itshape c} and {\itshape d} are then found by the method of undetermined coefficients. Substituting, we find that the first remainder has the form \inlineTFinmath {z (-3 c-2)+6 c-3 d+\frac{6 d+1}{z}+3,} \noindent and, setting the coefficients of the {\itshape z} and constant terms to zero, we get \inlineTFinmath{c=-\frac{2}{3}, d=-\frac{1}{3}}, so our first quotient is \inlineTFinmath{-\frac{2 {z^2}}{3}-\frac{z}{3}}, and the first remainder is \inlineTFinmath{-\frac{1}{z}}. We now apply the Euclidean algorithm to this remainder and \inlineTFinmath{b} to get \inlineTFinmath {b={q_2} {r_1}+{r_2} \hspace{1em}\Rightarrow \hspace{1em} {r_2}=b-{q_2} {r_1}.} \noindent Once again, the remainder \inlineTFinmath{{r_2}} must have degree less than that of \inlineTFinmath{{r_1}}. Since \inlineTFinmath{{r_1}} has degree zero, it follows that \inlineTFinmath{{r_2}} must be zero (recall that the degree of the zero polynomial is defined as -\(\infty \)). As before, we choose a quotient form so that both powers are matched, and find that the second remainder has the form \inlineTFinmath {\frac{c+3}{z}+\frac{d-6}{{z^2}}\cdot } \noindent We can set both terms, and thus the remainder, to zero to get \inlineTFinmath{c=\Mvariable{--3}, d=6}, so the second quotient is \inlineTFinmath{\frac{6}{z}-3}. Alternatively, we could have decided to eliminate the lowest two powers from the original expression, or even the lowest and the highest power. This gives three possible factorisations, each with different quotients. The gcd here is defined only up to a power of {\itshape z }: the gcd is 'greatest' in the sense that it has greatest possible Laurent degree. In this case, we found that the gcd was \inlineTFinmath{\frac{-1}{z}} (the last nonzero remainder). The other factorisation branches lead to a monomial gcd also. In this case the original Laurent polynomial {\itshape a} can be retrieved via \inlineTFinmath{a={q_1}b+\kappa }, where \(\kappa \) is the gcd, and the reader can easily verify that the three factorisations below are all valid: \inlineTFinmath {a=\bigg(-\frac{2 {z^2}}{3}-\frac{z}{3}\bigg)b-\frac{1}{z},} \inlineTFinmath {a=\bigg(-\frac{2 {z^2}}{3}-\frac{z}{6}\bigg)b-\frac{1}{2},} \inlineTFinmath {a=\bigg(-\frac{7 {z^2}}{12}-\frac{z}{6}\bigg)b-\frac{z}{4}\cdot } The gcd is the same for each factorisation up to a monomial factor, but the first quotients differ nontrivially. The last quotients differ only by monomial factors, because at this stage we have no freedom over which powers to eliminate. Note that multiplying \inlineTFinmath{a} by \inlineTFinmath{z} and \inlineTFinmath{b} by \inlineTFinmath{{z^2}} converts these Laurent polynomials into (ordinary) polynomials in \inlineTFinmath{z}. Polynomial division then gives \inlineTFinmath {(a z)=-2 {z^2}+3 z+1=\big(-\frac{2 z}{3}-\frac{1}{3}\big)(b {z^2})-1,} \noindent which is equivalent to the first factorisation above, providing a simple computational check. Similarly, the last factorisation can be obtained by dividing \inlineTFinmath{a/z} by \inlineTFinmath{b z} (since these are ordinary polynomials in \inlineTFinmath{1/z}), {\itshape i.e.}, \inlineTFinmath {\frac{a}{z}=-2+\frac{3}{z}+\frac{1}{{z^2}}=\big(-\frac{7}{12}-\frac{1}{6 z}\big)(b z)-\frac{1}{4}\cdot } \noindent Furthermore, there are fast algorithms for polynomial division and these can be used for these computations. However, the second factorisation cannot be obtained in this way. Later, we will factorise certain \inlineTFinmath{2\times 2} matrices whose entries are Laurent polynomials (specifically {\itshape z}-transforms as defined above) by executing the Euclidean algorithm on two of the entries. The factors are matrices whose entries are determined by the quotients resulting from the algorithm. Thus, using the different sets of quotients, we can obtain several different matrix factorisations. \Subsection*{2.4. Gr\ODoubleDot{}bner bases and polynomial reduction} The Euclidean algorithm can yield complicated algebraic expressions when carried out on {\itshape z}-transforms with symbolic filter coefficients. It is desirable to reduce such expressions as much as possible, using any relations that exist between the filter coefficients. A useful way of carrying out such simplifications is to use {\itshape Gr\ODoubleDot{}bner bases} and {\itshape polynomial reduction}. {\itshape Mathematica} has built-in functions to generate Gr\ODoubleDot{}bner bases and to carry out polynomial reduction. The formal definition of a Gr\ODoubleDot{}bner basis is somewhat involved, making use of several concepts in the theory of rings and multivariate polynomials (the interested reader is referred to [1]). The essential point is that Gr\ODoubleDot{}bner bases encode the information contained in the constraints on the filter coefficients, which can then be exploited to simplify factorisations by the technique of polynomial reduction. The utility of Gr\ODoubleDot{}bner bases may be demonstrated by considering an example. In the general case of D{\itshape N} orthogonal wavelets ({\itshape D} denotes Daubechies and {\itshape N} here denotes the number of non-zero filter coefficients), the following relations are satisfied by the filter coefficients \inlineTFinmath{{h_k}} [5]: \inlineTFinmath {{{\sum }_k}{h_k} = {\sqrt{2}}}, \inlineTFinmath {{{\sum }_k}{h_k}{h_{k-2m}} = 2 {{\delta }_{0,m}}} for \inlineTFinmath {m=0,1,2,..., \frac{N}{2}-1}, \inlineTFinmath {{{\sum }_k}{{(-1)}^k}{k^{p-1}}{h_k}=0} where \inlineTFinmath {p=1,2,...,\frac{N}{2}}. These conditions arise from applying certain approximation and orthonormality constraints. The equations may be solved for the coefficients of a given D{\itshape N} basis. In general the above conditions give \inlineTFinmath{N+1} equations in \inlineTFinmath{N} unknowns, however the system of equations is consistent; square normalisation can be derived from the other conditions. They can be solved exactly when \inlineTFinmath{N\leq 6}, but for \inlineTFinmath{N>6} they must be solved numerically. Solving for \inlineTFinmath{N=6} gives \inlineTFinmath {\MathBegin{MathArray}[c]{ll} {h_0}=\frac{{\sqrt{2}}}{32} \big(1+{\sqrt{10}}+{\sqrt{5+2 {\sqrt{10}}}}\big),&{h_1}=\frac{{\sqrt{2}}}{32} \big(5+{\sqrt{10}}+3 {\sqrt{5+2 {\sqrt{10}}}}\big), \\ {h_2}=\frac{{\sqrt{2}}}{32} \big(10-2 {\sqrt{10}}+2 {\sqrt{5+2 {\sqrt{10}}}}\big),&{h_3}=\frac{{\sqrt{2}}}{32} \big(10-2 {\sqrt{10}}-2 {\sqrt{5+2 {\sqrt{10}}}}\big), \\ {h_4}=\frac{{\sqrt{2}}}{32} \big(5+{\sqrt{10}}-3 {\sqrt{5+2 {\sqrt{10}}}}\big),&{h_5}=\frac{{\sqrt{2}}}{32} \big(1+{\sqrt{10}}-{\sqrt{5+2 {\sqrt{10}}}}\big). \MathEnd{MathArray}} A Gr\ODoubleDot{}bner basis is a set of polynomial expressions which are constructed using the zero conditions satisfied by the coefficients. An important feature is that the resulting Gr\ODoubleDot{}bner basis always has exactly the same collection of common roots as the original set of polynomials and is a canonical form from which many properties can conveniently be deduced. {\itshape Mathematica} has a built-in command for computing a Gr\ODoubleDot{}bner basis: \dispTFinmath{\MathBegin{MathArray}{l} g=\Mfunction{GroebnerBasis}\big[\big\{\sum _{i=0}^{5}{h_i}-{\sqrt{2}},\sum _{i=0}^{5}{{(-1)}^i} {h_i},\sum _{i=0}^{5}i {{(-1)}^i} {h_i}, \\ \noalign{\vspace{1.59375ex}} \hspace{3.em} \sum _{i=0}^{5}{i^2} {{(-1)}^i} {h_i},\sum _{i=0}^{3}{h_i} {h_{i+2}},\sum _{i=0}^{1}{h_i} {h_{i+4}}\big\},\Mfunction{Table}[{h_i},\{i,0,5\}]\big]\\ \MathEnd{MathArray}} \dispTFoutmath{\MathBegin{MathArray}{l} \big\{65536 h_{5}^{4}-8192 {\sqrt{2}} h_{5}^{3}-3072 h_{5}^{2}-96 {\sqrt{2}} {h_5}+9, \\ \noalign{\vspace{0.625ex}} \hspace{1.em} -256 h_{5}^{3}+32 {\sqrt{2}} h_{5}^{2}+6 {h_5}+3 {h_4}, \\ \noalign{\vspace{0.625ex}} \hspace{1.em} -4096 h_{5}^{3}+512 {\sqrt{2}} h_{5}^{2}+192 {h_5}+24 {h_3}-3 {\sqrt{2}}, \\ \noalign{\vspace{0.625ex}} \hspace{1.em} 8 {h_2}+16 {h_5}-3 {\sqrt{2}},4096 h_{5}^{3}-512 {\sqrt{2}} h_{5}^{2}-168 {h_5}+24 {h_1}-9 {\sqrt{2}}, \\ \noalign{\vspace{0.625ex}} \hspace{1.em} 2048 h_{5}^{3}-256 {\sqrt{2}} h_{5}^{2}-96 {h_5}+24 {h_0}-3 {\sqrt{2}}\big\}\\ \MathEnd{MathArray}} \noindent This basis will eliminate \inlineTFinmath{{h_0},{h_1},{h_2},{h_3}}, and \inlineTFinmath{{h_4}} from any expression, replacing it with terms in \inlineTFinmath{{h_5}}. Now consider the following expression: \dispTFinmath{\Mvariable{expr}={h_1} {h_5}+{h_0} {h_3}+{h_4} {h_2};} This can be reduced to an expression involving only a single coefficient, here \inlineTFinmath{{h_5}}, using polynomial reduction with the above Gr\ODoubleDot{}bner basis: \dispTFinmath{\Mvariable{reduced}=\Mfunction{PolynomialReduce}[\Mvariable{expr},g,\Mfunction{Table}[{h_i},\{i,0,5\}]]//\Mfunction{Last}} \dispTFoutmath{\frac{64}{3} {\sqrt{2}} h_{5}^{3}-\frac{25 h_{5}^{2}}{3}-\frac{5 {h_5}}{4 {\sqrt{2}}}-\frac{1}{64}} Observe that the form of these filter coefficients makes it difficult to give this reduced form of \inlineTFinmath{\Mvariable{expr}} using elementary algebra. Note also that we do not need to know the exact values of all of the coefficients, only their defining relationships. For example, it can be shown that finding the exact form of the D8 coefficients is equivalent to solving an irreducible sixth degree polynomial. However, polynomial reduction can still be carried out on the D8 coefficients, giving rise to {\itshape exact} expressions containing only {\itshape one} coefficient. This may be more satisfactory for numerical implementations, as we only need to substitute one numerical value into exact expressions. \Subsection*{2.5. Mallat's algorithm} We briefly describe here an established algorithm for computing the discrete wavelet transform. For a sampled signal, considered for simplicity to be of length \inlineTFinmath{l={2^n}} (this is easily generalized), we can use the relations satisfied by the wavelets to compute the wavelet expansion coefficients. We define matrices called {\itshape quadrature mirror filters} (QMF) in terms of the filter coefficients, and multiply by our data vector to produce a step of the wavelet transform. After each matrix multiplication, half of the vector is 'left alone' -- it is one component of the wavelet transform. Before the next step, a new QMF matrix of half the original dimensions is constructed. Eventually we reach a point where we are operating on a two component vector (only if \inlineTFinmath{l={2^n}}), and the result of this step is the last component of the discrete wavelet transform. This is called {\itshape Mallat's algorithm }[5]. More precisely, we denote our sampled signal by the vector \inlineTFinmath{{s^{[n]}}}; these sampled values correspond to the highest resolution level. With the coefficients \inlineTFinmath{{h_k}} (\inlineTFinmath{{g_k}}) forming a low(high)-pass filter, we define the QMF matrices as \dispTFNumberedEquationmath{{{(H)}_{k,l}}={h_{l-2 k}}, {{(G)}_{k,l}}={g_{l-2 k}}.} \noindent Then we have the relations; \dispTFNumberedEquationmath{{s^{[j]}}=H {s^{[j+1]}}, {d^{[j]}}=G {s^{[j+1]}},} \noindent where the vectors \inlineTFinmath{{s^{[j]}}} and \inlineTFinmath{{d^{[j]}}} are the 'signal' and 'difference' terms at resolution {\itshape j}. In other words, the components of these vectors are the coefficients in the expansion of a signal into the wavelet basis. This defines a step of the discrete wavelet transform in terms of multiplication by {\itshape H} and {\itshape G}, the QMF matrices. The (in the \inlineTFinmath{l={2^n}} case) one-component vector \inlineTFinmath{{s^{[0]}}} is the average or 'DC' term of our signal. After the first step the QMF matrices are redefined at half the dimensions, and we carry out the step (4) on the signal term found at the previous level. Thus we are able to find the expansion coefficients at {\itshape all} levels beginning with our sampled values \inlineTFinmath{{s^{[n]}}}. As a simple example of this, consider an explicit matrix implementation of the D4 wavelet transform. It is, of course, highly inefficient to implement Mallat's algorithm using matrix multiplication: this implementation is presented here only to introduce some notation and indicate the algorithm's inherent parallelisability. For an initial data vector \inlineTFinmath{{s^{[3]}}=\{{s_0},{d_0},{s_1},{d_1},{s_2},{d_2},{s_3},{d_3}\}} (the reason for the artificial structure of this vector will become clear later), the first step, \inlineTFinmath {\Big(\MathBegin{MathArray}[c]{c} {s^{[2]}} \\ {d^{[2]}} \MathEnd{MathArray}\Big)=\Bigg(\MathBegin{MathArray}[c]{cccccccc} {h_0}&{h_1}&{h_2}&{h_3}&0&0&0&0 \\ 0&0&{h_0}&{h_1}&{h_2}&{h_3}&0&0 \\ 0&0&0&0&{h_0}&{h_1}&{h_2}&{h_3} \\ {h_2}&{h_3}&0&0&0&0&{h_0}&{h_1} \\ {g_2}&{g_3}&0&0&0&0&{g_0}&{g_1} \\ {g_0}&{g_1}&{g_2}&{g_3}&0&0&0&0 \\ 0&0&{g_0}&{g_1}&{g_2}&{g_3}&0&0 \\ 0&0&0&0&{g_0}&{g_1}&{g_2}&{g_3} \MathEnd{MathArray}\Bigg).\Bigg(\MathBegin{MathArray}[c]{c} {s_0} \\ {d_0} \\ {s_1} \\ {d_1} \\ {s_2} \\ {d_2} \\ {s_3} \\ {d_3} \MathEnd{MathArray}\Bigg),} \noindent can be split into the parallel operations involving only {\itshape h}; \dispTFNumberedEquationmath{\MathBegin{MathArray}{l} {s^{[2]}}\leftarrow \Bigg(\MathBegin{MathArray}[c]{cccc} {h_0}&{h_2}&0&0 \\ 0&{h_0}&{h_2}&0 \\ 0&0&{h_0}&{h_2} \\ {h_2}&0&0&{h_0} \MathEnd{MathArray}\Bigg).\Bigg(\MathBegin{MathArray}[c]{c} {s_0} \\ {s_1} \\ {s_2} \\ {s_3} \MathEnd{MathArray}\Bigg)+ \\ \noalign{\vspace{3.16667ex}} \hspace{2.em} \Bigg(\MathBegin{MathArray}[c]{cccc} {h_1}&{h_3}&0&0 \\ 0&{h_1}&{h_3}&0 \\ 0&0&{h_1}&{h_3} \\ {h_3}&0&0&{h_1} \MathEnd{MathArray}\Bigg).\Bigg(\MathBegin{MathArray}[c]{c} {d_0} \\ {d_1} \\ {d_2} \\ {d_3} \MathEnd{MathArray}\Bigg)\\ \MathEnd{MathArray},} \noindent and {\itshape g}: \dispTFNumberedEquationmath{\MathBegin{MathArray}{l} {d^{[2]}}\leftarrow \Bigg(\MathBegin{MathArray}[c]{cccc} {g_2}&0&0&{g_0} \\ {g_0}&{g_2}&0&0 \\ 0&{g_0}&{g_2}&0 \\ 0&0&{g_0}&{g_2} \MathEnd{MathArray}\Bigg).\Bigg(\MathBegin{MathArray}[c]{c} {s_0} \\ {s_1} \\ {s_2} \\ {s_3} \MathEnd{MathArray}\Bigg)+ \\ \noalign{\vspace{3.16667ex}} \hspace{2.em} \Bigg(\MathBegin{MathArray}[c]{cccc} {g_3}&0&0&{g_1} \\ {g_1}&{g_3}&0&0 \\ 0&{g_1}&{g_3}&0 \\ 0&0&{g_1}&{g_3} \MathEnd{MathArray}\Bigg).\Bigg(\MathBegin{MathArray}[c]{c} {d_0} \\ {d_1} \\ {d_2} \\ {d_3} \MathEnd{MathArray}\Bigg)\\ \MathEnd{MathArray}.} The {\itshape lifting scheme} described in \S{}3 presents an alternative, highly efficient, method for computing the wavelet transform --- only involving simple multiplication and addition operations --- which is faster than Mallat's algorithm by a factor of 2 for large datasets [2]. Moreover, lifting is more general, applicable to {\itshape n}-dimensional transforms and wavelets on manifolds. \Section*{3. Lifting} In this section we outline the general theory of the {\itshape lifting scheme}. Further details can be found in{\itshape }[2], and broader information in [10,11,12,13]. The lifting scheme has been shown to be equivalent to the traditional method of implementing the wavelet transform via the QMF matrices [2]. It has several advantages over Mallat's algorithm, as will be shown below. It turns out that the procedure is equivalent to the factorisation of a matrix whose components are related to the {\itshape z}-transform of the filter. This matrix factorisation is non-unique, and we can find all factorisations through automation of the Euclidean algorithm for Laurent polynomials. \Subsection*{3.1. Polyphase representations} We define the {\itshape polyphase representation} of a filter; \dispTFNumberedEquationmath{h(z)={h_e}({z^2})+{z^{-1}}{h_o}({z^2}),} \noindent where {\itshape h} is the {\itshape z}-transform of the filter as defined in equation (1). Here \inlineTFinmath{{h_e}} (\inlineTFinmath{{h_o}}) is the {\itshape even (odd) polyphase component of h}, defined by \dispTFNumberedEquationmath{{h_e}({z^2})=\frac{h(z)+h(-z)}{2}, {h_o}({z^2})=\frac{h(z)-h(-z)}{2 {z^{-1}}}\cdot } \noindent We see that \inlineTFinmath{{h_e}} and \inlineTFinmath{{h_o}} contain the filter coefficients at even and odd locations respectively. We now define the {\itshape polyphase matrix}: \dispTFNumberedEquationmath{P(z)=\big(\MathBegin{MathArray}[c]{cc} {h_e}(z)&{g_e}(z) \\ {h_o}(z)&{g_o}(z) \MathEnd{MathArray}\big),} \noindent where \inlineTFinmath{{g_e}} and \inlineTFinmath{{g_o}} are defined in terms of the high pass filter {\itshape g} in analogously to (7). Note that the determinant of this matrix is always monomial [2]. The program requires it to be unity, but any polyphase matrix can be suitably adjusted to satisfy this requirement [2]. Here {\itshape z} is not a value to be substituted, but a symbolic entity representing a data shift. In general for a filter {\itshape f}, we have the rule: \dispTFNumberedEquationmath{{z^m} {f_l}\rightarrow {f_{l-m}},} \noindent where the direction of the shift is a choice of convention. The polyphase matrix provides an alternative method of executing the lifting scheme. For a data vector \inlineTFinmath{{s^{[j+1]}}}, we begin by performing the {\itshape lazy wavelet transform}, which simply involves subsampling the data at even and odd locations. That is; \dispTFNumberedEquationmath{\MathBegin{MathArray}{l} {{{s^{[j]}}}_l}\leftarrow {{{s^{[j+1]}}}_{2 l}}, \\ \noalign{\vspace{0.604167ex}}{{{d^{[j]}}}_l}\leftarrow {{{s^{[j+1]}}}_{2 l+1}}.\\ \MathEnd{MathArray}} In analogy with Mallat's algorithm, The vector \inlineTFinmath{({s^{[j]}},{d^{[j]}})} then premultiplies the polyphase matrix to give one stage of the discrete wavelet transform, bearing in mind that the components \inlineTFinmath{{s^{[j]}}}{\itshape }and \inlineTFinmath{{d^{[j]}}} are {\itshape themselves} vectors, and the polyphase matrix elements are Laurent polynomials containing the {\itshape z}-operator (shift). That is, a general step may be written as \dispTFNumberedEquationmath{({s^{[j]}},{d^{[j]}})\leftarrow ({s^{[j+1]}},{d^{[j+1]}}).\big(\MathBegin{MathArray}[c]{cc} {h_e}(z)&{g_e}(z) \\ {h_o}(z)&{g_o}(z) \MathEnd{MathArray}\big).} \noindent The result of the multiplication is the signal and detail at the next level, and the next step is to repeat the lazy wavelet transform (11) on \inlineTFinmath{{s^{[j]}}} and proceed as in (12). The length of \inlineTFinmath{{s^{[j]}}} and \inlineTFinmath{{d^{[j]}}} is reduced by a factor of 2 at each step, and thus the algorithm must always terminate after {\itshape n} steps. \Subsection*{3.2. The factorisation of the polyphase matrix} The factorisation of the polyphase matrix is obtained by using Theorem 7 of [2]. We first apply the Euclidean algorithm to divide \inlineTFinmath{{h_e}(z)} by \inlineTFinmath{{h_o}(z)}. That is, we start by setting \inlineTFinmath{{a_0}={h_e}(z)}, \inlineTFinmath{{b_0}={h_o}(z)}, and then iterate the scheme as in (2), \inlineTFinmath {\MathBegin{MathArray}{l} {a_i}={b_{i-1}}, \\ \noalign{\vspace{0.5ex}}{b_i}={a_{i-1}}-{a_i} {q_i},\\ \MathEnd{MathArray}} \noindent choosing quotients \inlineTFinmath{{q_i}} such that the degrees of the remainders \inlineTFinmath{{b_i}} are reduced at each step. The factorisation of the polyphase matrix is given in terms of these quotients by \dispTFNumberedEquationmath{P(z)=\bigg(\prod _{i=1}^{n/2-1}\big(\MathBegin{MathArray}[c]{cc} 1&{q_{2i-1}}(z) \\ 0&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ {q_{2i}}(z)&1 \MathEnd{MathArray}\big)\bigg).\big(\MathBegin{MathArray}[c]{cc} 1&{{\kappa }^2}s(z) \\ 0&1 \MathEnd{MathArray}\big).\bigg(\MathBegin{MathArray}[c]{cc} \kappa &0 \\ 0&\frac{1}{\kappa } \MathEnd{MathArray}\bigg),} \noindent where \(\kappa \) is a constant {\itshape scaling factor} which is simply the gcd, and \inlineTFinmath{s(z) }is a {\itshape lifting term}, found by comparison of the factorised and unfactorised forms of the polyphase matrix. We refer to the matrices in the initial product as {\itshape quotient matrices}, and to the last two matrices as the {\itshape lifting matrix} and the {\itshape scaling matrix} respectively. The alternating upper and lower triangular structure of this factorisation is what allows an 'in-place' implementation of the wavelet transform. Consider the effect of premultiplication of \inlineTFinmath{P(z)} in {\itshape factorised} form by a two-element row vector \inlineTFinmath{({s^{[j]}},{d^{[j]}})} as in (12). The first factor here,\inlineTFinmath{\big(\MathBegin{MathArray}[c]{cc} 1&{q_1}(z) \\ 0&1 \MathEnd{MathArray}\big)}, is upper triangular with diagonal entries equal to 1. The {\itshape lifting step} corresponding to this matrix factor is thus \dispTFNumberedEquationmath{({s^{[j]}},{d^{[j]}}).\big(\MathBegin{MathArray}[c]{cc} 1&{q_1}(z) \\ 0&1 \MathEnd{MathArray}\big) \equiv \big(\MathBegin{MathArray}[c]{c} {s^{[j]}}\leftarrow {s^{[j]}} \\ {d^{[j]}}\leftarrow {d^{[j]}}+{q_1}(z){s^{[j]}} \MathEnd{MathArray}\big) \cdot} Similar steps result from subsequent matrices in the factorisation. This illustrates the fact that the lifting method is in-place and trivially invertible because \inlineTFinmath{{{\big(\MathBegin{MathArray}[c]{cc} 1&x \\ 0&1 \MathEnd{MathArray}\big)}^{-1}}=\big(\MathBegin{MathArray}[c]{cc} 1&-x \\ 0&1 \MathEnd{MathArray}\big)} and \inlineTFinmath{{{\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ x&1 \MathEnd{MathArray}\big)}^{-1}}=\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ -x&1 \MathEnd{MathArray}\big)}. Note that in general the quotients contain powers of {\itshape z}, which will give rise to shifting of the lists. The lifting scheme can thus be written as a sequence of steps of the form (14) with alternating {\bfseries {\itshape s}} and {\bfseries {\itshape d}} symbols, followed by the scaling step. The initial values of \inlineTFinmath{s} and \inlineTFinmath{d} are the results of the lazy wavelet transform as in (11). When we reach the end of the factorisation, that is, when the scaling steps have been executed, we record the final values of \inlineTFinmath{{s^{[j]}}} and \inlineTFinmath{{d^{[j]}}}. As before, \inlineTFinmath{{d^{[j]}}} represents the high frequency component at scale {\itshape j}; this is one component of the wavelet transform. We conduct the lazy wavelet transform on the result of scaling \inlineTFinmath{{s^{[j]}}}, and follow the same lifting steps as before. These ideas are illustrated by the examples in \S{}5. The scaling matrix, \inlineTFinmath{\bigg(\MathBegin{MathArray}[c]{cc} \kappa &0 \\ 0&\frac{1}{\kappa } \MathEnd{MathArray}\bigg)}, is implemented as \inlineTFinmath{{s^{[j]}}\leftarrow \kappa {s^{[j]}}} and \inlineTFinmath{{d^{[j]}}\leftarrow \frac{1}{\kappa } {d^{[j]}}}. For implementation purposes it is preferable to have a factorisation with {\itshape constant} gcd. The gcd is always monomial [2], and, for a general gcd, we can factorise the scaling matrix into lifting steps followed by a constant scaling matrix. To see this, note that \dispTFNumberedEquationmath{\MathBegin{MathArray}{l} \big(\MathBegin{MathArray}[c]{cc} x&0 \\ 0&1/x \MathEnd{MathArray}\big)=\big(\MathBegin{MathArray}[c]{cc} 1&x-{x^2} \\ 0&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ -1/x&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&x-1 \\ 0&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ 1&1 \MathEnd{MathArray}\big)\equiv \\ \noalign{\vspace{1.16667ex}} \hspace{2.em} \big(\MathBegin{MathArray}[c]{cc} 1&0 \\ -1&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&1-1/x \\ 0&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ x&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&1/{x^2}-1/x \\ 0&1 \MathEnd{MathArray}\big),\\ \MathEnd{MathArray}} \noindent where {\itshape x} represents \inlineTFinmath{{z^{\pm 1}}}-- only in this case do the matrices in these factorisations have degree 1 entries. So for an arbitrary gcd, \inlineTFinmath{\alpha {x^m}, m>0}, we could factorise the scaling matrix as \inlineTFinmath {\MathBegin{MathArray}{l} \big(\MathBegin{MathArray}[c]{cc} \alpha {x^m}&0 \\ 0&1/(\alpha {x^m}) \MathEnd{MathArray}\big)={{\big(\MathBegin{MathArray}[c]{cc} x&0 \\ 0&1/x \MathEnd{MathArray}\big)}^m}.\big(\MathBegin{MathArray}[c]{cc} \alpha &0 \\ 0&1/\alpha \MathEnd{MathArray}\big)= \\ \noalign{\vspace{1.25ex}} \hspace{2.em} {{\big(\big(\MathBegin{MathArray}[c]{cc} 1&x-{x^2} \\ 0&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ -1/x&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&x-1 \\ 0&1 \MathEnd{MathArray}\big).\big(\MathBegin{MathArray}[c]{cc} 1&0 \\ 1&1 \MathEnd{MathArray}\big)\big)}^m}\big(\MathBegin{MathArray}[c]{cc} \alpha &0 \\ 0&1/\alpha \MathEnd{MathArray}\big),\\ \MathEnd{MathArray}} \noindent and still have alternating upper and lower triangular matrices with degree 1 entries followed by a constant scaling matrix. Similarly we could factorise the matrix using the second elementary factorisation in (15), or indeed a combination of the two. Hence there are \inlineTFinmath{{2^m}} ways of factorising the scaling matrix using these elementary factorisations, quite apart from the fact that we can choose the quotients in different ways. \Subsection*{3.3. The number of factorisations} It is important to consider the number of possible factorisations of a given wavelet transform. This is equivalent to the question of how many factorisations exist for two Laurent polynomials. In the following analysis, we only consider quotients of degree at most 1. Suppose \inlineTFinmath{{a_0}} and \inlineTFinmath{{b_0}} are our starting polynomials. We first consider the case where \inlineTFinmath{\LeftBracketingBar {a_0}\RightBracketingBar =\LeftBracketingBar {b_0}\RightBracketingBar +1}. Consider the first step of the Euclidean algorithm as in (2), \dispTFNumberedEquationmath{\MathBegin{MathArray}{l} {a_1}={b_0}, \\ \noalign{\vspace{0.5ex}}{b_1}={a_0}-{a_1} {q_1}.\\ \MathEnd{MathArray}} The aim is to choose \inlineTFinmath{{q_1}} so that the degree of {\itshape b} is reduced ({\itshape i.e.}, so \inlineTFinmath{\LeftBracketingBar {b_1}\RightBracketingBar <\LeftBracketingBar {b_0}\RightBracketingBar }). Since we are using degree 1 quotients, we have that \inlineTFinmath{\LeftBracketingBar {b_0}{q_1}\RightBracketingBar =\LeftBracketingBar {b_0}\RightBracketingBar +\LeftBracketingBar {q_1}\RightBracketingBar =\LeftBracketingBar {a_0}\RightBracketingBar }, that is, \inlineTFinmath{{a_0}} and \inlineTFinmath{{b_0}{q_1}} have the {\itshape same }number of terms. We must choose the quotient so that the highest (or equivalently the lowest) powers of \inlineTFinmath{{a_0}} and \inlineTFinmath{{b_0}{q_1}} are the same. If this were not so, there would be extra powers introduced by the subtraction, and we could not reduce the degree. Since a generic degree 1 quotient has the form \inlineTFinmath{{z^n}\big(c+\frac{d}{z}\big)}, this is simply a matter of choosing {\itshape n} appropriately. In order to reduce the degree, we may choose to eliminate the two highest or two lowest powers, or each extreme power, by choice of the quotient coefficients {\itshape c} and {\itshape d}. Clearly there is no point eliminating intermediate terms, since this will have no effect on the degree. Thus we have three possibilities at this step. After this step, we set \inlineTFinmath{{a_1}={b_0}} and proceed similarly to (16) to find \inlineTFinmath{{b_2}}. The degrees of {\itshape a} and {\itshape b} have each been reduced by one, so we proceed as in the first step. Thus we have three possibilities at each step {\itshape except} the last one. This is because at the last step, {\itshape a} has degree 1, and all options are {\itshape equivalent} -- there are only two terms that can be eliminated. Since we reduce the degree of {\itshape b} by 1 at each step, the number of steps in the Euclidean algorithm is \inlineTFinmath{\LeftBracketingBar b\RightBracketingBar +1}, and hence the number of factorisations in this case is \inlineTFinmath{{3^{\LeftBracketingBar b\RightBracketingBar }}}. Now let us consider the case where \inlineTFinmath{\LeftBracketingBar a\RightBracketingBar =\LeftBracketingBar b\RightBracketingBar }. Considering again (16), we see that the degree can be reduced if \inlineTFinmath{{q_1}} has degree 0 ({\itshape i.e.}, is monomial) under certain circumstances. We can cover all cases by defining a generic quotient as before. This time, however, \inlineTFinmath{{b_0}{q_1}} has one more term than \inlineTFinmath{{a_0}}. We can still eliminate each extreme power, but this time we must choose whether we want to match the highest {\itshape or} lowest powers in \inlineTFinmath{{a_0}} and \inlineTFinmath{{b_0}{q_1}} -- these choices are no longer equivalent. We eliminate the highest two powers or lowest two powers as before. The degree 0 quotients are arise when we set c or d to zero when eliminating terms. Hence there are now {\itshape four} possible elimination options at the first step. After the first step, we set \inlineTFinmath{{a_1}={b_0}} as above. This time, however, the degree of {\itshape b} has been reduced by 1, but the degree of {\itshape a} has {\itshape not} changed. Thus we are in the situation where \inlineTFinmath{\LeftBracketingBar a\RightBracketingBar =\LeftBracketingBar b\RightBracketingBar +1}, and we can proceed as above. The number of factorisations in this case is thus \inlineTFinmath{4\times {3^{\LeftBracketingBar b\RightBracketingBar -1}}}. If {\itshape a} and {\itshape b} have degrees differing by more than 1, we {\itshape cannot} execute the Euclidean algorithm using only quotients of degree at most 1. To see this, consider again the first step (16), and suppose that \inlineTFinmath{\LeftBracketingBar {a_0}\RightBracketingBar >\LeftBracketingBar {b_0}\RightBracketingBar +1}. For degree 1 quotients, we can cancel two terms from \inlineTFinmath{{a_0}}, so \inlineTFinmath{\LeftBracketingBar {b_1}\RightBracketingBar =\LeftBracketingBar {a_0}\RightBracketingBar -2>\LeftBracketingBar {b_0}\RightBracketingBar -1}, that is \inlineTFinmath{\LeftBracketingBar {b_1}\RightBracketingBar \geq \LeftBracketingBar {b_0}\RightBracketingBar } (since we are dealing with positive integers) so we have {\itshape not} reduced the degree of {\itshape b}. This is not a great limitation, since we are not dealing with {\itshape arbitrary} Laurent polynomials, but {\itshape even and odd components of polyphase filters}. Usually the degrees of the even and odd components will have degrees differing by 0 or 1. The only exception would be if one of the filter coefficients in the range \inlineTFinmath{{k_b}