Documentation

wt1d is the 1-dimensional discrete wavelet transform.

Usage: wt1d x nd jmin jmax

• x - input vector of size \( M = 2^{jmax} \).
• nd - Daubechies wavelet identifier \( N \), where \( 2 \leq N \leq 10 \).
• jmin - Last coarse vector has size \( 2^{jmin} \).
• jmax - input vector size is \( 2^{jmax} \).

The Daubechies wavelets are defined in TenLectures. The input vector is the initial detail vector. The first level of the transform splits this vector into two parts, coarse (coarse coefficients---see below) and detail (detail coefficients). Then the coarse vector is split into two parts, one coarse and the other detail (the previous detail vector is unchanged). This continues until the coarse vector is reduced to size \( 2^{jmin} \), where \( jmax > jmin \ge 0 \). The output vector returned by this function has the same size as the input vector, and it has the form \( (coarse, detail, d, d', d'',...) \), where coarse is the last coarse vector, detail is its companion, and d, d', d'', etc. are detail vectors from previous steps, increasing in size by factors of 2. The inverse wavelet transform iwt1d starts with this vector and works backwards to recover the original input vector. Of course, it needs jmin to know where to start. See Examples below.

Here are the fundamental relations derived in Compact1988 that define the discrete compactly supported wavelet transform and its inverse \[ a_k^{j+1} = \sum_n h(n-2k) a_n^j,\qquad d_k^{j+1} = \sum_n g(n-2k) a_n^j, \] where \( a_k^j \) and \( d_k^j \) are the approximation and detail coefficients, respectively. See Compact1988, pp.935-938. The reverse or inverse DWT is defined by \[ a_n^{j-1} = \sum_k h(n-2k) a_k^j + \sum_k g(n-2k) d_k^j. \] It is clear from this that decimation by two and convolution is involved at each step. More precisely, these equations and the corresponding Haskell implementation makes it clear that in the forward transform the filters are reversed, convolved with the input (prior level coefficients), and subsampled, while in the reverse transform, the prior level coefficients are upsampled, convolved with the filters (not reversed), and summed. Implementations in Python, R, Julia and Matlab can be found in Data2021 and NumericalTours.

Note that much of the analysis in Compact1988 is done where n varies over all integers. In the implementation below we work modulo the size of the input vector, which turns ordinary convolutions into circular convolutions in the usual way.

The \( N \)-th Daubechies wavelet filter h has support \( [0,2N-1] \) (even number of points). See Compact1988, Table 1, p.980, and TenLectures, Table 6.1, p.195 (higher precision). The novelty of this work was the discovery of invertible filters with compact support (the Shannon wavelet does not have compact support). This requires detailed Fourier analysis, from which it is deduced that the conjugate filter g can be chosen to satisfy \( g(k) = (-1)^k h(2p+1-k) \) for a convenient choice for p (see TenLectures, p.136). To define g with the same support as h, for the N-th wavelet, use p = N-1, so \( g(k) = (-1)^k h(2N-1-k), k=0,1,...,2N-1 \).

Following NumericalTours, we define the circular convolution cconv1d in such a way that the filter h is centered, and for this purpose a zero is prepended to both h and g (so these vectors have an odd number of points, with the understanding that their support does not change).

Examples:

>>> dist :: [Double] -> [Double] -> Double
>>> dist x y = sqrt (sum (zipWith (\x y -> (x-y)^2) x y))
>>> sig = deltaFunc 5 1024
>>> forward = wt1d sig 10 5 10 -- wavelet 10, jmin=5, jmax=10
>>> backward = iwt1d forward 10 5 10
>>> dist sig backward

1.2345e-12

>>> sig = deltaFunc 5 1024
>>> wavelet2 = iwt1d sig 2 0 10
>>> [rgraph| plot(wavelet2_hs,type='l') |]

>>> x = deltaFunc 100 1024
>>> y = deltaFunc 400 1024
>>> z = zipWith (\x y -> 2*x + y) x y
>>> wt = wt1d z 2 0 10
>>> [rgraph| plot(wt_hs,type='l') |]

iwt1d is the 1-dimensional inverse discrete wavelet transform.

Usage: iwt1d x nd jmin jmax

• x - The output from wt1d when the last coarse vector has size \( 2^{jmin} \), where \( 0 \leq jmin \lt jmax \).
• nd - Daubechies wavelet identifier \( N \), where \( 2 \leq N \leq 10 \).
• jmin - Last coarse vector has size \( 2^{jmin} \).
• jmax - Output vector has size is \( 2^{jmax} \).

Reverses the steps taken in wt1d above.

cconv1d is the 1-dimensional circular convolution (with centering).

Usage: cconv1d h x        

• x - input signal
• h - filter to convolve with.

Let N = length x, and M = length h, and assume \( M \leq N \). The convolution with centering offset p is defined by \[ y_i = \sum_{j=0}^{M-1} h_j x_{i-j+p},\qquad i=0,1,2,...,N-1, \] where \( p = (M-1)/2 \) when M is odd. The idea is to have each \( y_i \) equal a weighted average of values of \( x_j \) for \( j \) near \( i \), because \( i \) is in the middle of the support of \( h \), approximately.

For this to be well-defined for the specified range of \( i \), \( x_i \) must be extended periodically outside of the range \( 0 \leq i \leq N-1 \), like this... \[ x_{J} \cdots x_{N-1} x_0 x_1 \cdots x_{N-1} x_0 x_1 \cdots x_{p-1}, \] where \( J \) is determined by the condition that we need \( M-1-p \) extra elements on the left, so we have \( (N-1) - J + 1 = (M-1) - p \), and \( J = N-M+1+p \). Since subscripts begin with 0, we must drop the first \( J \) elements of \( x \) to get the padding on the left.

Consequently, padding on the left is determined by dropping \( N-M+1+p \) elements of \( x \), and padding on the right is determined by taking the first p elements of \( x \). The circular convolution is then found by breaking the extended list into sublists of length \( M \), and taking the inner produce of each of these sublists with \( h \) reversed (see source code).

Examples:

>>> cconv1d [1..4] [1..4]

[28,26,20,26]

cconv1dnc standard circular convolution without centering. Same as cconv1d with \( p = 0 \). No padding on the right is needed in this case.

Examples:

>>> cconv1dnc [1..4] [1..4]

[26,28,26,20]

R equivalent: convolve(1:4,1:4,type='circular',conj=FALSE)

conv1d is the 1-dimensional conventional convolution (no FFT).

Usage: conv1d h x        

• x - input signal
• h - filter to convolve with.

Examples:

>>> conv1d [1..4] [1..4]

[1,4,10,20,25,24,16]

R equivalent: convolve(1:4,rev(1:4), type='open')
Python: np.convolve([1,2,3,4], [1,2,3,4]
Octave/Matlab: conv(1:4,1:4)

deltaFunc creates a list of zero's except for one element that equals 1

Usage: deltaFunc k n

• n - length of the list
• k - position that equals 1 (zero-based, k < n).