Synchrosqueezing in Python
Synchrosqueezing is a powerful reassignment method that focuses time-frequency representations, and allows extraction of instantaneous amplitudes and frequencies. Friendly overview.
Features
- Continuous Wavelet Transform (CWT), forward & inverse, and its Synchrosqueezing
- Short-Time Fourier Transform (STFT), forward & inverse, and its Synchrosqueezing
- Wavelet visualizations and testing suite
- Generalized Morse Wavelets
- Ridge extraction
- Fastest wavelet transforms in Python1, beating MATLAB
1: feel free to open Issue showing otherwise
Installation
pip install ssqueezepy
. Or, for latest version (most likely stable):
pip install git+https://github.com/OverLordGoldDragon/ssqueezepy
GPU & CPU acceleration
Multi-threaded execution is enabled by default (disable via os.environ['SSQ_PARALLEL'] = '0'
). GPU requires CuPy >= 8.0.0
and PyTorch >= 1.8.0 installed (enable via os.environ['SSQ_GPU'] = '1'
). pyfftw
optionally supported for maximum CPU FFT speed.
See Performance guide.
Benchmarks
Code. Transforms use padding, float32
precision (float64
supported), and output shape
(300, len(x))
, averaged over 10 runs. pyfftw
not used, which'd speed 1-thread & parallel further. Benched on author's i7-7700HQ, GTX 1070.
len(x) -transform |
1-thread CPU | parallel | gpu | pywavelets | scipy | librosa |
---|---|---|---|---|---|---|
10k-cwt | 0.126 | 0.0462 | 0.00393 | 3.58 | 0.523 | - |
10k-stft | 0.108 | 0.0385 | 0.00534 | - | 0.118 | 0.0909 |
10k-ssq_cwt | 0.372 | 0.148 | 0.00941 | - | - | - |
10k-ssq_stft | 0.282 | 0.147 | 0.0278 | - | - | - |
160k-cwt | 2.99 | 1.25 | 0.0367 | 12.7 | 10.7 | - |
160k-stft | 1.66 | 0.418 | 0.0643 | - | 1.93 | 1.38 |
160k-ssq_cwt | 8.38 | 3.16 | 0.0856 | - | - | - |
160k-ssq_stft | 4.65 | 2.48 | 0.159 | - | - | - |
Examples
1. Signal recovery under severe noise
2. Medical: EEG
3. Testing suite: CWT vs STFT, reflect-added parallel linear chirp
4. Ridge extraction: cubic polynom. F.M. + pure tone; noiseless & 1.69dB SNR
5. Testing suite: GMW vs Morlet, reflect-added hyperbolic chirp (extreme time-loc.)
6. Higher-order GMW CWT, reflect-added parallel linear chirp, 3.06dB SNR
Introspection
ssqueezepy
is equipped with a visualization toolkit, useful for exploring wavelet behavior across scales and configurations. (Also see explanations and code)
How's it work?
In a nutshell, synchrosqueezing exploits redundancy of a time-frequency representation to sparsely localize oscillations, by imposing a prior. That is, we assume x
is well-captured by AM-FM components, e.g. based on our knowledge of the underlying process. We surpass Heisenberg's limitations, but only for a subset of all possible signals. It's also akin to an attention mechanism.
Convolve with localized, analytic kernels
compute phase transform, then combine oscillations with a shared rate
Minimal example
import numpy as np
import matplotlib.pyplot as plt
from ssqueezepy import ssq_cwt, ssq_stft
from ssqueezepy.experimental import scale_to_freq
def viz(x, Tx, Wx):
plt.imshow(np.abs(Wx), aspect='auto', cmap='turbo')
plt.show()
plt.imshow(np.abs(Tx), aspect='auto', vmin=0, vmax=.2, cmap='turbo')
plt.show()
#%%# Define signal ####################################
N = 2048
t = np.linspace(0, 10, N, endpoint=False)
xo = np.cos(2 * np.pi * 2 * (np.exp(t / 2.2) - 1))
xo += xo[::-1] # add self reflected
x = xo + np.sqrt(2) * np.random.randn(N) # add noise
plt.plot(xo); plt.show()
plt.plot(x); plt.show()
#%%# CWT + SSQ CWT ####################################
Twxo, Wxo, *_ = ssq_cwt(xo)
viz(xo, Twxo, Wxo)
Twx, Wx, *_ = ssq_cwt(x)
viz(x, Twx, Wx)
#%%# STFT + SSQ STFT ##################################
Tsxo, Sxo, *_ = ssq_stft(xo)
viz(xo, np.flipud(Tsxo), np.flipud(Sxo))
Tsx, Sx, *_ = ssq_stft(x)
viz(x, np.flipud(Tsx), np.flipud(Sx))
#%%# With units #######################################
from ssqueezepy import Wavelet, cwt, stft, imshow
fs = 400
t = np.linspace(0, N/fs, N)
wavelet = Wavelet()
Wx, scales = cwt(x, wavelet)
Sx = stft(x)[::-1]
freqs_cwt = scale_to_freq(scales, wavelet, len(x), fs=fs)
freqs_stft = np.linspace(1, 0, len(Sx)) * fs/2
ikw = dict(abs=1, xticks=t, xlabel="Time [sec]", ylabel="Frequency [Hz]")
imshow(Wx, **ikw, yticks=freqs_cwt)
imshow(Sx, **ikw, yticks=freqs_stft)
Also see ridge extraction README.
Interesting use cases (with code)
- Identify abrupt changes in audio -
ssq_cwt
andssq_stft
used together to solve an ML problem without ML
Feel free to share yours here.
Learning resources
- Continuous Wavelet Transform, & vs STFT
- Synchrosqueezing's phase transform, intuitively
- Wavelet time & frequency resolution visuals
- Why oscillations in SSQ of mixed sines? Separability visuals
- Zero-padding's effect on spectrum
DSP fundamentals: I recommend starting with 3b1b's Fourier Transform, then proceeding with DSP Guide chapters 7-11. The Discrete Fourier Transform lays the foundation of signal processing with real data. Deeper on DFT coefficients here, also 3b1b.
Contributors (noteworthy)
- David Bondesson: ridge extraction (
ridge_extraction.py
;examples/
:extracting_ridges.py
,ridge_extraction/README.md
)
How to cite
Short form:
John Muradeli, ssqueezepy, 2020. GitHub repository, https://github.com/OverLordGoldDragon/ssqueezepy/. DOI: 10.5281/zenodo.5080508
BibTeX:
@article{OverLordGoldDragon2020ssqueezepy,
title={ssqueezepy},
author={John Muradeli},
journal={GitHub. Note: https://github.com/OverLordGoldDragon/ssqueezepy/},
year={2020},
doi={10.5281/zenodo.5080508},
}
References
ssqueezepy
was originally ported from MATLAB's Synchrosqueezing Toolbox, authored by E. Brevdo and G. Thakur [1]. Synchrosqueezed Wavelet Transform was introduced by I. Daubechies and S. Maes [2], which was followed-up in [3], and adapted to STFT in [4]. Many implementation details draw from [5]. Ridge extraction based on [6].
- G. Thakur, E. Brevdo, N.-S. Fuฤkar, and H.-T. Wu. "The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications", Signal Processing 93:1079-1094, 2013.
- I. Daubechies, S. Maes. "A Nonlinear squeezing of the Continuous Wavelet Transform Based on Auditory Nerve Models".
- I. Daubechies, J. Lu, H.T. Wu. "Synchrosqueezed Wavelet Transforms: a Tool for Empirical Mode Decomposition", Applied and Computational Harmonic Analysis 30(2):243-261, 2011.
- G. Thakur, H.T. Wu. "Synchrosqueezing-based Recovery of Instantaneous Frequency from Nonuniform Samples", SIAM Journal on Mathematical Analysis, 43(5):2078-2095, 2011.
- Mallat, S. "Wavelet Tour of Signal Processing 3rd ed".
- D. Iatsenko, P. V. E. McClintock, A. Stefanovska. "On the extraction of instantaneous frequencies from ridges in time-frequency representations of signals".
License
ssqueezepy is MIT licensed, as found in the LICENSE file. Some source functions may be under other authorship/licenses; see NOTICE.txt.