Synchrosqueezing Wavelet Transform explanation?

Question

How does Synchrosqueezing Wavelet Transform work, intuitively? What does the "synchrosqueezed" part do, and how is it different from simply the (continuous) Wavelet Transform?

Answer 1

Synchrosqueezing is a powerful reassignment method. To grasp its mechanisms, we dissect the (continuous) Wavelet Transform, and how its pitfalls can be remedied. Physical and statistical interpretations are provided.

If unfamiliar with CWT, I recommend this tutorial. SSWT is implemented in MATLAB as wsst, and in Python, ssqueezepy. (-- All answer code)

---

Begin with CWT of a pure tone:

A straight line in the time-frequency (rather, time-scale) plane, for our fixed-frequency sinusoid over all time - fair. ... except is it a straight line? No, it's a band of lines, seemingly centered about some maximum, likely the "true scale". Zooming,

makes this more pronounced. Let's plot rows within this zoomed band, one by one:

and all superimposed, each for samples 0 to 127 (horizontal zoom):

Notice anything interesting? They all have the same frequency. It isn't particular to this sinusoid, but is how CWT works in correlating wavelets with signals.

It appears much of information "repeats"; there is redundancy. Can we take advantage of this? Well, if we just assume that all these adjacent bands actually stem from one and the same band, then we can merge them into one - and this, in a nutshell, is what synchrosqueezing does. Naturally it's more nuanced, but the underlying idea is that we sum components of the same instantaneous frequency to obtain a sharper, focused time-frequency representation.

Here's that same CWT, synchrosqueezed:

Now that is a straight line.

---

How's it work, exactly?

We have an idea, but how exactly is this mathematically formulated? Motivated by speaker identification and Empirical Mode Decomposition, SSWT builds upon the modulation model:

$$ f(t) = \sum_{k=1}^{K} A_k(t) \cos(\phi_k (t)), \tag{1} $$

where $A_k(t)$ is the instantaneous amplitude and

$$ \omega_k(t) = \frac{d}{dt}(\phi_k(t)) \tag{2} $$

the instantaneous frequency of component $k$, where we seek to find $K$ such "components" that sum to the original signal. More on this below, "MM vs FT".

At this stage, we only have the CWT, $W_f(a, b)$ (a=scale, b=timeshift); how do we extract $\omega$ from it? Revisit the zoomed pure tone plots; again, the $b$-dependence preserves the original harmonic oscillations at the correct frequency, regardless of $a$. This suggests we compute, for any $(a, b)$, the instantaneous frequency via

$$ \omega(a, b) = -j[W_f(a, b)]^{-1} \frac{\partial}{\partial b}W_f(a, b), \tag{3} $$

where we've taken the log-derivative, $f' / f$. To see why, we can show that CWT of $f(t)=A_0 \cos (\omega_0 t)$ is:

$$ W_f(a, b) = \frac{A_0}{4 \pi} \sqrt{a} \overline{\hat{\psi}(a \omega_0)} e^{j b \omega_0} \tag{4} $$

and thus partial-diffing w.r.t. $b$, we extract $\omega_0$, and the rest in (3) gets divided out. ("But what if $f$ is less nice?" - see caveats).

Finally, equipped with $\omega (a, b)$, we transfer the information from the $(a, b)$-plane to a $(\omega, b)$ plane:

$$ \boxed{ S_f (\omega_l, b) = \sum_{a_k\text{ such that } |\omega(a_k, b) - w_l| \leq \Delta \omega / 2} W_f (a_k, b) a_k^{-3/2}} \tag{5} $$

with $w_l$ spaced apart by $\Delta w$, and $a^{-3/2}$ for normalization (see "Notes").

And that's about it. Essentially, take our CWT, and reassign it, intelligently.

---

So where are the "components"? -- Extracted from high-valued (ridge) curves in the SSWT plane; in the pure tone case, it's one line, and $K=1$. More examples; we select a part of the plane and invert over it as many times as needed.

