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)
$$
ssqueezepy
, with a bettercwt
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$