Introduction to discrete wavelet transforms

This text summarizes key wavelet facts as a convenience for the hasty reader. See, for example, [Mal99, SN96] or [JlCH01] for excellent detailed introductions to the topic. This text is partially based material from [WBHG22].

The fwt relies on convolution operations with filter pairs.

Fig. 1 Overview of the fwt computation.

Fig. 1 illustrates the process. \(\mathbf{h}_\mathcal{A}\) denotes the analysis low-pass filter. \(\mathbf{h}_\mathcal{D}\) the analysis high pass filter. \(\mathbf{f}_\mathcal{A}\) and \(\mathbf{f}_\mathcal{D}\) the synthesis filer pair. \(\downarrow_2\) denotes downsampling with a factor of two, \(\uparrow_2\) means upsampling. In machine learning terms, the analysis transform relies on stride two convolutions. The synthesis or inverse transform on the right works with stride two transposed convolutions. \(\mathbf{H}_{k}\) and \(\mathbf{F}_{k}\) with \(k \in [\mathcal{A}, \mathcal{D}]\) denote the corresponding convolution operators.

\[\mathbf{x}_s * \mathbf{h}_k = \mathbf{c}_{k, s+1}\]

with \(k \in [\mathcal{A}, \mathcal{D}]\) and \(s \in \mathbb{N}_0\) the set of natural numbers, where \(\mathbf{x}_0\) is equal to the original input signal \(\mathbf{x}\). At higher scales, the fwt uses the low-pass filtered result as input, \(\mathbf{x}_s = \mathbf{c}_{\mathcal{A}, s}\) if \(s > 0\). The dashed arrow indicates that we could continue to expand the fwt tree here. ptwt.conv_transform.wavedec() implements this transformation. ptwt.conv_transform.waverec() provides the inverse functionality visible on the right side of Figure Fig. 1.

The wavelet packet transform (pwt) additionally expands the high-frequency part of the tree. The figure below depicts the idea.

Fig. 2 Scematic drawing of the full wpt in a single dimension. Compared to Fig. 1, the high-pass filtered side of the tree is expanded, too.

Whole expansion is not the only possible way to construct a wavelet packet tree. See [JlCH01] for a discussion of other options. In Fig. 1 and Fig. 2, capital letters denote convolution operators. These may be expressed as Toeplitz matrices [SN96]. The matrix nature of these operators explains the capital boldface notation. Coefficient subscripts record the path that leads to a particular coefficient. ptwt.packets.WaveletPacket() provides this functionality for single dimensional inputs.

The two-dimensional transform

This toolbox provides two dimensional input processing functionality. We construct filter quadruples from the original filter pairs to process two-dimensional inputs. The process uses outer products [VYP18]:

\[\mathbf{h}_{a} = \mathbf{h}_\mathcal{A}\mathbf{h}_\mathcal{A}^T, \mathbf{h}_{h} = \mathbf{h}_\mathcal{A}\mathbf{h}_\mathcal{D}^T, \mathbf{h}_{v} = \mathbf{h}_\mathcal{D}\mathbf{h}_\mathcal{A}^T, \mathbf{h}_{d} = \mathbf{h}_\mathcal{D}\mathbf{h}_\mathcal{D}^T\]

With \(a\) for approximation, \(h\) for horizontal, \(v\) for vertical, and \(d\) for diagonal [LGW+19].

With the four filters we are now able to compute,

\[\mathbf{x}_s *_2 \mathbf{h}_k = \mathbf{c}_{k, s+1}\]

with \(k \in [a, h, v, d]\) and \(s \in \mathbb{N}_0\) the set of natural numbers, where \(\mathbf{x}_0\) is equal to the original input image \(\mathbf{X}\). \(*_2\) indicates two dimensional-convolution. Computations at subsequent scales work exclusively with approximation coefficients \(c_{a, s}\) as inputs. The figure below illustrates the process.

Fig. 3 Two-dimensional wavelet transform computation diagram. \(\mathbf{X}\) and \(\hat{\mathbf{X}}\) denote input image and reconstruction respectively.

ptwt.conv_transform_2.wavedec2() and ptwt.conv_transform_2.waverec2() support forward and backward transforms respectively. Potential further decomposition of all coefficient leads us to wavelet packets.

Fig. 4 Two-dimensional wavelet packet transform computation diagram. Dashed lines indicate potential full expansion of the tree.

Fig. 4 illustrates the computation of a full two-dimensional wavelet packet tree. At higher scales, all resulting coefficients from previous scales serve as inputs. The four filters repeatedly convolved with all outputs to build the full tree. The inverse transforms work analogously. ptwt.packets.WaveletPacket2D() provides this functionality. We refer to the standard literature [JlCH01, SN96] for an extended discussion.

Compared to the FWT, the high-frequency half of the tree is subdivided into more bins, yielding a fine-grained view of the entire spectrum. We always show analysis and synthesis transforms to stress that all wavelet transforms are lossless. Synthesis transforms reconstruct the original input based on the results from the analysis transform.

Common wavelets and their properties

A key property of the wavelet transform is its invertibility. Additionally, we expect an alias-free representation. Standard literature like [SN96] formulates the perfect reconstruction and alias cancellation conditions to satisfy both requirements. For an analysis filter coefficient vector \(\mathbf{h}\) the equations below use the polynomial \(H(z) = \sum_n h(n)z^{-n}\). We construct \(F(z)\) the same way using the synthesis filter coefficients in \(\mathbf{f}\). To guarantee perfect reconstruction the filters must respect

\[H_\mathcal{A}(z)F_\mathcal{A}(z) + H_\mathcal{D}(-z)F_\mathcal{D}(z) = 2z^{-l}.\]

Similarly

\[F_\mathcal{A}(z)H_\mathcal{A}(-z) + F_\mathcal{D}(z)H_\mathcal{D}(-z) = 0\]

guarantees alias cancellation.

Filters that satisfy both equations qualify as wavelets. Lets consider i.e. a Daubechies wavelet and a Symlet:

Fig. 5 Visualization of the Symlet 6 filter coefficients.

Fig. 6 Visualization of the Daubechies 6 filter coefficients.

Fig. 5 and Fig. 6 visualize the Daubechies and Symlet filters of 6th degree. Compared to the Daubechies Wavelet family, their Symlet cousins have more mass at the center. Fig. 5 illustrates this fact. Large deviations occur around the fifth filter in the center, unlike the Daubechies’ six filters in Fig. 6. Consider the sign patterns in Fig. 6. The decomposition highpass (orange) and the reconstruction lowpass (green) filters display an alternating sign pattern. This behavior is a possible solution to the alias cancellation condition. To understand why substitute \(F_\mathcal{A}(z) = H_\mathcal{D}(-z)\) and \(F_\mathcal{D} = -H_\mathcal{A}(-z)\) into the perfect reconstruction condition [SN96]. \(F_\mathcal{A}(z) = H_\mathcal{D}(-z)\) requires an opposing sign at even and equal signs at odd powers of the polynomial.