1
$\begingroup$

PyWavelets (1) takes index of max DFT magnitude, (2) adds 1 to it, (3) divides by domain, which is the range of input values to the wavelet ("support"). I don't see where steps 2 and 3 come from, or why 1 doesn't suffice.

MATLAB's docs suggest only step 1, also here - but closed source, can't verify. If "center frequency" is frequency that correlates the strongest with the wavelet at some reference scale, that's by definition what the DFT finds, so why the extra steps?

$\endgroup$
2
  • $\begingroup$ I've really enjoy reading the last few posts/answers you've made. I've been using pywavelets over the last year in a lot of my work. I sadly think that it isn't being maintained as frequently as it used to, so your PR may go unanswered. It would be really great if the team running it can keep up the good work. Wavelets are being used a lot more in the field of research i'm in compared to the last 10 years and I just want to avoid having to use matlab if possible :) haha. $\endgroup$
    – IsmailE
    Oct 31, 2020 at 15:29
  • $\begingroup$ Comment appreciated; I'm currently actively developing ssqueezepy, with a better cwt than any other Python library I know of; expect first pre-release today or tomorrow. More wavelets and other features are also planned. -- Agreed on MATLAB, closed-source stinks. $\endgroup$ Oct 31, 2020 at 15:33

1 Answer 1

3
$\begingroup$

Short version: DFT's bin indices are input length-dependent; "center frequency" is measured relative to the function generating the wavelet. Generated length can vary, so must be accounted for to yield 'correct' $f_c$. This requries injecting information from knowledge of the function, as DFT is blind to "absolute duration"; PyWavelets does this via domain, the range of input values to the function. domain=16, so DFT's bin indices are x16 shifted; undo by dividing.


The answer's in DFT. Take Real Morlet, with the $\sigma >=5$ approximation:

$$ \psi(t) = e^{-t^2/2}\cos(\sigma t) \tag{1} $$

Suppose we evaluate it for $\sigma=5$, $t\in[-8,8]$; what's it look like?

So what is its "center frequency", or frequency of a sinusoid with which the wavelet correlates most strongly? Try tweaking $\beta$ here. Spoiler: sigma; the exponential is almost ignored, just slightly shifting optima locations (but strongest correlate still is $\sigma$):

As there's no $2 \pi$ in $(1)$, 5 is the angular frequency, so $f=5 / 2\pi=0.796$, or $f^{-1}=1.26$; these match the red cosine, which cycles once every 1.26 units, or 0.796 times every 1 unit.

But in the discrete realm, we've no notion of "time" unless explicitly labeled; all we have is samples. What happens if we take the DFT of red? Do we peak at $5$, or $5 / 2 \pi$? Neither; both of these numbers are "per second", and DFT doesn't understand that. All it knows is that a cosine completing one full cycle from left to right ... completes one full cycle from left to right. It calls this $k=1$.

Then we seek $k$ such that the cosine matches most closely with red. Since red makes 0.796 cycles every 1 unit, and there are 16 units, then it'll make 0.796 * 16 = (1/f) * domain = k cycles in total. That's the index of max(abs(DFT)); we just swap f and k in the equation, inputting k to estimate f.

All said of red above also applies to blue; their peaks match. Putting it all together:

f_center = (t[-1] - t[0]) / k_max = domain / k_max

But this looks different from PyWavelets

Because PyWavelets looked ahead. First, we never expect our wavelet to be flat, so we eliminate the possibility of k_max=0 in case the wavelet averages nonzero:

argmax(abs(fft(psi)[1:]))  # _k_max - 1

But this shifts all k, such that original k=1 is now k=0; so we do +1:

index = argmax(abs(fft(psi_pad))[1:])) + 1  # _k_max

What if $k > N/2$, where N == len(psi)? Negative frequency. This is fine, except we seek the counts of cycles, which are non-negative, so we just find its positive. Suppose N=101; then, there is no Nyquist frequency, only max positive and negative frequencies (equal) at k=50, 51.

  • _k_max = 100. This is negative k=1. To obtain positive: N - _k_max = 101 - 100 = 1.
  • _k_max = 50. This is the positive max; do nothing.
  • _k_max = 51. This is the negative max; make it (+): N - _k_max = 101 - 51 = 50.

This explains

if index > len(psi) / 2:
    index = len(psi) - index  # now we get actual `k_max`

And that's all. Except it's not; PyWavelets does domain / (index - 1), and +2 instead of +1, and -index + 2. The behavior is identical for most, but not all, possible len(psi) and _k_max; grid search. Consider (len_psi, _k_max) = (101, 50). pywt finds k_max=51, which is a negative frequency.


How accurate is it?

Not even within two significant figures. Authors may have mistaken higher precision for higher precision, when in reality k_max gets scaled down by the exact amount psi gets scaled up. (Exception for precision < 5, for which estimation greatly degrades):

Can we do better? -- Yes; above right is the result. See here. Opened a PR, also on k_max.


More rigorously, but briefly: center frequency isn't necessarily at max frequency component.

From Mallat 4.1, if $\psi$ is a wavelet, then $|\psi|^2$ is interpreted as a probability distribution per

$$ \int_{-\infty}^{\infty} |\psi (t)|^2 dt = 1 $$

centered in time at

$$ u = \int_{-\infty}^{\infty} t |\psi (t)|^2 dt $$

(i.e. expected value, $\mathbb{E}[|\psi(t)|^2]$) and in frequency ("center frequency") at

$$ \omega_c = \frac{1}{2\pi} \int_{-\infty}^{\infty} \omega |\hat{\psi} (\omega)|^2 d\omega $$

Corresponding variance relations are given, from which we define Heisenberg boxes for time-frequency resolution:

If a wavelet's FT happens to be symmetric about the peak (as with most used wavelets), then center is at peak's location.


Addendum: 'relative to what' is $f_c$? How to interpret?

Relative to "the function", surely, but how to interpret this in terms of signal data we work with? Like with DFT, to find "physical" values, we need to know the sampling frequency, $f_s$. Physical $f_c$ has units of Hertz, so if we know the number of samples that comprise one second, we can relate "frequency relative to samples" (DFT) to "frequency relative to seconds" (more here). Thus:

$$ f_c \left[ \frac{\text{cycles}}{\text{second}} \right] = \left( f_c \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right) $$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.