How is wavelet center frequency computed?

Question

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?

Answer 1

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) $$