---

Modulation Model vs Fourier Transform:

What's $(1)$ all about, and why not just use FT? Consider a pendulum oscillating with fixed period and constant damping, and its FT:

$$ s(t) = e^{-t} \cos (25t) u(t)\ \Leftrightarrow\ S(\omega) = \frac{1 + j\omega}{(1 + j\omega)^2 + 625} $$

What does the Fourier Transform tell us? Infinitely many frequencies, but at least peaking at the pendulum's actual frequency. Is this a sensible physical description? Hardly (only in certain indirect senses); the problem is, FT uses fixed-amplitude complex sinusoid frequencies as its building blocks (basis functions, or "bases"), whereas here we have a variable amplitude that cannot be easily represented by constant frequencies, so FT is forced to "compensate" with all these additional "frequencies".

This isn't limited to amplitude modulation; the less sinusoidal or non-periodic the function, the less meaningful its FT spectrum (though not always). Simple example: 1Hz triangle wave, multiple FT frequencies. Frequency-modulation suffers likewise; more intuition here.

These are the pitfalls the Modulation Model aims to address - by decoupling amplitude and frequency over time from the global signal, rather than assuming the same (and constant!) amplitude and frequency for all time.

Meanwhile, SSWT - perfection:

---

Is synchrosqueezing magic?

We seem to gain a lot by ssqueezing - an apparently perfect frequency resolution, violating Heisenberg's uncertainty, and partial noise cancellation ("Notes"). How can this be?

A prior. We assume $f(t)$ is well-captured by the $A_k(t) \cos(\phi_k (t))$ components, e.g. based on our knowledge of the underlying physical process. In fact we assume much more than that, shown bit later, but the idea is, this works well on a subset of all possible signals:

Indeed, there are many ways synchrosqueezing can go awry, and the more the input obeys SSWT's assumptions (which aren't too restrictive, and many signals naturally comply), the better the results.

---

What are SSWT's assumptions? (when will it fail?)

This is a topic of its own (which I may post on later), but briefly, the formulation's as follows. Firstly note that we must somehow restrict what $A(t)$ and $\psi(t)$ can be, else, for example, $A(t)$ can simply cancel out the cosine and become any other function. More precisely, the components are to be such that:

More info in ref 2.

---

How would it be implemented? There's now Python code, clean & commented. Regardless, worth noting:

1. For very small CWT coefficients, phase is unstable (just like for DFT), which we work around by zeroing all such coefficients below a given threshold.
2. For any frequency row/bin $w_l$ in SSWT plane, we reassign from $W_f(a, b)$ based on what's closest to $w_l$ according to $\omega (a, b)$, and for log-scaled CWT we use log-distance.

---

Summary:

SSWT is a time-frequency analysis tool. CWT extracts the time-frequency information, and synchrosqueezing intelligently reassigns it - providing a sparser, sharper, noise-robust, and partly denoised representation. The success of synchrosqueezing is based in and explanied by its prior; the more the input obeys assumptions, the better the results.

---

Notes & caveats:

• What if $f$ isn't nice in $\omega(a, b)$ example? Valid question; in practice, the more the function satisfies aforementioned assumptions, the less of a problem this is, as authors demonstrate through various lemmas.
• In the SSWT of damped pendulum, I cheated a little by extending signal's time to $(-2, 6)$; this is only to prevent boundary effects, which is a CWT phenomenon that can be remedied; here's directly 0 to 6.
• Partial noise cancellation? Indeed; see pg 536 of ref 1.
• What's the $a^{-3/2}$ in $(5)$? Synchrosqueezing effectively inverts $W_f$ onto the reassigned plane, using one-integral iCWT.
• "Fourier bad?" My earlier comparison is prone to criticism. To be clear, FT is the most solid and general-purpose basis that we have for a signals framework. But it's not an all-purpose-best; depending on context, other constructions are more meaningful and more useful.

---

Where to learn more?

