Abstract
The discrete wavelet transform can be decomposed into a finite sequence of simple filtering steps (lifting steps). This decomposition corresponds to a factorization of the polyphase matrix of the wavelet filters into elementary matrices. The decomposition asymptotically reduces the computational complexity of the transform by a factor two. The scheme can be generalized for integer-to-integer and other non-linear transforms.
Filters and Polynomials
The widely-used $z$-transform notation is used for the description of underlying wavelet filters. The $z$-transform of a FIR filter $h$ is a Laurent polynomial $h(z)$ given by $$ h(z) = \sum h_k z^{-k}. $$ Exact division of the polynomials is not possible in general. However, division with remainder is possible. Take two polynomials $a(z)$ and $b(z) \neq 0$, then there always exists a polynomial $q(z)$ (the quotient) and a polynomial $r(z)$ (the remainder) so that $$ a(z) = b(z) q(z) + r(z). $$ If $b(z)$ is a monomial, then $r(z) = 0$ and the division is exact. The long division of Laurent polynomials is not necessarily unique. This fact is particularly useful in the factoring algorithm.
Wavelet Transforms
The figure below shows the general block scheme of a discrete wavelet transform.
Discrete Wavelet Transform
The forward transform uses two analysis FIR filters $\tilde{h}$ and $\tilde{g}$ followed by subsampling. The inverse transform first upsamples and then uses two synthesis filters $h$ and $g$. The polyphase representation of a filter $h$ is given by $$ h(z) = h_e(z^2) + z^{-1} h_o(z^2), $$ where $h_e$ contains the even coefficients, and $h_o$ contains the odd coefficients. The polyphase matrix is assembled as $$ P(z) = \begin{bmatrix} h_e(z) & g_e(z) \\ h_o(z) & g_o(z) \end{bmatrix}. $$ The $\tilde{P}(z)$ is defined similarly.
Polyphase Representation of Discrete Wavelet Transform
The discrete wavelet transform is now represented in the figure above.
Factoring Algorithm
Any pair of complementary filters $(h,g)$ can be factored into lifting steps. The fundamental part is to run the Euclidean algorithm starting from $h_e(z)$ and $h_o(z)$. Thus one have $$ \begin{bmatrix} h_e(z) \\ h_o(z) \end{bmatrix} = \prod_{i=1}^{n} \begin{bmatrix} q_i(z) & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} K \\ 0 \end{bmatrix}, $$ where $K$ is a constant. One can always assume that $n$ is even. Given a filter $h$, one can always find a complementary filter $g^0$ by letting $$ P^0(z) = \begin{bmatrix} h_e(z) & g^0_e(z) \\ h_o(z) & g^0_o(z) \end{bmatrix} = \prod_{i=1}^{n} \begin{bmatrix} q_i(z) & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} K & 0 \\ 0 & 1/K \end{bmatrix}. $$ Now observe that $$ \begin{bmatrix} q_i(z) & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} 1 & q_i(z) \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ q_i(z) & 1 \end{bmatrix}. $$ Using the first equation in case $i$ is odd and the second one in case $i$ is even leads to $$ P^0(z) = \prod_{i=1}^{n/2} \begin{bmatrix} 1 & q_{2i-1}(z) \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ q_{2i}(z) & 1 \end{bmatrix} \begin{bmatrix} K & 0 \\ 0 & 1/K \end{bmatrix}. $$ Since multiple filters $g$ complementary to a filter $h$ exist, the $g$ can be obtained from $g^0$ with one lifting step. Finally, the lifting scheme is given by $$ P(z) = \prod_{i=1}^{m} \begin{bmatrix} 1 & s_i(z) \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ t_i(z) & 1 \end{bmatrix} \begin{bmatrix} K & 0 \\ 0 & 1/K \end{bmatrix}. $$ In other words, every discrete wavelet transform can be obtained by starting with so-called lazy wavelet, followed by $m$ pairs of lifting steps, followed with a scaling. Finally, the dual polyphase matrix is given by $$ \tilde{P}(z) = \prod_{i=1}^{m} \begin{bmatrix} 1 & 0 \\ -s_i(z^{-1}) & 1 \end{bmatrix} \begin{bmatrix} 1 & -t_i(z^{-1}) \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 1/K & 0 \\ 0 & K \end{bmatrix}. $$
References