Initialize repository

.gitignore

@@ -0,0 +1,129 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/

README.md

@@ -0,0 +1,99 @@
+# Synchrosqueezing in Python
+Synchrosqueezing Toolbox ported to Python, authored by Eugene Brevdo, Gaurav Thakur. Original: https://github.com/ebrevdo/synchrosqueezing
+
+**Reviewers needed**; the repository is in a development stage - details below.
+
+## Features
+  - Forward & inverse CWT- and STFT-based Synchrosqueezing
+  - Forward & inverse discretized Continuous Wavelet Transform (CWT)
+  - Forward & inverse discretized Short-Time Fourier Transform (STFT)
+  - Phase CWT & STFT
+
+## Reviewers needed
+An eventual goal is a merged Pull Request to PyWavelets ([relevant Issue](https://github.com/PyWavelets/pywt/issues/258)). Points of review include:
+  1. **Correctness**; plots generated from CWT transforms resemble those in the publications, but not entirely
+  2. **Completeness**; parts of code in the original Github are missing or incorrect, and I'm unable to replace them all
+  3. **Unit tests**; I'm not familiar with Synchrosqueezing itself, so I don't know how to validate its various functionalities
+  4. **Licensing**; unsure how to proceed here; [original's](https://github.com/ebrevdo/synchrosqueezing/blob/master/LICENSE) says to "Redistributions in binary form must reproduce the above copyright notice" - but I'm not "redistributing" it, I'm distributing my rewriting of it
+  5. <s>**Code style**</s>; I'm aware PyWavelets conforms with PEP8 (but I don't), so I'll edit PR code accordingly
+
+## Review To-do:
+
+**Correctness**:
+  - [ ] 1. Example 1
+  - [ ] 2. Example 2
+
+**Completeness**:
+  - [ ] 1. `freqband` in `synsq_cwt_inv` and `synsq_stft_inf` is defaulted to an integer, but indexed into as a 2D array; the two have nearly identical docstrings, and reference the same equation, but the equation appears completely irrelevant to both.
+  - [ ] 2. `quadgk` has been ported as quadpy's [`quad`](https://github.com/nschloe/quadpy/blob/master/quadpy/line_segment/_tools.py#L16) (linked its wrapped function), which does not support infinite integration bounds, and has [trouble](https://github.com/nschloe/quadpy/issues/236) with computing `synsq_stft_inv`'s integral. Needs a workaround.
+  - [ ] 3. As seen in examples, the y-axis shows "scales", not frequencies - and the relation between the two is neither clear nor linear; it also isn't linear w.r.t. `len(t)`, `nv`, or `fs`. Publications show frequencies instead.
+
+**Unit tests**: Whatever suffices for PyWavelets will suffice for me
+
+## Implementation To-do:
+  One checkmark = code written; two = reviewed
+
+| Status | Toolbox name | Repository name | Repository file |
+| --- | --- | --- | --- |
+| [ ] [**x**] | [`synsq_cwt_fw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/synsq_cwt_fw.m) | `synsq_cwt_fwd` | [synsq_cwt.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/synsq_cwt.py) |
+| [ ] [**x**] | [`synsq_cwt_iw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/synsq_cwt_iw.m) | `synsq_cwt_inv` | [synsq_cwt.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/synsq_cwt.py) |
+| [ ] [**x**] | [`synsq_stft_fw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/synsq_stft_fw.m) | `synsq_stft_fwd` | [synsq_stft.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/synsq_stft.py) |
+| [ ] [**x**] | [`synsq_stft_iw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/synsq_stft_iw.m) | `synsq_stft_inv` | [synsq_stft.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/synsq_stft.py) |
+| [ ] [**x**] | [`synsq_squeeze`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/synsq_squeeze.m) | `synsq_squeeze` | [wavelet_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/wavelet_transforms.py) |
+| [ ] [**x**] | [`synsq_cwt_squeeze`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/synsq_cwt_squeeze.m) | `synsq_cwt_squeeze` | [wavelet_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/wavelet_transforms.py) |
+| [ ] [**x**] | [`phase_cwt`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/phase_cwt.m) | `phase_cwt` | [wavelet_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/wavelet_transforms.py) |
+| [ ] [**x**] | [`phase_cwt_num`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/phase_cwt_num.m) | `phase_cwt_num` | [wavelet_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/wavelet_transforms.py) |
+| [ ] [**x**] | [`cwt_fw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/cwt_fw.m) | `cwt_fwd` | [wavelet_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/wavelet_transforms.py) |
+| [ ] [ ]     | [`cwt_iw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/cwt_iw.m) |
+| [ ] [**x**] | [`stft_fw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/stft_fw.m) | `stft_fwd` | [stft_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/stft_transforms.py) |
+| [ ] [**x**] | [`stft_iw`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/stft_iw.m) | `stft_inv` | [stft_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/stft_transforms.py) |
+| [ ] [**x**] | [`phase_stft`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/phase_stft.m) | `phase_stft` | [stft_transforms.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/stft_transforms.py) |
+| [ ] [**x**] | [`padsignal`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/padsignal.m) | `padsignal` | [utils.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/utils.py) |
+| [ ] [**x**] | [`wfiltfn`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/wfiltfn.m) | `wfiltfn` | [utils.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/utils.py) |
+| [ ] [ ] | [`wfilth`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/wfilth.m) |
+| [ ] [**x**] | [`synsq_adm`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/synsq_adm.m) | `synsq_adm` | [utils.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/utils.py) |
+| [ ] [**x**] | [`buffer`](https://www.mathworks.com/help/signal/ref/buffer.html) | `buffer` | [utils.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/utils.py) |
+| [**x**] [**x**] | [`est_riskshrink_thresh`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/est_riskshrink_thresh.m) | `est_riskshrink_thresh` | [utils.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/utils.py) |
+| [**x**] [**x**] | [`p2up`](https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/p2up.m) | `p2up` | [utils.py](https://github.com/OverLordGoldDragon/synchrosqueezing_python/blob/master/synchrosqueezing/utils.py) |
+
+There are more unlisted (see original repo), but not all will be implemented, in particular GUI implementations.
+
+## Differences w.r.t. original
+
+ - **Renamed variables/functions**; more Pythonic & readable
+ - **Removed unused arguments / variables**
+ - **Improved nan/inf handling**
+ - **Added examples**; original repo lacks any examples in README
+ - **Indexing / var references**; MATLAB is 1-indexed, and handles object reference / assignment, and 'range' ops, differently
+ - **Edited docstrings**; filled missing info, & few corrections
+ - **Moved functions**; each no longer has its own file, but is grouped with other relevant functions
+ - **Code style**; grouped parts of code as sub-functions for improved readability; indentation for vertical alignment; other
+ - **Performance**; this repo may work faster or slower, as Numpy arrays are faster than C arrays, but some of original funcs use MEX-optimized code with no Numpy equivalent. Also using dense instead of sparse matrices (see below).
+
+ **Other**:
+  - Dense instead of sparse matrices for `stft_fwd` in [stft_transforms.py](https://github.com/OverLordGoldDragon/ssqueezepy/blob/master/synchrosqueezing/stft_transforms.py), as Numpy doesn't handle latter in ops involved
+
+
+
+## Examples
+
+See [examples.py](https://github.com/OverLordGoldDragon/ssqueezepy/blob/master/examples.py). Links to: [paper [1]](https://sci-hub.se/https://doi.org/10.1016/j.sigpro.2012.11.029), [paper[2]](https://arxiv.org/pdf/0912.2437.pdf). `_inv` methods (reconstruction, inversion) have been omitted as they involve `freqband`.
+
+**EX 1:** Paper [1], pg. 1086
+
+Only real components shown; imaginary are nearly identical, sometimes sign-opposite.
+
+<img src="https://user-images.githubusercontent.com/16495490/73328411-176e5380-4273-11ea-8bc1-502dfd107444.png" width="600">
+<img src="https://user-images.githubusercontent.com/16495490/73328308-b777ad00-4272-11ea-9b7a-89cca0042d1e.png" width="600">
+<img src="https://user-images.githubusercontent.com/16495490/73328370-e9890f00-4272-11ea-9418-9ff741de66f8.png" width="600">
+<img src="https://user-images.githubusercontent.com/16495490/73328515-8350bc00-4273-11ea-883b-1b253d4cd603.png" width="550">
+
+synsq-CWT (`synsq_cwt_fwd`) appears to produce strongest agreement with paper (FIG 4), while none of STFT yield any resemblance of anything in the papers. It's also unclear whether `synsq_squeeze` was used for "synsq" in FIG 4 instead.
+
+**EX 2:** Paper [2], pg. 18
+
+<img src="https://user-images.githubusercontent.com/16495490/73328745-518c2500-4274-11ea-8140-489ee2403118.png" width="600">
+<img src="https://user-images.githubusercontent.com/16495490/73328804-8a2bfe80-4274-11ea-8c04-6a3b29791e52.png" width="600">
+<img src="https://user-images.githubusercontent.com/16495490/73328894-d8410200-4274-11ea-9731-9537ea9588bf.png" width="550">
+
+Similar situation as EX 1; again CWT has close resemblance, and STFT is in a separate reality. The two apparent discrepancies w/ CWT are: (1) slope of the forking incline, steeper in FIG. 3; (2) position of horizontal line, lower in FIG. 3. As for the black lines in FIG 3, they seem to be the (manual) "markings" mentioned under the figure in the paper.

examples.py

@@ -0,0 +1,111 @@
+# -*- coding: utf-8 -*-
+# Papers:
+# [1] https://sci-hub.se/https://doi.org/10.1016/j.sigpro.2012.11.029
+# [2] https://arxiv.org/pdf/0912.2437.pdf
+import numpy as np
+import matplotlib.pyplot as plt
+
+from ssqueezepy import synsq_cwt_fwd, synsq_stft_fwd
+from ssqueezepy import cwt_fwd, stft_fwd
+
+#%%
+def viz_y(y, y1, y2):
+    _, axes = plt.subplots(1, 2, sharey=True, figsize=(11, 3))
+    axes[0].plot(y1)
+    axes[1].plot(y2)
+    plt.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0, hspace=0)
+    plt.show()
+
+    plt.plot(y)
+    plt.gcf().set_size_inches(14, 4)
+    plt.show()
+
+def viz_s(s, sN, s1, s2, s3):
+    _, axes = plt.subplots(1, 3, sharey=True, figsize=(10, 3))
+    axes[0].plot(t, s1)
+    axes[1].plot(t, s2)
+    axes[2].plot(t, s3)
+    plt.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0, hspace=0)
+    plt.show()
+
+    _, axes = plt.subplots(2, 1, sharex=True, figsize=(10, 5))
+    plt.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0, hspace=0)
+    axes[0].plot(t, s)
+    axes[1].plot(t, sN)
+    plt.show()
+
+def _get_norm(data, norm_rel, norm_abs):
+    if norm_abs is not None and norm_rel != 1:
+        raise ValueError("specify only one of `norm_rel`, `norm_abs`")
+
+    if norm_abs is None:
+        vmax = np.max(np.abs(data)) * norm_rel
+        vmin = -vmax
+    else:
+        vmin, vmax = norm_abs
+    return vmin, vmax
+
+def _make_plots(*data, cmap='bwr', norm=None, titles=None):
+    vmin, vmax = norm or None, None
+
+    _, axes = plt.subplots(1, len(data), sharey=True, figsize=(11, 4))
+    plt.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=.1, hspace=0)
+
+    for i, x in enumerate(data):
+        axes[i].imshow(x, cmap=cmap, aspect='auto', vmin=vmin, vmax=vmax)
+        axes[i].invert_yaxis()
+        axes[i].set_title(titles[i], fontsize=14, weight='bold')
+    plt.show()    
+
+def viz_gen(x, cmap='bwr', norm_rel=1, norm_abs=None):
+    vmin, vmax = _get_norm(np.real(x), norm_rel, norm_abs)
+    _make_plots(np.real(x), np.imag(x), cmap=cmap, 
+                norm=(vmin, vmax), titles=("Real", "Imag"))
+
+def viz_gen2(x1, x2, cmap='bwr', norm_rel=1, norm_abs=None):
+    vmin, vmax = _get_norm(np.real(x1), norm_rel, norm_abs)
+    _make_plots(np.real(x1), np.real(x2), cmap=cmap,
+                norm=(vmin, vmax), titles=("s", "sN"))    
+
+OPTS = {'type': 'bump', 'mu':1}
+#%%
+"""Paper [2], pg. 18"""
+t = np.linspace(0, 12, 1000)
+y1 = np.cos(8*t)
+y2 = np.cos(t**2 + t + np.cos(t))
+y  = y1 + y2
+viz_y(y, y1, y2)
+#%%
+Tx_yc, *_ = synsq_cwt_fwd(y, fs=len(t)/t[-1], nv=32, opts=OPTS)
+Tx_ys, *_ = synsq_stft_fwd(t, y, opts=OPTS)
+#%%
+viz_gen(Tx_yc, cmap='bwr', norm_abs=(-5e-5, 5e-5))
+viz_gen(Tx_ys, cmap='bwr', norm_abs=(-5e-5, 5e-5))
+#%%
+"""Paper [1], pg. 1086"""
+t = np.linspace(0, 10, 2048)
+s1 = (1 + .2*np.cos(t)) * np.cos(2*np.pi*(2*t   + .3*np.cos(t)))
+s2 = (1 + .3*np.cos(t)) * np.cos(2*np.pi*(2.4*t + .3*np.sin(t) + .5*t**1.2)
+                                 ) * np.exp(-t/15)
+s3 = np.cos(2*np.pi*(5.3*t + 0.2*t**1.3))
+N = np.random.randn(len(t)) * np.sqrt(2.4)
+s  = s1 + s2 + s3
+sN = s + N
+viz_s(s, sN, s1, s2, s3)
+#%%
+# feed "denoised" (actually noiseless) signal as noted on pg. 1086 of [1]
+Tx_sc,  *_ = synsq_cwt_fwd(s,  fs=len(t)/t[-1], nv=32, opts=OPTS)
+Tx_sNc, *_ = synsq_cwt_fwd(sN, fs=len(t)/t[-1], nv=32, opts=OPTS)
+Tx_ss,  *_ = synsq_stft_fwd(t, s,  opts=OPTS)
+Tx_sNs, *_ = synsq_stft_fwd(t, sN, opts=OPTS)
+#%%
+viz_gen2(Tx_sc, Tx_sNc, cmap='bwr', norm_abs=(-5e-5, 5e-5))
+viz_gen2(Tx_ss, Tx_sNs, cmap='bwr', norm_abs=(-5e-5, 5e-5))
+#%%
+Wx_s,  *_ = cwt_fwd(s,  'bump',        opts=OPTS)
+Wx_sN, *_ = cwt_fwd(sN, 'bump',        opts=OPTS)
+Sx_s,  *_ = stft_fwd(s,  dt=t[1]-t[0], opts=OPTS)
+Sx_sN, *_ = stft_fwd(sN, dt=t[1]-t[0], opts=OPTS)
+#%%
+viz_gen2(Wx_s, Wx_sN, cmap='bwr', norm_abs=(-1.3, 1.3))
+viz_gen2(Sx_s, Sx_sN, cmap='bwr', norm_abs=(-1.3, 1.3))

ssqueezepy/__init__.py

@@ -0,0 +1,8 @@
+from synsq_cwt import *
+from synsq_stft import *
+from wavelet_transforms import *
+from stft_transforms import *
+from utils import *
+
+
+__version__ = '0.80'