The Daubechies4 wavelet and the forward lifting transform (23) 

Obtain the Daubechies 4 (D4, also called db2) wavelet coefficients:

In[22]:=

Out[22]=

Create two corresponding Laurent polynomials:

In[23]:=

Out[23]=

Compute the Laurent GCD and list of quotients:

In[24]:=

Out[24]=

Create a list of 2×2 matrices that represent the operation of reducing a pair of Laurent polynomials by their quotient:

In[25]:=

Out[26]=

Take their product:

In[27]:=

Out[27]=

Applying this product to the vector comprised of the GCD and 0 recovers the original pair of Laurent polynomials:

In[28]:=

Out[28]=

Invert the reducing matrices:

In[29]:=

Out[30]=

Form the reverse-product of these inverses:

In[31]:=

Out[31]=

Check that this is in fact the inverse:

In[32]:=

Out[32]=

This inverse product matrix provides the Bezout multipliers (top row) and syzygy (bottom row), as can be seen by its action on the input pair:

In[33]:=

Out[33]=

Create the quadrature mirror filter from the D4 coefficients:

In[34]:=

Out[34]=

A certain upper triangular matrix is needed to recover these Laurent polynomials:

In[35]:=

Out[35]=

Form the augmented set of 2×2 matrices required to create the wavelet polyphase matrix:

In[36]:=

Out[37]=

Compute the polyphase matrix:

In[38]:=

Out[38]=

Show that the second column is indeed the quadrature mirror filter:

In[39]:=

Out[39]=

The same matrix can be obtained from WaveletFilterCoefficients:

In[40]:=

Out[41]=

Create a function to recover the discrete wavelet transform (DWT) from these 2×2 matrix factors. For this, powers of the Laurent variable z become shifts with positive exponents being leftward, and indices are modulo the length of the sequence (that is, using periodic wrapping):

In[42]:=

Work with {1,2,…,16} split into odds and evens (the so-called "Lazy Wavelet Transform"):

In[43]:=

Compute the DWT:

In[44]:=

Out[44]=

Compare to the forward lift from the LiftingFilterData for the D2 wavelet:

In[45]:=

Out[45]=

Get the same result directly from LiftingWaveletTransform:

In[46]:=

Out[46]=

Invert the process to recover the original values, using Riffle at the end to undo the Lazy Wavelet Transform:

In[47]:=

Out[48]=

This is equivalent to using InverseWaveletTransform on the DiscreteWaveletData object:

In[49]:=

Out[49]=

The conventional way to use the quotient matrices and their inverses to form a lifting filter actually makes them alternating upper and lower triangular. For simplicity of exposition, this step was omitted.

Compressing an image with the Daubechies10 wavelet (12) 

Create the even and odd Laurent polynomials that arise from the Daubechies10 wavelet coefficients:

In[50]:=

Out[52]=

Compute an extended GCD with quotients coefficients that are reasonably close in magnitude to unity:

In[53]:=

Out[53]=

Another good set of quotients for the D10 Laurent pair is given in the reference by Maslen and Abbott; the latter happens to be the one provided by WaveletFilterCoefficients:

In[54]:=

Out[55]=

Create the quotient matrices and the scaling value (a factor of z² is required for this particular factorization):

In[56]:=

Out[64]=

Set up functions that transform from Laurent polynomials to vector operations:

In[65]:=

Pad a square image to have a power of 2 dimension, removing the average pixel value from each row:

In[66]:=

Out[72]=

Compute the DWT on rows and then compute the DWT on the columns of that result:

In[73]:=

Create a sparse array, chopping values larger than 0.1 to reduce size:

In[74]:=

Out[75]=

Create the inverse operations:

In[76]:=

Perform one inversion and compare to the first level of the previous discrete transform:

In[77]:=

Out[78]=

Perform a second inversion and compare to the input:

In[79]:=

Out[81]=

Recover a reasonable approximation to the original image:

In[82]:=

Out[82]=