Can't accurately compute C_d for an arbitrary wavelet.

[aaren]

Example code:

from wavelets import WaveletAnalysis
wa = WaveletAnalysis()

# hard coded Morlet value
print wa.C_d
print wa.wavelet.C_d

# override hard coding
wa.wavelet.C_d = None
# Now C_d is explicitly computed
print wa.C_d

The explicitly computed version is giving me a value of (0.13708 + 0j), rather than the 0.776 that we expect for a Morlet with w0 = 6 (which is the default).

[aaren]

@nabobalis @Cadair see section 3.i of Torrence and Compo for the maths.
I've implemented the function in wavelets.WaveletAnalysis.compute_Cdelta here
but maybe I've got something wrong.

[nabobalis]

From what I can see either, the value in the paper is wrong or either Y_00 or W_D are returning the wrong values.

For the morlet, is the frequency_rep function calculating the fourier transform of the Morlet wavelet as is given http://en.wikipedia.org/wiki/Morlet_wavelet ?

[aaren]

It is implemented as in Table 1 of TC98 (paper here if you need it. ) I'm not sure how this maps to the wikipedia version.

[nabobalis]

Yeah, it is the same as far as I can tell as the paper version. I can't find an issue with the code.

[aaren]

Are we talking about the same thing? I'm comparing this code, i.e.

def frequency_rep(self, w, s=1.0):
        """Frequency representation of morlet.

        s - scale
        w - angular frequency
        """
        x = w * s
        # heaviside mock
        Hw = 0.5 * (np.sign(x) + 1)
        return np.pi ** -.25 * Hw * np.exp((-(x - self.w0) ** 2) / 2)

with the fourier transform on wikipedia:

What I can't see is how the heaviside function is represented on wikipedia.

[nabobalis]

Yeah, I do mean that. The last part of the equation vanishes for sigma
greater than 5. So it does simplify.
I guess I shouldn't be going to wikipedia anyway.

Regardless, the code follows the paper. So I'm not sure what the problem is
at the moment.

On 31 August 2014 13:51, aaren notifications@github.com wrote:

Are we talking about the same thing? I'm comparing this code
https://github.com/aaren/wavelets/blob/master/wavelets/wavelets.py#L207,
i.e.

def frequency_rep(self, w, s=1.0):
"""Frequency representation of morlet.
s - scale w - angular frequency """
x = w * s
# heaviside mock
Hw = 0.5 * (np.sign(x) + 1)
return np.pi ** -.25 * Hw * np.exp((-(x - self.w0) ** 2) / 2)

with the fourier transform on wikipedia:

[image: wikipedia version]
https://camo.githubusercontent.com/98674993726b55c71cecd0dab15e2eb5131a4ebc/687474703a2f2f75706c6f61642e77696b696d656469612e6f72672f6d6174682f342f342f662f34346666363266356232343461396663653363353662396566653164313432302e706e67

What I can't see is how the heaviside function is represented on wikipedia.

—
Reply to this email directly or view it on GitHub
#1 (comment).

[aaren]

I made a mistake in the code. I'm still not getting the exact same value as TC98, but at least it's a lot closer (getting 0.724 now).

[aaren]

The actual computation of Cdelta is a double summation of the fourier wavelet over scale and frequency. This is sensitive to the values of N, dj and s0 (dt cancels out).

It looks like N just needs to be big enough (>1000), whereas dj and s0 can be a bit smaller than default.

s0 is the smallest scale. The default choice for s0 is to find where this corresponds to the equivalent fourier period being equal to 2*dt, which comes out as ~1.9. Setting s0=1 improves the Cdelta calculation a lot, giving us Cdelta = 0.778 (compared with 0.776 in TC98).

[aaren]

08baf66 nearly fixes the issue, giving Cdelta = 0.776 with s0 = 1. However setting s0 to lower than the default (based on the equivalent fourier period) changes the wavelet variance (TC98 equation 14), meaning that variance (energy) is no longer conserved through the wavelet transform.

The wavelet variance is dependent on the value of s0 for two reasons:

1. s0 sets the scaling for all of the scales
2. s0 sets the uppermost scale sJ via J = log2(N dt / s0) / dj

The scales are defined as sj = s0 * 2 ^ (j * dj), j=0, 1, ..., J (TC98 eq 9-10).

We are told to choose s0 such that the equivalent fourier period is about 2 * dt.