|↑↑ Home||↑ Science|
Among many other things, Haskell is quite well suited for modelling signal processing, so well it has served as the basis for a circuit description language. All the more surprising that no cabal package for wavelet analysis and transform seems to exist. Wavelet analysis can be expressed as a "pyramid algorithm" made up of identical stages. At each stage, the data stream is processed by a detail filter and a smoothing filter, both reducing the data rate by half and resulting in wavelet and smoothing coefficients, the latter of which are processed further in the next stage. The recursive nature of this transform (and its inverse) makes it a good fit for Haskell.
My implementation mostly follows the book of Percival and Walden on the topic. For example, I have adopted their modifications of the least-asymmetric Daubechies wavelets. In addition I have made some design choices of my own: For wavelets longer than two (all except Haar wavelets), some authors seem to divider the data stream up into fixed-length blocks and make the convolutions necessary for the transform wrap around at their boundaries. I do not do this, as it would defeat the recursiveness of the algorithm and can lead to jumps in the data resulting from the inverse transform (synthesis). A consequence of this decision is that a wavelet transform generates coefficients corresponding to data before the start of the original series, and a following inverse transform will synthesise leading almost-zero values. To handle this, two inverse transform functions exist: idwt cuts off the spurious prepended values, while idwtsynth, which is intended for synthesis, does not. Another design choice is to align wavelets at each stage by their middle coefficients, so the same wavelet coefficients from different scales (stages) correspond to the same position in the original data stream.
As a purely organisational matter, I have separated the transform itself from the way the resulting coefficients are stored. This allows to choose a convenient data structure by choosing or creating a WaveletPacker.
data WaveletFilter t
Wavelet filter type containing the wavelet linear filter, its corresponding smoothing filter and their common length.
|MakeWaveletFilter [t] [t] Int|
|Show t => Show (WaveletFilter t)|
data WaveletPacker t c
A data type for storing wavelet coefficients after performing the transform. The first two parameters are a start value and a function for adding the list of wavelet coefficients from a stage of the transform to the stored set. This function will be applied to the smoothing coefficient list and then to the detail coefficients with rising frequency. The other two type parameters are functions for unpacking the stored coefficients and recreating the lists of coefficients for each stage. The third should return the detail coefficients from highest to lowest frequencies with each successive call. The fourth should return the smoothing coefficients.
|MakeWaveletPacker c ([t] -> c -> c) (c -> ([t], c)) (c -> [t])|
Least-asymmetric Daubechies wavelets in even sizes from 4 to 20. These transforms are a good compromise between good phase behaviour, separation of frequency bands and lack of artifacts. The order of the coefficients has been reversed as recommended by Percival and Walden.
Best localised Daubechies wavelets in even sizes from 4 to 20. These transforms are optimised to be nearly zero-phase. The order of the coefficients has been reversed as recommended by Percival and Walden.
Coifman wavelets in mutiple-of-6 sizes from 6 to 30. The order of the coefficients has been reversed as recommended by Percival and Walden.
Create wavelet object from smoothing filter coefficients.
Verify the orthonormality relation of a wavelet filter. Actually the judgement has still to be made by the calling context, as the relation will not be fulfilled exactly for reasons of numerical errors. This function merely computes the RMS deviation of the normality condition (first return value) and orthogonality condition (second).
This packer interleaves wavelet coefficients in one list so that the n-th level coefficient occurs every 2^n entries starting with the 2^(n-1)th (1-based) (n in [1 .. order of the transform] starting with the highest frequency; the smoothing coefficient takes the remaining place at the end of each block of 2^order coefficients).
This packer stores the wavelet coefficients corresponding to different frequencies / levels of detail in separate lists. The highest-frequency coefficients are in the first list. Note that the number of coefficients describing a given stretch of the original time series differs by a factor of two between neighbouring coefficient lists.
Discrete wavelet transform. The first argument is the number of stages of the transform. The result will be (stages + 1) sets of coefficients, including one set of smoothed values. The second argument is the type of wavelet, the third the way to wrap the resulting coefficients. The last argument is the series of data.
The data are interpreted as part of an infinitely long series, not (as is sometimes done) a set of blocks with implied cyclical boundary conditions. To allow reversing the transform, enough zeros are prepended to the data so that the first coefficients from the last stage corresponds to minimal overlap between the filter and its input data (i.e. one or two values).
Inverse discrete wavelet transform reversing the transform of the dwt function. Arguments and operation are the same as for idwtsynth, only the zero values prepended to the data by dwt are removed. One approximate zero value will remain after the original data for wavelets of lengths 4*n+2, because the result of the inverse transform (before removing the prefix) is always even in length.
Inverse discrete wavelet transform suitable for synthesis. The first three arguments are the same as for dwt: The number of stages of the transform, the wavelet type and the coefficient packer. The last argument is the result of the transform with the same packer. The result contains the zero values that dwt has prepended to create a reversible result (which may in fact differ slightly from zero for numerical reasons). This allows to use this function for synthesis (where the first values will usually be non-zero). If you do not want them, use idwt instead.
TOS / Impressum