How to compute the wavelet approximation of a function?

Question

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.

Answer 1

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$

Answer 2

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[0]) + len(fourier_coef[1]),len(chebyshev_coef), 1 + len(haar_coef[1]))
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)