This answer delves deeper into low-level aspects of the phase transform to better understand the wavy phenomenon; complements main answer.
2. How wavy is w?
Recall, the phase transform is: imag(dWx / Wx / (2*pi)). Let's focus on dWx / Wx, and assume frequency-domain differentiation, i.e.
Wx = ifftshift(ifft(Psih * xh))
dWx = ifftshift(ifft(Psih * xh * 1j * xi / dt))
Plot that same row as earlier, of each Wx and dWx, and their ratio:
Individually, they seem to be complex sinusoids of different amplitudes, and in ratio, they produce real and imaginary components that also appear sinusoidal, but with differing offsets (not we discard the real part anyway). What's really happening?
Suppose for a moment that these are complex sinusoids; then, in general, we have:
$$
\begin{align}
\frac{a \exp{(j(\omega_1 n/N + \phi_{01})})}
{b \exp{(j(\omega_2 n/N + \phi_{02}))}}
&= \left|\frac{a}{b}\right| \exp{\left(j \left( \frac{n}{N}(\omega_1 - \omega_2) + (\phi_{01} - \phi_{02}) \right) \right)} \\
&= \left|\frac{a}{b}\right| \exp{\left(j \left( \frac{n}{N} \Delta \omega + \Delta \phi \right) \right)} \tag{1}
\end{align}
$$
The result is another complex sinusoid, whose amplitude is ratio of the two's, and frequency and phase offset the difference of the two's. But a cisoid's real and imaginary components never have different means; this is enough to rule Wx and dWx out as pure cisoids. Then what are they? Let's plot their moduli, and rescale to plot on same graph:
Interestingly, |Wx| and |dWx| themselves appear sinusoidal, characteristic of an amplitude-modulated cisoid. Note also the moduli are of same frequency, but different amplitudes and offsets (means). The ratio of such sines is another sine - and since the real component is much smaller than imaginary, majority of w (MAPE = .03%) is simply:
$$
\left| \frac{a}{b} \right| \sin\left(\frac{n}{N}\Delta \omega + \Delta \phi \right) =
\left| \frac{a}{b} \right| \sin (\Delta \phi) =
\left| \frac{a}{b} \right| \approx w \tag{2}
$$
i.e. $|a/b|$ ($\Delta \omega = 0$ if the two cisoids are of same frequency, and $\Delta \phi = \pi/2$ from differentiation). But why are the moduli of the same frequency and different amplitudes and offsets? Why are they sine-looking at all?
3. AM sinusoids
Necessary background; what's a cisoid with sinusoidal amplitude? A cisoid is just a sinusoid with real/imag component 90-deg out of phase, so it suffices to look at just a sinusoid. It'll be a carrier wave at frequency $f_c$, modulated by a modulator (envelope) at another frequency $f_m$:
Now, a central relation underpinning this mystery and much of synchrosqueezing: AM sine = sum of sines:
$$
\boxed{ \cos(A)\cos(B) = [\cos(A - B) + \cos(A + B)]/2 } \tag{3}
$$
I'll refer back to this. For now, we make a prediction: each cisoid, Wx and dWx, isn't pure, but is a sum of cisoids. First answer continues.
5. Two-band phase transform, numerically
We explore Eq. $(4)$ by inspecting interesting cases numerically:
$$
\frac{dW_x}{W_x} = \frac{a\exp{(j(\phi_1 + \pi/2))} + b\exp{(j(\phi_2 + \pi/2))}}{c\exp{(j\phi_1)} + d\exp{(j\phi_2)}}, \tag{4} \\
\phi_i = \omega_i \frac{n}{N}, \ \ n=[0, 1, ..., N - 1]
$$
$|f2| / |f1| = $ modulus ratio, or ratio of cisoid amplitudes:
The exact relation appears to naturally depend on $f_1, f_2$, spanning from what appears to be 2*abs(cos) when f1 = f2, to 1 + cos(...) when f2 >> f1 (and -> 1 as f2 / f1 -> inf). Predicting imag() of ratio of two of these isn't straightforward, but predicting abs() very much is: it's just ratio of respective abs(), i.e. of right side plots. As shown earlier, this ratio is, to a very close approximation, the entire phase transform w.
Moreover, the right-side plots should look familiar: they're similar to the row's plots of Wx and dWx earlier. In fact, that's exactly what they are, just at different frequencies and number of samples. What we've shown is, the sinusoidal moduli at "same" frequencies but different offsets result from the pure-real Dirac DFT coefficients; a modulus is then the ratio of sum of complex sinusoids at different relative amplitudes but same frequencies.
Note that while "different relative amplitudes but same frequencies" does imply a different "mean frequency" (or center frequency) for each cisoid, the ratio doesn't seem to care, instead retaining the same frequency for all $a, b, c, d$ in $(4)$:
But of course, nothing beats a proof.
6. Proof: frequency preservation in phase transform
We can show the imaginary part of the ratio $(4)$, and thus the complete phase transform, to be exactly:
$$
\boxed{w =
\frac{\cos(n/N(\omega_1 - \omega_2))\cdot (ad + bc) + (ac + bd)}
{\cos(n/N(\omega_1 - \omega_2))\cdot 2cd + c^2 + d^2}
} \tag{5}
$$
Simplifying to barebones,
$$
w = k_0 \frac{\cos(\phi) + k_1}{\cos(\phi) + k_2} \tag{6}
$$
Both numerator and denominator are a cosine plus a constant, each individually having the same frequency + dc, and thus in ratio having that same frequency for all $a, b, c, d, \omega_1, \omega_2$. However; "period" is a more apt term, since the frequency here is not that of a sinusoid. This period is exactly:
$$
T_w = |f_2 - f_1|^{-1} \tag{7}
$$
This concludes the proof. Numeric confirmation in next section.
7. Critical behavior
- Constant: happens if $a/b=c/d$.
- Explosions: $w\rightarrow \infty$ as $c\rightarrow d$, if $a \neq b$. Approaching from left, explodes to $+\infty$ if $a < b$, else to $-\infty$, and vice versa from right. Greater $|a-b| \rightarrow$ faster explosion.
- Constant + undefined singularities: happens if $c=d$.
- Exact sinusoid: impossible. The reassignment pattern will never be exactly sine.
#2 is the explanation for dot-like spots in plots of w at boundary between wavy and flat; it's when the even-symmetric wavelet is centered exactly between the two frequencies. Conveniently, the explosion points get removed by default as part of w[np.abs(Wx) < eps]. In our running example, the wavelet was never exactly centered as such, but was still close; let's visualize raw w (no cleanups):
Still, earlier we saw the row of w looking very sinusoidal; how far off is it? Plotting overlayed:
Very close, but difference is vastly more apparent in respective spectra. You can play around with various $a, b, c, d$ in this Desmos, defaulted to give same ratios as w[80], ratio's multiplied by 10 for clarity. Plotting again overlayed vs true sinusoid:
Based on all above, we should also be able to predict that rows of w will retain the same frequency in neighborhood, but vary slightly in amplitude. And surely enough:
