Discrete wavelet transform

In numerical analysis and functional analysis, a discrete wavelet transform (DWT) is any wavelet transform for which the wavelets are discretely sampled. As with other wavelet transforms, a key advantage it has over Fourier transforms is temporal resolution: it captures both frequency "and" location information (location in time).
Examples.
Haar wavelets.
The first DWT was invented by Hungarian mathematician Alfréd Haar. For an input represented by a list of formula_1 numbers, the Haar wavelet transform may be considered to pair up input values, storing the difference and passing the sum. This process is repeated recursively, pairing up the sums to prove the next scale, which leads to formula_2 differences and a final sum.
Daubechies wavelets.
The most commonly used set of discrete wavelet transforms was formulated by the Belgian mathematician Ingrid Daubechies in 1988. This formulation is based on the use of recurrence relations to generate progressively finer discrete samplings of an implicit mother wavelet function; each resolution is twice that of the previous scale. In her seminal paper, Daubechies derives a family of wavelets, the first of which is the Haar wavelet. Interest in this field has exploded since then, and many variations of Daubechies' original wavelets were developed.
The dual-tree complex wavelet transform (DℂWT).
The dual-tree complex wavelet transform (ℂWT) is a relatively recent enhancement to the discrete wavelet transform (DWT), with important additional properties: It is nearly shift invariant and directionally selective in two and higher dimensions. It achieves this with a redundancy factor of only formula_3, substantially lower than the undecimated DWT. The multidimensional (M-D) dual-tree ℂWT is nonseparable but is based on a computationally efficient, separable filter bank (FB).
Others.
Other forms of discrete wavelet transform include the LeGall-Tabatabai (LGT) 5/3 wavelet developed by Didier Le Gall and Ali J. Tabatabai in 1988 (used in JPEG 2000), the Binomial QMF developed by Ali Naci Akansu in 1990, the set partitioning in hierarchical trees (SPIHT) algorithm developed by Amir Said with William A. Pearlman in 1996, the non- or undecimated wavelet transform (where downsampling is omitted), and the Newland transform (where an orthonormal basis of wavelets is formed from appropriately constructed top-hat filters in frequency space). Wavelet packet transforms are also related to the discrete wavelet transform. Complex wavelet transform is another form.
Properties.
The Haar DWT illustrates the desirable properties of wavelets in general. First, it can be performed in formula_4 operations; second, it captures not only a notion of the frequency content of the input, by examining it at different scales, but also temporal content, i.e. the times at which these frequencies occur. Combined, these two properties make the Fast wavelet transform (FWT) an alternative to the conventional fast Fourier transform (FFT).
Time issues.
Due to the rate-change operators in the filter bank, the discrete WT is not time-invariant but actually very sensitive to the alignment of the signal in time. To address the time-varying problem of wavelet transforms, Mallat and Zhong proposed a new algorithm for wavelet representation of a signal, which is invariant to time shifts. According to this algorithm, which is called a TI-DWT, only the scale parameter is sampled along the dyadic sequence 2^j (j∈Z) and the wavelet transform is calculated for each point in time.
Applications.
The discrete wavelet transform has a huge number of applications in science, engineering, mathematics and computer science. Most notably, it is used for signal coding, to represent a discrete signal in a more redundant form, often as a preconditioning for data compression. Practical applications can also be found in signal processing of accelerations for gait analysis, image processing, in digital communications and many others.
It is shown that discrete wavelet transform (discrete in scale and shift, and continuous in time) is successfully implemented as analog filter bank in biomedical signal processing for design of low-power pacemakers and also in ultra-wideband (UWB) wireless communications.
Example in image processing.
Wavelets are often used to denoise two dimensional signals, such as images. The following example provides three steps to remove unwanted white Gaussian noise from the noisy image shown. Matlab was used to import and filter the image.
The first step is to choose a wavelet type, and a level N of decomposition. In this case biorthogonal 3.5 wavelets were chosen with a level N of 10. Biorthogonal wavelets are commonly used in image processing to detect and filter white Gaussian noise, due to their high contrast of neighboring pixel intensity values. Using these wavelets a wavelet transformation is performed on the two dimensional image.
Following the decomposition of the image file, the next step is to determine threshold values for each level from 1 to N. Birgé-Massart strategy is a fairly common method for selecting these thresholds. Using this process individual thresholds are made for N = 10 levels. Applying these thresholds are the majority of the actual filtering of the signal.
The final step is to reconstruct the image from the modified levels. This is accomplished using an inverse wavelet transform. The resulting image, with white Gaussian noise removed is shown below the original image. When filtering any form of data it is important to quantify the signal-to-noise-ratio of the result. In this case, the SNR of the noisy image in comparison to the original was 30.4958%, and the SNR of the denoised image is 32.5525%. The resulting improvement of the wavelet filtering is a SNR gain of 2.0567%.
It is important to note that choosing other wavelets, levels, and thresholding strategies can result in different types of filtering. In this example, white Gaussian noise was chosen to be removed. Although, with different thresholding, it could just as easily have been amplified.
To illustrate the differences and similarities between the discrete wavelet transform with the discrete Fourier transform, consider the DWT and DFT of the following sequence: (1,0,0,0), a unit impulse.
The DFT has orthogonal basis (DFT matrix):
while the DWT with Haar wavelets for length 4 data has orthogonal basis in the rows of:
Preliminary observations include:
The DWT demonstrates the localization: the (1,1,1,1) term gives the average signal value, the (1,1,–1,–1) places the signal in the left side of the domain, and the 
(1,–1,0,0) places it at the left side of the left side, and truncating at any stage yields a downsampled version of the signal:
The DFT, by contrast, expresses the sequence by the interference of waves of various frequencies – thus truncating the series yields a low-pass filtered version of the series:
Notably, the middle approximation (2-term) differs. From the frequency domain perspective, this is a better approximation, but from the time domain perspective it has drawbacks – it exhibits undershoot – one of the values is negative, though the original series is non-negative everywhere – and ringing, where the right side is non-zero, unlike in the wavelet transform. On the other hand, the Fourier approximation correctly shows a peak, and all points are within formula_10 of their correct value, though all points have error. The wavelet approximation, by contrast, places a peak on the left half, but has no peak at the first point, and while it is exactly correct for half the values (reflecting location), it has an error of formula_11 for the other values.
This illustrates the kinds of trade-offs between these transforms, and how in some respects the DWT provides preferable behavior, particularly for the modeling of transients.
Definition.
One level of the transform.
The DWT of a signal formula_12 is calculated by passing it through a series of filters. First the samples are passed through a low pass filter with impulse response formula_13 resulting in a convolution of the two:
The signal is also decomposed simultaneously using a high-pass filter formula_15. The outputs give the detail coefficients (from the high-pass filter) and approximation coefficients (from the low-pass). It is important that the two filters are related to each other and they are known as a quadrature mirror filter.
However, since half the frequencies of the signal have now been removed, half the samples can be discarded according to Nyquist’s rule. The filter output of the low-pass filter formula_13 in the diagram above is then subsampled by 2 and further processed by passing it again through a new low- pass filter formula_13 and a high- pass filter formula_15 with half the cut-off frequency of the previous one,i.e.:
This decomposition has halved the time resolution since only half of each filter output characterises the signal. However, each output has half the frequency band of the input, so the frequency resolution has been doubled.
With the subsampling operator formula_21
the above summation can be written more concisely.
However computing a complete convolution formula_25 with subsequent downsampling would waste computation time.
The Lifting scheme is an optimization where these two computations are interleaved.
Cascading and filter banks.
This decomposition is repeated to further increase the frequency resolution and the approximation coefficients decomposed with high and low pass filters and then down-sampled. This is represented as a binary tree with nodes representing a sub-space with a different time-frequency localisation. The tree is known as a filter bank.
At each level in the above diagram the signal is decomposed into low and high frequencies. Due to the decomposition process the input signal must be a multiple of formula_1 where formula_27 is the number of levels.
For example a signal with 32 samples, frequency range 0 to formula_28 and 3 levels of decomposition, 4 output scales are produced:
Relationship to the mother wavelet.
The filterbank implementation of wavelets can be interpreted as computing the wavelet coefficients of a discrete set of child wavelets for a given mother wavelet formula_29. In the case of the discrete wavelet transform, the mother wavelet is shifted and scaled by powers of two
formula_30
where formula_31 is the scale parameter and formula_32 is the shift parameter, both which are integers.
Recall that the wavelet coefficient formula_33 of a signal formula_34 is the projection of formula_34 onto a wavelet, and let formula_34 be a signal of length formula_37. In the case of a child wavelet in the discrete family above,
formula_38
Now fix formula_31 at a particular scale, so that formula_40 is a function of formula_32 only. In light of the above equation, formula_42 can be viewed as a convolution of formula_34 with a dilated, reflected, and normalized version of the mother wavelet, formula_44, sampled at the points formula_45. But this is precisely what the detail coefficients give at level formula_31 of the discrete wavelet transform. Therefore, for an appropriate choice of formula_47 and formula_48, the detail coefficients of the filter bank correspond exactly to a wavelet coefficient of a discrete set of child wavelets for a given mother wavelet formula_29.
As an example, consider the discrete Haar wavelet, whose mother wavelet is formula_50. Then the dilated, reflected, and normalized version of this wavelet is formula_51, which is, indeed, the highpass decomposition filter for the discrete Haar wavelet transform.
Time complexity.
The filterbank implementation of the Discrete Wavelet Transform takes only O("N") in certain cases, as compared to O("N" log "N") for the fast Fourier transform.
Note that if formula_48 and formula_47 are both a constant length (i.e. their length is independent of N), then formula_54 and formula_55 each take O("N") time. The wavelet filterbank does each of these two O("N") convolutions, then splits the signal into two branches of size N/2. But it only recursively splits the upper branch convolved with formula_48 (as contrasted with the FFT, which recursively splits both the upper branch and the lower branch). This leads to the following recurrence relation
which leads to an O("N") time for the entire operation, as can be shown by a geometric series expansion of the above relation.
As an example, the discrete Haar wavelet transform is linear, since in that case formula_47 and formula_48 are constant length 2.
The locality of wavelets, coupled with the O("N") complexity, guarantees that the transform can be computed online (on a streaming basis). This property is in sharp contrast to FFT, which requires access to the entire signal at once. It also applies to the multi-scale transform and also to the multi-dimensional transforms (e.g., 2-D DWT).
Code example.
In its simplest form, the DWT is remarkably easy to compute.
The Haar wavelet in Java:
public static int[] discreteHaarWaveletTransform(int[] input) {
 // This function assumes that input.length=2^n, n>1
 int[] output = new int[input.length];
 for (int length = input.length / 2; ; length = length / 2) {
 // length is the current length of the working area of the output array.
 // length starts at half of the array size and every iteration is halved until it is 1.
 for (int i = 0; i < length; ++i) {
 int sum = input[i * 2] + input[i * 2 + 1];
 int difference = input[i * 2] - input[i * 2 + 1];
 output[i] = sum;
 output[length + i] = difference;
 if (length == 1) {
 return output;
 //Swap arrays to do next iteration
 System.arraycopy(output, 0, input, 0, length);
Complete Java code for a 1-D and 2-D DWT using Haar, Daubechies, Coiflet, and Legendre wavelets is available from the open source project: JWave.
Furthermore, a fast lifting implementation of the discrete biorthogonal CDF 9/7 wavelet transform in C, used in the JPEG 2000 image compression standard can be found here (archived 5 March 2012).
Example of above code.
This figure shows an example of applying the above code to compute the Haar wavelet coefficients on a sound waveform. This example highlights two key properties of the wavelet transform: