Synchrosqueezed STFT phase transform

Question

Synchrosqueezed STFT paper, Synchrosqueezing based Recovery of Instantaneous Frequency from Nonuniform Samples, defines STFT phase transform as:

$$ w(u, \xi) = \xi - \frac{\partial_tS_{g}f(u, \xi)}{j2\pi S_gf(u, \xi)} $$

Which differs from synchrosqueezing original paper's - A Nonlinear Squeezing of the Continuous Wavelet Transform Based on Auditory Nerve Models (per CWT) two-fold:

1. The $\xi\ - W$, where original's is simply $W$
2. The $2\pi$, absent in original's

Why the differences; what motivates STFT's definition? And why use 'modified' STFT?

Answer 1

Why 'modified' STFT?

Authors merely state it's "convenient for our purposes", which isn't much. My exploring reveals: extreme numeric instability in reassignment and component inversion of 'standard'. In fact, most STFT implementations do not compute the 'standard' (clarified around $(2)$):

$$ S_f[u, k] = \sum_{n=-\infty}^{\infty} f[n] g[n - u] e^{-j 2\pi k n / N} \tag{1} $$

Default STFT is extremely unstable along columns

The default (i.e. as implemented by most libraries) STFT is extremely sensitive to imperfections or mere finite-ness in the phase transform: for a given column, values in rows plot into wavelet-like oscillations, meaning upon summation (ssqueezing), exclusion of even a single value greatly impacts the final value in a given Tx location. This seems to be a major (but unstated) motivation for ssqueezing according to the Lebesgue measure, rather than invertible summation.

I've verified this with a perfectly set up pure sinusoid (no shadow frequencies); perhaps the remedy is in also using a perfect window, but this many perfections are hardly realistic (or generalizable). Since Lebesgue sets all values in Sx to 1, exclusion of any one value is no big deal, and we see a dramatic improvement:

lebesgue's main strength is also its main weakness: values away from ridge (STFT max row) are (as should be) much smaller, but are counted as equals, which exacerbates noise energy and general perturbations. A simple compromise, then, is to sum abs(Sx) - which seems to work ideally:

Problem is, we may lose invertibility, but not necessarily by a lot (e.g. Griffin-Lim).

Example: superimposed inverted exponential chirp

"Modulated" (modified) sums well both directly and via abs(), while default creates a mess. Key is the last plots, showing Sx values for a single column. This is the direction we synchrosqueeze, or sum, along; summing directly, we alternate between locally close peaks of opposite signs, keeping the sum total small. Adding or removing just one peak can amount to doubling the net-sum and flipping its sign.

To contrary, modulated is a lot smoother; no big deal if a few samples get misassigned. How's it work? For once, the answer's in fft(fft(x)); aside unimportant differences, DFT of DFT of x is x, and x is the windowed signal:

A centered window thus corresponds to fft(x) that has high-frequency oscillations, whereas a window 'centered' about n = 0 (wrapping around) bears low frequencies, thus the smooth plot.

'Standard' STFT is also modulated

DFT is taken always with window centering a data segment, meaning the window is fixed relative to the DFT complex sinusoids (cisoids), or the cisoids move together with the window when they should remain zero-phased at t=0 (Eq. $(1)$).

Difference is, this "modified" variant has a much more favorable shift factor, and is in exact agreement with its claimed mathematical formulation:

$$ S^{\text{mod}}_f[u, k] = \sum_{n=-\infty}^{\infty} f[n] g[n - u] e^{-j 2\pi k (n - u) / N} \tag{2} $$

Answer 2

TL;DR:

1. CWT's $\partial_t$ is taken with respect to wavelet timeshifts, STFT's with respect to window time samples, with CWT's FT taken along rows, and STFT's along columns.
2. The $2\pi$ is per assuming a different 'oscillatory factor' of Fourier Transform's kernel, i.e. $e^{-j2\pi f}$ vs. $e^{-j\omega}$. Numeric implementations can work with either.
3. "Why modified STFT" is a separate question, addressed in second answer.

---

Motivating phase transform

To derive CWT's phase transform, original paper uses $f(t) = A \cos(\Omega t)$ as a motivating example, or for our purposes equivalently $A e^{j\Omega t}$, then works backwards from the result of $\partial_t W_f / W_f$ to define $w$ such that $\Omega$ is extracted. The idea is, this $f(t)$ is the basis function (building block, atomic unit) of the Fourier Transform, so results here will transfer to any other $h(t)$ the more $h(t)$'s (D)FT is "well-behaved" (more specifically, $h(t)$'s "components" obey synchrosqueezing assumptions detailed in followup paper), which spans a large variety of real-world signals.

---

Deriving phase transform

The SSQ STFT paper considers $e^{j\omega t}$ directly, but omits a derivation; I will derive it here for both the 'standard' and 'modified' STFT, along their derivatives, from which the phase transform will arise naturally.

Let $f(t) = e^{j\xi_0 t}$, and assume a real, even-symmetric ($g=g^*,\ \hat{g} = \hat{g^*}$), window $g(t)$. We have:

$$ \begin{align} S_gf(u, \xi) & = \int_{-\infty}^{\infty} x(t) g(t-u) e^{-j\xi t} dt \\ & = \int_{-\infty}^{\infty} e^{j\xi_0t} g(t-u) e^{-j\xi t} dt \\ & = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) \hat{g}(\omega) e^{-ju\omega} d\omega \\ & = \hat{g}(\xi - \xi_0) e^{-ju(\xi - \xi_0)} \\ \end{align} $$

