Function Approximation by Wavelet method If $ f(x)\in{L^{2}(R)}$ , then $ f(x) $ can be projected into the subspace $ V_{j}$ as,
$$
P_{j}f(x) = \sum\limits_{k\in{Z}}c_{j,k}\phi_{j,k}(x)
$$
where $ k\in{Z}$. The projection equation is reasonable because the domian is in real axis. However, in numerical simulation application, only limited samples of function $ f(x) $ ($ f(x_{1}),f(x_{2}),\cdots,f(x_{n})$) are available, my question is how to get the scaling funciton coefficients $c_{j,k}$ numerically based on the given sampled values of $ f(x) $ ? How many number of $ k$ are used during the projection? 
I do some preparation for these problems. I read the Kevin Amaratunga's paper titled "wavelet-Galerkin solutions for one-dimension PDE" and his PHD dissteration, as well as his MATLAB code for approximation of a functions by the scaling function and its translates. I found that the theory in references are different from MATLAB code he shared. 
By the way, I did not find a reference that states the problem (in numerical simulation sense) clear. 
If you are intersted in the problem, you can also check the MATLAB code of Kevin Amaratunga by the link
http://m2matlabdb.ma.tum.de/search_result.jspSearch=13&SearchCategory=AuthorWithFiles&AuthorName=Kevin+Amaratunga
 A: Ok so one thing you have to realize about wavelet analysis is that it exists at the intersection of a wide variety of fields. So, in general, the practice of wavelet analysis will look very different from the theory of wavelets. I'm going to assume that you have a solid understanding of the theory but if any of the steps are unclear, just let me know and I will fill in the details. Anyways, one important relationship that bridges this gap is the dilation equation given by 
$$
\varphi(t)=\sum_{n\in \mathbb{Z}}h_{n}\sqrt{2} \varphi(2t-n)
$$
where $h_{n}$ is the scaling function coefficient but also note that we can write the wavelet function $\psi$ in terms of the scaling function via the relationship
$$
\psi(t)=\sum_{n \in \mathbb{Z}} h^{*}_{k}\sqrt{2}\varphi(2t-n)
$$
where $h^{*}_{n}$ is the wavelet function coefficient which is related to the scaling function coefficient via the relationship
$$
h^{*}_{n}=(-1)^{n}h_{1-n}
$$ 
which then allows us to generate our usual wavelet mother function 
$$\psi_{j,k}(t)=2^{j/2}\psi(2^{j}t-k)$$
Now to see how this works in practice, go back to the original dilation equation 
$$
\varphi(t)=\sum_{n \in \mathbb{Z}} h_{n}\sqrt{2}\varphi(2t-n)
$$
and then dilate and translate $t$ so that
$$
\varphi(2^{j}t-k)=\sum_{n \in \mathbb{Z}} h_{n} \sqrt{2} \varphi(2^{j+1}t-2k-n)
$$
then let $m=2k+n$ and re-write as 
$$
\phi(2^{j}t-k)=\sum_{m \in \mathbb{Z}}h_{m-2k}\sqrt{2}\varphi(2^{j+1}t-m)
$$
Then if we let $\mathcal{V}_{j}=\text{Span}_{k}\{2^{j/2}\varphi(2^{j}t-k)\}$ then
$$
f(t) \in \mathcal{V}_{j+1} \implies f(t)=\sum_{k \in \mathbb{Z}} c_{j+1,k} 2^{(j+1)/2} \varphi(2^{j+1}t-k)$$
but if you go one scale lower in resolution, you need to add the wavelet function for the details lost. So at resolution $j$, we have
$$
f(t)=\sum_{k \in \mathbb{Z}} c_{j,k}2^{j/2}\varphi(2^{j}t-k) + \sum_{k \in \mathbb{Z}} d_{j,k}2^{j/2} \psi(2^{j}t-k)
$$
and we have that 
$$
c_{j,k}=\langle f, \varphi_{j,k} \rangle = \int_{-\infty}^{\infty} f(t)2^{j/2}\psi(2^{j}t-k) dt
$$
and then using our previous relationships, we get 
$$
c_{j,k}=\sum_{m \in \mathbb{Z}}h_{m-2k} \int_{-\infty}^{\infty} f(t) \, 2^{(j+1)/2} \varphi(2^{j+1}t-m) \, dt
$$
so finally we get
$$
c_{j,k}=\sum_{m \in \mathbb{Z}}h_{m-2k}c_{j+1,m}
$$
and doing likewise for the wavelet function, we get
$$d_{j,k}=\sum_{m \in \mathbb{Z}}h^{*}_{m-2k}c_{j+1,m}$$
Now, we change gears and consider the space of square summable sequences $\ell^{2}(\mathbb{Z})$. We also will wade out into the murky waters of digital signal processing terminology and notation as will so be ready for some bizarre notation as well. Ok so consider a sequence which we will denote by $x[n]$ which is convolved another sequence $h[n]$ so that we arrive at an output sequence $y[n]$, given by
$$
y[n]=\sum^{N-1}_{k=0}h[k]x[n-k]
$$
and then notice how familiar this looks to our results from the previous derivation if we re-write them as 
$$
c_{j}[k]=\sum_{m}h[m-2k]c_{j+1}[m]
$$ and
$$
d_{j}[k]=\sum_{m}h_{1}[m-2k]c_{j+1}[m]
$$
If you can understand how we get from the theory to this point, you are ready to start reading about filter banks, downsampling, upsampling, etc. Basically, the DWT can be informaly thought of as a bandpass filter with the wavelet scaling functions being the high pass filter and the scaling function coefficients functioning as a low-pass filter. This perspective is much easier to work with (though much less beautiful) than the mathematical theory. 
Also, if you want to work with wavelets in MATLAB, I recommend forking over the 25 or so dollars for the Wavelet toolbox. It will make your life much easier than trying to start from scratch. I hope this answers your question and I skipped and/or omitted a great deal of detail to avoid writing a book so if you are unsure about anything let me know and I'll try my best to clarify.
Edit:
Well, you always have to sum over all of the integer indexes for the value of $k$, to see this, imagine you given a specific resolution, say $j_{0}$, consider the fact that $L^{2}=\mathcal{V}_{j_{0}} \oplus \mathcal{W}_{j_{0}} \oplus \mathcal{W}_{j_{0}+1} \oplus \cdots $ and $f \in L^{2}(\mathbb{R})$ you would have
$$
f(t)=\sum^{\infty}_{k=-\infty} c_{j_{0},k} \, \varphi_{j_{0},k}(t) + \sum^{\infty}_{j=j_{0}} \sum^{\infty}_{k=-\infty} d_{j,k}\psi_{j,k}(t)$$
and for $f \in L^{2}([a,b])$, we have
$$
f(t)=\sum^{2^{j}b-1}_{k=2^{j}a} c_{j_{0},k} \, \varphi_{j_{0},k}(t) + \sum^{\infty}_{j=j_{0}} \sum^{2^{j}b-1}_{k=2^{j}a} d_{j,k}\psi_{j,k}(t)$$
and for a finite sequence $\{a_{n}\}^{N}_{n=-N} \in \ell^{2}$, the values for k are just $k=-N$ summing through to $N$.
