4
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ You show the code Total[Flatten[Last[#] & /@ dwt[Automatic]]^2]/31 and then go on to say "note that /31! Where did that come from?" Is it not there because you typed it in? If not, where did the whole expression come from? I note 1/31 is the table iteration step you used in an earlier line of code. $\endgroup$ Commented Aug 20, 2014 at 7:34
  • $\begingroup$ Sorry @m_goldberg, what I meant is that that 1/31 does not appear in the Parseval's identity I mentioned, but I had to put it to make the identity work. And unfortunately it is not as simple as the iteration step, as I used wavelets with bigger support and I was not able to find such correction factor. $\endgroup$ Commented Aug 20, 2014 at 11:34

1 Answer 1

Reset to default
3
$\begingroup$

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

Mathematica graphics

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

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.