# Connecting DWT matrix operations to Daubechies D4 scaling $\left(\phi(x)\right)$ and wavelet $\left(\psi(x)\right)$ functions by plotting

I'd like to duplicate the plot in the Wikipedia entry for Daubechies D4 scaling and wavelet functions:

I don't have Matlab or Mathematica, and my question is not about hacking the curve, but rather understanding the recursion that generates it. So I wonder if it'd be OK to ask for some pseudo-code implementable in Python or R, basically following the steps in the cascade algorithm, or some breakdown of the steps involved. Ideally I'd like to be able to write some ad hoc code (possibly in R, Python or Octave) to get an idea of the complexity of the process, and the difference between the differences between the discrete transform and the continuous plots above.

I see that there are initial conditions, which are resolved as an eigenvalue problem, and coefficients already worked-up for D4; yet, there is no closed formula, and an iterative numerical process is carried out.

This may involve:

1. Finding function $$H_0,$$ such that $$\hat\phi(2\gamma)=H_0(\gamma)\hat \phi(\gamma),$$ with $$H_0\in L^2(0,1),$$ and $$\hat \phi(\cdot)$$ standing for the Fourier transform.
2. Finding coefficients $$c_k,$$ such that

\begin{align} \frac{1}{\sqrt 2} D^{-1}\phi &=\sum_{k\in\mathbb Z} c_k T_{-k}\phi\\[2ex] \frac{1}{\sqrt 2}\frac{1}{\sqrt 2}\phi\left(\frac x 2\right) &= \sum_{k\in \mathbb Z} c_k\phi(x+k) \end{align}

with the dilation operator $$D_{k} f=\frac{1}{\sqrt k}f\left(\frac x k\right)$$; and the translation operator $$T_{k}f=f(x -k).$$

1. Define $$H_1(\gamma):= H_0\left(\gamma +\frac 1 2\right) e^{-w\pi i\gamma}$$
2. Define the wavelet function $$\psi$$ by its Fourier transform

$$\hat \psi(2\gamma) = H_1(\gamma)\hat\phi(\gamma)$$

All of it rather complicated, and especially puzzling in its assumption of $$\phi,$$ apropos of which Strang mentions:

We never discover the exact value $$\phi(\sqrt 2)$$. It is amazing to compute with a function we do not know — but the applications only require the c's.

After the comment by @Maxtron (thank you!) and using this post on SE, I was able to replicate the plots above of D4 scaling and wavelet functions - well... with some change of sign of the wavelet function:

import numpy as np
import pywt
[phi, psi, x] = pywt.Wavelet('db2').wavefun(level=12)

import pylab
pylab.plot(x, -psi,'-r', label='Wavelet fx, ' r'$$\psi$$')
pylab.plot(x, phi, '-b', label='Scaling fx, ' r'$$\varphi$$')
pylab.axhline(y=0, color='k')
pylab.legend(loc='lower right')
pylab.title('Daubechies D4')
pylab.show()


• There's Magma code on the Wikimedia Commons page of that plot. – Olli Niemitalo Nov 5 '18 at 17:36
• Have you considered pywavelets in Python? link – Maxtron Nov 5 '18 at 22:55
• @Maxtron With some sign swap, pywavelets did the trick; I am really after understanding the process of iteration, though... – Antoni Parellada Nov 6 '18 at 1:09
• @AntoniParellada I can redirect you to a MATLAB implementation of wavelets which should give you good intuition and motivation for using wavelets. Check out this link – Maxtron Nov 6 '18 at 1:57
• @Maxtron Thanks for the link to the Jupyter NB document - I wonder if you were meaning to send me to the Haar wavelet page, as opposed to the Daubechie's wavelet. – Antoni Parellada Nov 6 '18 at 2:43