Can't accurately compute C_d for an arbitrary wavelet.
aaren opened this issue · 11 comments
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_dThe 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).
@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.
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 ?
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.
Yeah, it is the same as far as I can tell as the paper version. I can't find an issue with the code.
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.
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/687474703a2f2f75706c6f61642e77696b696d656469612e6f72672f6d6174682f342f342f662f34346666363266356232343461396663653363353662396566653164313432302e706e67What 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).
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).
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).
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:
- s0 sets the scaling for all of the scales
- 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.
I started coding up my own implementation of an inverse continuous wavelet transform. While searching for my own solution to this same problem, I found your repo.
I think the number quoted in the T&C paper is "wrong".
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).
As you say, in the real, live, finite sample size case Cdelta is dependent on N, dj, and s0. If you use T&C eq. (5) to set up your w_k's, dt will affect the numerical integration step size (in the sum), even though it explicitly cancels. To reproduce the T&C result you need to use the exact settings they used.
If you look in their fortran code, you'll see they set:
N = 504
dt = 1/4
dj = 1/4
s0 = 1/4
jtot = 44
These are not the same values stated in T&C section 3.f. They also comment:
Note: for accurate reconstruction and wavelet-derived variance do not pad with zeroes, set s0=dt (for Paul set s0=dt/4), and use a large "jtot" (even though the extra scales will be within the cone of influence).
For plotting purposes, it is only necessary to use s0=2dt (for Morlet) and "jtot" from Eqn(10) Torrence&Compo(1998).
It appears that in practice T&C are over computing for numerical stability.
Anyways, using their values, I can reproduce Cdelta = 0.77574... for the w0=6 Morlet. Following thier comment and setting s0 = 1/16 for the m=4 Paul wavelet, I get Cdelta = 1.1302.... Actually, I used jtot = 35, because my scale generator sets that. I assume that's the cause of my small difference with T&C's Table 2.
The main point is that T&C's quoted numbers are not general. I think it's best to compute Cdelta and any other empirically determined constants for a particular wavelet basis as needed and ignore the T&C values in their table 2.
Thanks for the nice code! Here I wish to share what I figured out about Cdelta after some struggling and frustration.
The TC98 Cdelta values are indeed resolution dependent. Idealised values of Cdelta can be computed with
def phi_lnk(lnk):
k = np.exp(lnk)
return self.wavelet.frequency(k).real
Cd = scipy.integrate.quad(phi_lnk, -10, 10)[0] / self.wavelet.time(0).real * (2.*np.pi)**0.5 / np.log(2.)
The ln(2) factor arises from the summation of the 2-based self.scales in the Cdelta computation, which has not been normalised in the Cdelta definition in TC98 or this code.
For complex wavelet functions whose Fourier transform is specified to be 0 for k<=0, such as Morlet, the resulting value needs to be further divided by 2.
The resulting Cdelta value for Morlet is 0.7784..., for Marr/Ricker is 3.6163..., both slightly larger from their TC98 value.
Hope this helps!