# How to compute the wavelet approximation of a function?

For the function $f(x)=x$, how to compute the wavelet approximation using Haar basis?

I'm new to wavelet, I'm looking for a package which will do something like this

from mpmath import *
mp.dps = 15; mp.pretty = True
nb_coeff = 3
interval = [0, 10]
fct = lambda x : x
fourier_coef = fourier(fct, interval, nb_coeff)
fourier_series_apx = lambda x: fourierval(fourier_coef, interval, x)


Can this be done with pywavelets on an Haar basis?
Or, is there another library which will produce the functionality above for Haar wavelet?

Edit
Suspecting that their maybe no package which provides the functionality above, I tried to roll my own.
My implementation currently doesn't work and I asked a question here.

• How do you want to represent $f$ and its approximation in your code? As far as I can see it, pywavelets uses discrete representations.
– Jan
Oct 23, 2013 at 11:45
• @Jan I want to represent $f$ trough its truncated Haar series. I do this to learn about wavelet and see how well they approximate different function. I though maybe pywavelets had this functionality, but you are probably right that it only focus on discrete transform.
– user729
Oct 23, 2013 at 12:02
• I was more on, whether $f$ is given as a discrete vector of values or, e.g., as a symbolic expression. If it is a vector, you can use pywavelets (forward DWT, truncation, backward DWT). If it is an expression, you may want to sample it and use pywavelets. This will give you an idea of how well the approximation will be. You can also implement the expansion into the (truncated) Haar basis by yourself using numerical or symbolic integration. Though, there should by modules around that have this functionality.
– Jan
Oct 23, 2013 at 12:21
• I have had a look at your implementation. Conceptually you seem right. However, the case $j=1$ and, maybe all others, is probably wrong by a scaling factor, unless $N=10$, i.e. your length of the interval. You should check that the generated basis functions are indeed orthonormal. How do you compute the integrals? Your plot looks to me like aliasing. I have once (hard) coded the first $64$ Haar wavelets on the interval $(0,te)$. It is Matlab but maybe it helps you a bit.
– Jan
Oct 24, 2013 at 9:19
– Jan
Oct 25, 2013 at 14:13

This is pure haar scaling function approximation of $f(t)$ $$f(t)=\sum_{k=-\infty}^{\infty}a_k\phi(2^jt-k)$$ where $$\phi(t) = \begin{cases}1 \quad; 0 \leq t < 1,\\0\quad;\mbox{otherwise.}\end{cases}$$ is haar scaling function, $j$ scaling aparameter, $k$ translation parameter and the coefficients $$a_k=\int f(t)\phi(2^jt-k)dt$$ So on our case if $f(t)=t$, for $k\in \mathbb{Z}$ and for $j=0$ we have to compute $a_k=\int f(t)\phi(t-k)dt$. We know $\phi (t-k)=1$ for $t\in [k,k+1)$ and zero elsewhere. So we only integrate over $[k,k+1)$. $$a_k=\int_{k}^{k+1} t \phi(t-k)dt=\frac{2k+1}{2}$$ Hence we have $$f(t)=\frac12\sum_{k=-\infty}^{\infty}(2k+1)\phi(t-k)$$ The following script plots haar approximation versus given function $f(x)=x$

import numpy as np
import matplotlib.pyplot as plt

def init(x,j,k):
y=np.piecewise(x,[x < 0, x > 1],[lambda x: 0, lambda x: 0, lambda x:1])
t=(x+k)/2**j
ii=np.where(y==1)
return y[ii], t[ii]

a=-4 #lower limit of sum
b=4 #upper limit of sum
j=0
k=np.arange(a*2**j,b*2**j)
y=(2*k+1)/2 # the coefficients a_k
T=np.arange(0,2,0.01)
[p,dummy]=init(T*2**j,j,0)
phi=np.ones((len(k),len(p)))
t=np.zeros((len(k),len(p)))

for i in range (1,len(k)+1):
[p,t1]=init(T*2**j,j,i-1)
t[i-1,:]=t1+a

for l in range (1,len(k)+1):
plt.plot(t[l-1,:],2**j*y[l-1]*phi[l-1,:],'r') Haar approximation of $f(x)=cos(x)$ for $j=2$ • I see there must be some mistake in my script. I will try to find it. Thank you.
– user729
Oct 27, 2013 at 14:14

Thanks to Jan. I made my implementation work.
The code below compare: Haar vs Fourier vs Chebyshev. from __future__ import division
from mpmath import *

# --------- Haar wavelet approximation of a function
# algorithm from : http://fourier.eng.hmc.edu/e161/lectures/wavelets/node5.html
# implementation only handle [0,1] for the moment: scaling and wavelet fcts need to be periodice

phi = lambda x : (0 <= x < 1) #scaling fct
psi = lambda x : (0 <= x < .5) - (.5 <= x < 1) #wavelet fct
phi_j_k = lambda x, j, k : 2**(j/2) * phi(2**j * x - k)
psi_j_k = lambda x, j, k : 2**(j/2) * psi(2**j * x - k)

def haar(f, interval, level):
c0 = quadgl(  lambda t : f(t) * phi_j_k(t, 0, 0), interval  )

coef = []
for j in xrange(0, level):
for k in xrange(0, 2**j):
djk = quadgl(  lambda t: f(t) * psi_j_k(t, j, k), interval  )
coef.append( (j, k, djk) )

return c0, coef

def haarval(haar_coef, x):
c0, coef = haar_coef
s = c0 * phi_j_k(x, 0, 0)
for j, k ,djk in coef:
s += djk * psi_j_k(x, j, k)
return s

# --------- to plot an Haar wave
interval = [0, 1]
plot([lambda x : phi_j_k(x,1,1)],interval)

# ---------- main
# below is code to compate : Haar vs Fourier vs Chebyshev

nb_coeff = 5
interval = [0, 1] # haar only handle [0,1] for the moment: scaling and wavelet fcts need to be periodice

fct = lambda x : x

haar_coef = haar(fct, interval, nb_coeff)
haar_series_apx = lambda x : haarval(haar_coef, x)

fourier_coef = fourier(fct, interval, nb_coeff)
fourier_series_apx = lambda x: fourierval(fourier_coef, interval, x)
chebyshev_coef = chebyfit(fct, interval, nb_coeff)
chebyshev_series_apx = lambda x : polyval(chebyshev_coef, x)

print 'fourier %d chebyshev %d haar %d' % ( len(fourier_coef) + len(fourier_coef),len(chebyshev_coef), 1 + len(haar_coef))
print 'error:'
print 'fourier', quadgl(  lambda x : abs( fct(x) - fourier_series_apx(x) ), interval  )
print 'chebyshev', quadgl(  lambda x : abs( fct(x) - chebyshev_series_apx(x) ), interval  )
print 'haar', quadgl(  lambda x : abs( fct(x) - haar_series_apx(x) ), interval  )

plot([fct, fourier_series_apx, chebyshev_series_apx, haar_series_apx], interval)

• Looks good. Although, I wonder why the wavelet apprx. is not symmetric. Have you included the $0$th wavelet, that is constant? Or do you have (#number of basis functions) different from $2^N$?
– Jan
Oct 25, 2013 at 14:10
• @Jan The size of my basis was less than $2**N$. I corrected the code and changed the image. You really have good eyes. Thank you.
– user729
Oct 25, 2013 at 14:35