where Plancherel-Parseval's Theorem was applied, and similarly

$$ \begin{align} S^{\text{mod}}_g f(u, \xi) & = \int_{-\infty}^{\infty} x(t) g(t-u) e^{-j\xi (t-u)} dt \\ & = \int_{-\infty}^{\infty} e^{j\xi_0t} g(t-u) e^{-j\xi (t-u)} dt \\ & = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) \hat{g}(\omega) e^{-ju(\omega - \xi)}d\omega \\ & = \hat{g}(\xi - \xi_0) e^{ju\xi_0} \\ \end{align} $$

The derivatives are:

$$ \begin{align} S_{g'} f(u, \xi) & = \int_{-\infty}^{\infty} x(t) g'(t-u) e^{-j\xi t} dt \\ & = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) j\omega\hat{g}(\omega) e^{-ju\omega} d\omega \\ & = j(\xi - \xi_0)\hat{g}(\xi - \xi_0) e^{-ju(\xi - \xi_0)} \\ \end{align} $$

$$ \begin{align} S_{g'}^{\text{mod}} f(u, \xi) & = \int_{-\infty}^{\infty} x(t) g'(t-u) e^{-j\xi (t-u)} dt \\ & = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) j\omega\hat{g}(\omega) e^{-ju(\omega - \xi)}d\omega \\ & = j(\xi - \xi_0)\hat{g}(\xi - \xi_0) e^{ju\xi_0} \\ \end{align} $$

Then, for both cases, to retrieve $\xi_0$, we consider the ratio

$$ \begin{align} \frac{\partial/\partial t( S_{g}f(u, \xi) )}{S_gf(u, \xi)} = \frac{S_{g'}f(u, \xi)}{S_gf(u, \xi)} = j(\xi - \xi_0) \end{align} $$

which leads us to define

$$ \begin{align} & \boxed{ w(u, \xi) = \xi - \frac{S_{g'}f(u, \xi)}{jS_gf(u, \xi)} } \\ &= \xi - \frac{j(\xi - \xi_0)}{j} = \xi_0 \\ \end{align} $$

Note the phase transform being same for both is one consequence of the division "reducing influence of windowing".

---

Verifying derivation

$S_g f(u, \xi)$ expression can be verified with an example, comparing against Wolfram's output (set normalization=1, oscillatory factor=-1):

$$ \text{Let}\ f(t) = e^{jt},\ g(t) = e^{-\pi t^2}, \text{ and thus}\ \hat{g}(\xi) = e^{-\xi^2 / (4\pi)}.\ \text{Then}, \\ $$ $$ \begin{align} S_gf(u, \xi) &= \hat{g}(\xi - 1) e^{-ju(\xi - 1)} \\ &= \exp{\left( -\frac{1}{4\pi}(\xi - 1)^2 - ju(\xi - 1) \right)} \\ &= \exp{\left(\frac{1}{4\pi}( (\xi - 1)(-j4\pi u - \xi + 1 )\right)} \\ \end{align} $$

More completely, ssqueezepy's ssq_stft on a variety of signals yields expected behavior.

---

Why the $2\pi$?

The TL;DR pretty much says it, which can be verified in the original MATLAB Toolbox's paper. In code, with frequency-domain differentiation, the $2\pi$ makes it so we define $j \omega$ to span $[0, ..., \pi]$ and $(-\pi, ..., 0)$, and without uses $0.5$ instead of $\pi$, reinterpreting the derivative's dimensionality as linear cyclic fraction of sampling frequency, rather than radian. See code.

---

Why the $\xi\ -$?

CWT's $\partial W_f$ takes $\partial / \partial b (\psi((t-b)/a))$, where $b = $ timeshift, whereas STFT's takes $\partial / \partial t (g(t))$. This is per how each is computed:

1. CWT takes DFT along rows; for every scale (row), we use frequency-domain convolution to correlate the wavelet with the signal at every timeshift $b$ (column).
2. STFT takes DFT along columns; for every timeshift (column), we apply DFT to correlate the windowed signal with cisoid bases at every frequency (row).

Following the derivations from motivating examples reveals the $\xi \ -$ stemming from this difference.

The definitions can be interchanged or made the same for both; why exactly they differ I'm unsure. My bet is, CWT is best computed for every timeshift (hop_length = 1), and the scales are already sparse, whereas in STFT we often apply large hop_length (small overlap), so convolution theorem is very compute-wasteful.