Wavelet expansion of given function (as opposed to DWT of given data)

Question

To illustrate, say I take the following function (the reason for the square root will be apparent soon)

dist = NormalDistribution[0.5,0.1];
f[x_] := Sqrt[PDF[dist,x]];

Then I can grab the Haar wavelet to calculate coefficients like

data = Table[Evaluate@f[x], {x,0,1,1/31}]; (* 32 equally spaced points *)
dwt = DiscreteWaveletTransform[data, HaarWavelet[], 5];

If I then take squares of coefficients, because of

$$ \sum _{j=J_o}^{J_1} \sum _{k\in \mathbb{Z}} \psi _{j,k}^2+\sum _{z\in \mathbb{Z}} \varphi _{J_0,k}^2\simeq \int f^2 (x) \, dx = \left\| f\right\| ^2 $$

I should get $\simeq 1.0$ as $f$ is, in this case, the square root of a probability density function. And it kind-of works; to wit:

Total[Flatten[Last[#] & /@ dwt[Automatic]]^2]/31

returns 1.0 (dwt[Automatic] gives you the coefficients for the inverse transform only; see here). This should be very close to NIntegrate[Evaluate@f[x]^2, {x, 0, 1}], which is 0.999999.

However, note that "/31"! Where did that come from? This is: that 1/31 does not appear in the Parseval's identity but it is needed for the calculation to work (!)

And if I change HaarWavelet[] by DaubechiesWavelet[], say, then this data sampling no longer works.

I am even taking care of the fact that the Daubechies family has bigger support and for that I am using the following function to determine the number of points required as:

dataSize[wave_, 1] := 
 Length[WaveletFilterCoefficients[wave, "PrimalLowpass"]];
dataSize[wave_, n_Integer /; n > 1] := 
 2 dataSize[wave, n - 1] - 
  Length[WaveletFilterCoefficients[wave, "PrimalLowpass"]] + 2;

What I am doing wrong?

(Why this? You may ask. For example, one could use wavelet expansions to calculate the norm (in Hilbert space) between two functions, or a function and some other approximation; instead of attempting NIntegrate. By the way, I am actually interested in dimensions greater than 1)

Note: I am only interested in orthogonal wavelets; Parseval's formula should be just ok as far as I know.

UPDATE

I modified dist to NormalDistribution[0.5,0.5]. In this case, one needs a longer interval {x,-2,3} to effectively cover $f$. In this case, the squared coefficients are added like 5Total[Flatten[Last[#] & /@ dwt[Automatic]]^2]/31 (note the factor 5.)

This points me to the fact that the $x_i$ sample points are usually considered as in the [0,1] interval (if I remember correctly), therefore the final answer needs to be rescaled by 5. Still... work in progress.

Answer 1

I found the issue and the reason.

First, the reason for that 1/31. As explained in documentation for DiscreteWaveletTransform (see "Properties & Relations")

The energy norm is conserved for orthogonal wavelet families:
In[1]:= data = RandomReal[1, {100}];
In[2]:= dwt = DiscreteWaveletTransform[data, Padding -> 0.];
In[3]:= Norm[data] == Norm[Flatten[Last /@ dwt[Automatic]]]
Out[3]= True

Note that in In[3], Norm[data] is $\sqrt{\sum d_i^2}$. Given that $d_i~=~f(x_i)$ is a sampling point for $f$, $\int f^2(x) \delta x$ ~ $\sum f^2(x_i) \Delta x$ = $\sum d_i^2 \Delta x$; as in my case $\Delta x = \frac{1}{31}$ this explains the factor I didn't understand.

The other was an issue calculating

Total[Flatten[Last[#] & /@ dwt[Automatic]]^2]/31

when the wavelet was not Haar's wavelet. The reference above also solves this. Note that Padding -> 0. option. The default padding is periodic, which treats the end points as if they were in a circle, not good in this case. This option therefore ensures the transform does not extend the data as periodic beyond the provided points (maybe my dataSize function is incorrect, not sure). With Padding -> 0. in DiscreteWaveletTransform and using $\Delta x$, it works.

Here another example, with a bigger region of interest, [-2,3] and using Daubechies wavelet with n=2.

Module[
 {dist, wave, f, data, dwt, level = 7},
 wave = DaubechiesWavelet[2];
 dist = NormalDistribution[0.5, 0.5];
 f[x_] := Sqrt[PDF[dist, x]];
 data = Table[f[x], {x, -2, 3, 5/(2^level - 1)}];
 Print[Plot[Evaluate@f[x], {x, -2, 3}, PlotRange -> {0, 1}, 
   PlotLabel -> Style["f", FontSlant -> "Italic"]]];
 Print[ListPlot[data, PlotLabel -> "Sampling"]];
 Print["Lenght:", Length[data]];
 Print["Aprox. integral from data: ", 5 Total[data^2 ]/2^level];
 dwt = DiscreteWaveletTransform[data, wave, level, Padding -> 0.];
 Print["NIntegrate: ", NIntegrate[Evaluate@(f[x]^2), {x, -2, 3}]];
 Print["Aprox. integral from coefficients: ", 
  5 Total[Flatten[Last[#] & /@ dwt[Automatic]]^2]/2^level];
 ]

which produces

Note in this case $\Delta x$ is $\frac{5}{\text{num_samples}}$ = $\frac{5}{2^{level}}$.