Computing The Wavelet Transform

Subchapters

Computing The Wavelet Transform

Syntax	`wavelet(`m`,` x`)`
See Also	iwavelet , fspec

Description
Returns the forward wavelet transform of x where x is a real, double-precision, or complex column vector with length equal to a power of 2 (i.e., 2, 4, 8, ...). The integer scalar m specifies the number of wavelet coefficients. If m is 2, the Haar coefficients are used. If m is 4 or 6, the Daubechies coefficients are used. These are the only values of m that are currently implemented.

The vector c of length m is used below to denote the wavelet coefficients. The smoothing operator S applied to a vector y is a vector S(y) with half the number of elements as y. The i-th element of S(y) is defined by



     S (y)  = c  y     + c  y   + . . . + c  y

      i        1  2i-1    2  2i            m  2i+m-2

If an index of y in the summation above exceeds the length of the vector y, the vector is extended circularly; that is, if n is the length of y and j is greater than n, y(j) is defined to be y(j-n). The detail operator D applied to a vector y is a vector D(y) with half the number of elements as y. The i-th element of D(y) is defined by



     D (y)  = c  y     - c    y   + . . . - c  y

      i        m  2i-1    m-1  2i            1  2i+m-2

where y is extended circularly when its index exceeds its length.

The wavelet coefficients are chosen so that the i-th element of D(y) is 0 if the m elements between y(2i-1) and y(2i+m-2) are constant. If m > 4, the i-th element of D(y) is 0 if the m elements satisfy an affine function of i (an affine function has a linear and constant term). If m > 6, the i-th element of D(y) is 0 if the m elements satisfy a quadratic function of i.

We use n to denote the length of x and we use r to denote the return value. The last n/2 elements of r contain D(x); that is,



     D(x) = {r(n/2+1), r(n/2+2), . . . , r(n)}

The previous n/4 elements of r contain D(S(x)); that is,



     D(S(x)) = {r(n/4+1), r(n/4+2), . . . , r(n/2)}

This pattern continues with the details of S(S(x)) stored starting at index n/8+1 of the return value. This procedure is repeated until the last smoothing value is in r(1) and the last detail is in r(2).

Example
We can compute the wavelet transform of a linear vector x by entering




     x = real(seq(16))

     m = 4

     r = wavelet(m, x)

We can then print the high-order details by entering




     format real "f6.3"

     print r.row(9, 8)

to which O-Matrix will respond

Note that all of the details are 0 except for the last one. This is because for 1 < i < n/2, the sequence



     x(2i-1), x(2i), x(2i+1), x(2i+2)

is affine (n is 16 for our example). But when i = n/2, the sequence is



     x(n-1), x(n), x(1), x(2)

which is not linear (due to the circular manner in which x is extended).

For i < n/2, the smoothing values are given by



     S (x) = c  x    + c  x   + c  x     + c  x

      i       1  2i-1   2  2i    3  2i+1    4  2i+2

It follows that these smoothing values are an affine function of the index i for i < n/2. The next order details, D(S(x)), are given by



     D (S(x))  = c  S (x)  - c  S (x)  +  c  S (x)  +  c  S (x)

      i           4  2i-1     3  2i        2  2i+1      4  2i+2

Thus when 2i+2 < n/2 these details are 0; that is the first two elements of D(S(x)) are 0. We can print D(S(x)) by entering




     print r.row(5, 4)

to which O-Matrix will respond