The refernced papers are a good source, so are MATLAB's wsst and cwt docs and ssqueezepy's source code. I may also write further Q&A's, which you can be notified of by subbing this thread.

---

References:

1. A Nonlinear Squeezing of the CWT Based on Auditory Nerve Models - I. Daubechies, S. Maes. Excellent origin paper with succinct intuitions.
2. Synchrosqueezed Wavelet Transforms: a tool for Empirical Mode Decomposition - I. Daubechies, J. Lu, H.T. Wu. Good followup paper with examples.
3. The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications - G. Thakur, E. Brevdo, et al. Further exploration of robustness properties and implementation details (including threshold-setting).

Answer 2

Low-level intuition can be obtained by inspecting the phase transform, visually. Answer complements and is complemented by this one. (-- Answer code)

---

We consider a pure sinusoidal tone; ideas extend naturally to more complex signals.

Band of lines concentrated about a maximum in CWT, and a perfect line in ssqueezed right about f=8. Next, the phase transform:

(See "Notes" on the 8.067.) A region of equal values spanning the same range of scales as the CWT, as expected (with 'small' values dropped as usual, see "Notes"). This is the core representation to synchrosqueezing: it says all CWT values in the band get remapped (summed) to the same frequency row.

There is, however, one more step: converting from values of frequencies to reassign to to indices to reassign to. We first define our range of frequencies to remap to, then translate every frequency value in w to an index in k:

• Given ssq_freqs ranging from 2 to 48, if w[5, 2] == 5, then k[5, 2] == np.where(ssq_freqs == 5)
• Or if no exact match, then index of ssq_freqs value closest to 5. E.g. [4.71, 4.92, 5.13] --> index=1.
• k.shape == w.shape == Wx.shape

k is thus the index-equivalent of w. I also has a smoothing effect, in that small perturbations in adjacent rows (scales) per discretization imperfections can get mapped to the same k, as if there were no perturbations ("Notes"). The heatmap of k will look identical to that of w in this case (except now (min, max) == (0, len(ssq_freqs) - 1)).

The synchrosqueezing step is now simply:

for i in range(k.shape[0]):
    for j in range(k.shape[1]):
        Tx[k[i, j], j] += Wx[i, j] * norm

Note that we ssqueeze along rows (scales): k[i,] dictates which row k[i,] of Tx does the row i of Wx gets summed to, for every j.

---

SSqueezing, visually

First show k, with y-axis relabeled to use indices instead of scales (arbitrary, helps next step):

Note the heatmap max here is 111, as k doesn't have any greater value, since we're only reassigning to one row/frequency (ssq_freqs[111] == 8.1). We pick a timestep (column) arbitrarily - let it be j=0 - and ssqueeze:

(For full accuracy, we should also be drawing a bunch of blue arrows (lines) from every other row of CWT to row 0 of SSQ CWT - and, be summing Wx (complex) rather than abs(Wx)). This is repeated for every other j, and thus completes the synchrosqueezing of the CWT.

---

Notes

1. First, instead of zeroing the phase transform per w[np.abs(Wx) < eps], I zeroed per w[np.abs(aWx) < np.abs(Wx).mean()], which is a lot more. This is to keep the example clean, as other behavior involves complications (covered in complement), and the contributions of those excluded bands are negligible anyway. Important to note, however, that neither Wx nor w abruptly drop off to zero, but span an exponentially decaying continuum.

2. Second simplification was pretending that w is perfectly constant over the red region; instead setting colors to range within w's magnitude:

Notice the wavy pattern across timesteps, and the small exponential variation across scales. Why wavy? - covered in complement. Again, this variation vanishes in k (in this case).

3. Why 8.067 and not 8.00? -- recall we added an extra sample; this adds a small fraction of an extra cycle, thus raising frequency. Why the extra sample (i.e. endpoint=True)? Because it's the only way to obtain the equivalent of endpoint=False post-padding, and not doing so pollutes the frequency-domain representation with leakage frequencies, thus ruining picturesque demos.