1
$\begingroup$

This question came up in the context of the ssqueezepy library. As a basic experiment I did compute the synchrosqueezed wavelet transform of three basic signals:

  1. A sine of 440 Hz.
  2. A sine of 880 Hz.
  3. A signal that mixes (1) and (2).

The result looks like this:

enter image description here

Full reproduction code: here

Basically the transform manages to perfectly localize signals (1) and (2), but for the mixed signal there is a surprising oscillation pattern. Zooming in a bit:

enter image description here

In general the issue depends on the choice of mu of the underlying wavelet. In this particular example, increasing mu slightly can help to make the transform of the mixed signal non-oscillating.

Nonetheless I'm wondering why synchrosqueezing leads to this oscillating interpretation of the signal in the first place?

$\endgroup$
1
$\begingroup$

This was interesting to figure out. The key lies in the phase transform, and how CWT interacts with own derivative upon insufficient component separation. Relevant are, and I'll be answering, the following:

  1. What causes the wavy pattern?
  2. Is the wavy pattern truly a sine, or a lookalike?
  3. Why is only one band wavy, not both?
  4. What can be done about it?

Answer split in two to separate meat from supplement, but much of insight's in latter; sections labeled in recommended reading order.

TL;DR: skip to Conclusion. Answer code.


0. How does the phase transform work?

A helpful requisite for following sections, covered here.


1. What causes the wavy pattern?

Begin by visualizing the phase transform, w:

enter image description here

The lower frequency's CWT and SSQ are nearly that of a pure tone, but the higher frequency is all wiggly. Let us plot one such wiggly row:

enter image description here

This sure looks sinusoidal. But is it? Stay tuned.

First: is the wavy SSQ explained by the w? Absolutely: w is wavy across timesteps, meaning for every column, we reassign CWT rows to a different SSQ row, in a wavy/sinusoidal pattern (see [complement] on reassignment). Our answer thus lies in understanding why w is wavy to begin with. Second answer begins here.


4. Inseparability: pre-transform perspective

Recall, the CWT and dCWT operations both take place in the frequency domain; we reach the lowest level of analysis by inspecting it directly. With xh = fft(pad(x)), psih = freq-domain wavelet, and dpsih = freq-domain wavelet with the derivative (* 1j * xi / dt), we plot both, overlapped and rescaled to fit in same graph, as well as the result of the multiplication - for our same row, and one on higher scale:

enter image description here

First, note for our original row scales[80] = 5.08, we indeed have a sum of two pure cisoids (analytic wavelet -> all negative frequencies are zero, and Morlet has pure-real DFT; cosine input -> pure real DFT; thus a DFT coefficient multiplies its complex exponential basis without cancellation), one at higher frequency having much greater amplitude.

Next, compare the cisoid amplitudes for psih vs dpsih; higher frequency dominates more for dpsih than for psih for both example scales. Coincidence? Well, they both do * xh, thus we expect the dpsih wavelet itself to skew toward right. Plotting,

enter image description here

this is confirmed. Does this happen for every scale? Yes: * xi = * linspace(0, pi, N) - taking the derivative involves multiplying by a line, thus skewing the entire wavelet toward higher frequencies, also breaking even symmetry (though not a lot).

All puzzle pieces are in. The skewness creates a difference in relative amplitudes at the two frequencies w.r.t. the non-derivative. The general form for dWx / Wx is then (still assuming purely real fft(Wx)):

$$ \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] $$

($+\pi/2$ per the derivative). Let's first inspect interesting cases numerically, then derive an exact relation as proof (second answer).


8. What can be done about it?

Strictly speaking, nothing, unless the wavelet decays completely to zero eventually (untrue for many wavelets); as long as there's more than one frequency in the entire DFT of the input, there will be some form of interference. However, practically, such interference can be made negligible by:

  1. Ensuring bands are sufficiently separated
  2. Reducing wavelet frequency width (for Morlet, increase mu)

Let's inspect Morlet mu=5 vs mu=20:

enter image description here

mu=20 is clearly narrower; to get a better idea, single out a few rows:

enter image description here

Overlapped at same center frequency (~same row):

enter image description here

We thus predict mu=20 will show no wiggliness we can see. In fact, mu=20 is overkill, we'll settle for mu=10:

enter image description here

Can we also guess how to make both bands wiggly for small mu? Note, increasing scales lowers the wavelet's center frequency (shifts it left), but also reduces its frequency width; thus our wiggle odds improve toward lower scales (or higher frequencies). This turned out easier said than done; after some experimentation, I've had limited success, but maybe it can be done better:

enter image description here


9. Conclusion

  1. What causes the wavy pattern? The frequency-domain wavelet and its derivative intercept the input bands at slightly different ratios per latter being skewed. The resulting phase transform is (imag. part of) a ratio of sums of complex sinusoids, which is of form $(\cos(\phi) + k_1)/(\cos(\phi) + k_2)$. The reassignment pattern dictated by w can then have a sine-like pattern across rows, which manifests as a sine-looking synchrosqueezed CWT.
  2. Is the wavy pattern truly a sine, or a lookalike? Latter; it can never be exactly sine, but close to. The pattern is characterized by $(\cos(f_1 - f_2) + k_1)/(\cos(f_1 - f_2) + k_2)$, and has a period of $|f_1 - f_2|^{-1}$, where $f_1,f_2=$ input frequencies.
  3. Why is only one band wavy, not both? With increasing scale (and decreasing center frequency), the frequency-domain wavelet moves to the left while shrinking in width; the resulting cisoid proportions work out such that dominantly sine reassignment pattern is toward the higher input frequency, where the wavelet is wider, but the lower one isn't devoid of such a pattern either.
  4. What can be done about it? Separate input frequencies more, or reduce wavelet frequency width across all scales (for Morlet wavelet, raise mu).
$\endgroup$
1
$\begingroup$

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:

enter image description here

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:

enter image description here

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

enter image description here

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:

enter image description here

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

enter image description here

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

  1. Constant: happens if $a/b=c/d$.
  2. 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.
  3. Constant + undefined singularities: happens if $c=d$.
  4. Exact sinusoid: impossible. The reassignment pattern will never be exactly sine.

enter image description here

#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):

enter image description here

Still, earlier we saw the row of w looking very sinusoidal; how far off is it? Plotting overlayed:

enter image description here

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:

enter image description here

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:

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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