Wavelet transform in Haskell
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.
A few words on design choices of the implementation
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.
Haddock documentation of the Wavelets module
Data types
data WaveletFilter t
Wavelet filter type containing the wavelet linear filter, its
corresponding smoothing filter and their common length.
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.
Wavelet types
haar :: Floating t => WaveletFilter t
Haar wavelet, the simplest wavelet, which has length 2.
daubechies_la :: Floating t => Int -> WaveletFilter 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.
daubechies_bl :: Floating t => Int -> WaveletFilter t
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.
coiflet :: Floating t => Int -> WaveletFilter t
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.
Creating custom wavelets
verifywavelet :: (Floating t, Eq t) => WaveletFilter t -> (t, t)
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).
Wavelet coefficient packers
wp_interleave :: Num t => WaveletPacker t [t]
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).
wp_separate :: Num t => WaveletPacker t [[t]]
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.
Transform functions
dwt :: Num t => Int -> WaveletFilter t -> WaveletPacker t c -> [t] -> c
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).
idwt :: Num t => Int -> WaveletFilter t -> WaveletPacker t c -> c -> [t]
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.
idwtsynth :: Num t => Int -> WaveletFilter t -> WaveletPacker t c -> c -> [t]
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