How is wavelet center frequency computed?
Short version: DFT's bin indices are input length-dependent; &quot;center frequency&quot; 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 &quot;absolute duration&quot;; 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 &gt;=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 &quot;center frequency&quot;, 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 &quot;time&quot; 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 &quot;per second&quot;, 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 &gt; 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 &gt; 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 &lt; 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 (&quot;center frequency&quot;) 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 &quot;the function&quot;, surely, but how to interpret this in terms of signal data we work with? Like with DFT, to find &quot;physical&quot; 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 &quot;frequency relative to samples&quot; (DFT) to &quot;frequency relative to seconds&quot; (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)
$$