<p align="center"><img src="https://user-images.githubusercontent.com/16495490/99882586-faa86f80-2c3a-11eb-899c-b3984e98b1c7.png" width="300"></p>

Synchrosqueezing in Python

<!-- [](https://travis-ci.com/OverLordGoldDragon/ssqueezepy) [](https://coveralls.io/github/OverLordGoldDragon/ssqueezepy) -->

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 Python¹, 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 | - | - | -

Questions?

See here.

Examples

1. Signal recovery under severe noise

2. Reconstruction (differentiable CWT)

<img src="https://media.githubusercontent.com/media/OverLordGoldDragon/ssqueezepy/master/examples/imgs/reconstruction.gif" width="700">

Code

3. Medical: EEG

<img src="https://user-images.githubusercontent.com/16495490/99880110-c88f1180-2c2a-11eb-8932-90bf3406a20d.png"> <img src="https://user-images.githubusercontent.com/16495490/150314341-df5c3092-4bef-4895-99ed-2765504329fd.png">

4. Testing suite: CWT vs STFT, reflect-added parallel A.M. linear chirp

<img src="https://github.com/OverLordGoldDragon/ssqueezepy/assets/16495490/c89727db-1bb3-4cf0-ac8c-c524cba75b2d">

5. Ridge extraction: cubic polynom. F.M. + pure tone; noiseless & 1.69dB SNR

<img src="https://user-images.githubusercontent.com/16495490/107919540-f4e5d000-6f84-11eb-9f86-dbfd34733084.png">

More

6. Testing suite: GMW vs Morlet, reflect-added hyperbolic chirp (extreme time-loc.)

<img src="https://github.com/OverLordGoldDragon/ssqueezepy/assets/16495490/8c41d5f2-4bdd-4537-8d82-6d5a5c0315d3">

7. Higher-order GMW CWT, reflect-added parallel linear chirp, 3.06dB SNR

<img src="https://user-images.githubusercontent.com/16495490/107921072-66bf1900-6f87-11eb-9bf5-afd0a6bbbc4d.png">

More examples

Introspection

ssqueezepy is equipped with a visualization toolkit, useful for exploring wavelet behavior across scales and configurations. (Also see explanations and code)

<p align="center"> <img src="https://media.githubusercontent.com/media/OverLordGoldDragon/ssqueezepy/master/examples/imgs/anim_tf_morlet20.gif" width="700"> </p> <img src="https://raw.githubusercontent.com/OverLordGoldDragon/ssqueezepy/master/examples/imgs/morlet_5vs20_tf.png"> <img src="https://user-images.githubusercontent.com/16495490/107297978-e6338080-6a8d-11eb-8a11-60bfd6e4137d.png">

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

<img src="https://raw.githubusercontent.com/OverLordGoldDragon/StackExchangeAnswers/main/SignalProcessing/Q78512%20-%20Wavelet%20Scattering%20explanation/cwt.gif" width="650">

compute phase transform, then combine oscillations with a shared rate

<img src="https://user-images.githubusercontent.com/16495490/150680428-4a651934-85c6-45e8-8a19-c9b4165e5381.png" width="700"> <hr>

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)

1. Identify abrupt changes in audio - ssq_cwt and ssq_stft used together to solve an ML problem without ML

Feel free to share yours here.

Learning resources

1. Continuous Wavelet Transform, & vs STFT
2. Synchrosqueezing's phase transform, intuitively
3. Wavelet time & frequency resolution visuals
4. Why oscillations in SSQ of mixed sines? Separability visuals
5. 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)

Asking questions

Open an Issue, and follow t