Merge edfb5f3 into ac2fd4b

.gitignore

@@ -1,3 +1,6 @@
+# local testing outputs
+tests/wavanim*
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

.travis.yml

@@ -29,7 +29,7 @@ before_install:
 script:
   - >
     pycodestyle --max-line-length=85 
-    --ignore=E221,E241,E225,E226,E402,E722,E741,E272,E266,E302,E731,E702,E201,E129,E203,W503,W504 
+    --ignore=E221,E241,E225,E226,E402,E722,E741,E272,E266,E302,E731,E702,E201,E129,E203,E202,W503,W504 
     ssqueezepy
   - pytest -s --cov=ssqueezepy tests/
 

README.md

@@ -15,17 +15,17 @@ Synchrosqueezing is a powerful _reassignment method_ that focuses time-frequency
 ## Features
   - Continuous Wavelet Transform (CWT), forward & inverse, and its Synchrosqueezing
   - Short-Time Fourier Transform (STFT), forward & inverse, and its Synchrosqueezing
-  - Clean code with explanations and learning references
-  - Wavelet visualizations
-
-### Coming soon
+  - Wavelet visualizations and testing suite
   - Generalized Morse Wavelets
-
+  - Ridge extraction
+
+
 ## Installation
 `pip install ssqueezepy`. Or, for latest version (most likely stable): 
 
 `pip install git+https://github.com/OverLordGoldDragon/ssqueezepy`
 
+
 ## Examples
 
 ### 1. Signal recovery under severe noise
@@ -38,6 +38,25 @@ Synchrosqueezing is a powerful _reassignment method_ that focuses time-frequency
 
 <img src="https://user-images.githubusercontent.com/16495490/104537035-9f8b6b80-5632-11eb-9fa4-444efec6c9be.png">
 
+### 3. Testing suite: CWT vs STFT, reflect-added parallel linear chirp
+
+<img src="https://user-images.githubusercontent.com/16495490/107452011-e7ce7880-6b61-11eb-972f-8aa5ea093dc8.png">
+
+### 4. 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">
+
+### 5. Testing suite: GMW vs Morlet, reflect-added hyperbolic chirp (extreme time-loc.)
+
+<img src="https://user-images.githubusercontent.com/16495490/107903903-d9b69880-6f63-11eb-9478-8ead016cf6f8.png">
+
+### 6. 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](https://overlordgolddragon.github.io/test-signals/)
+
+
 ## Introspection
 
 `ssqueezepy` is equipped with a visualization toolkit, useful for exploring wavelet behavior across scales and configurations. (Also see [explanations and code](https://dsp.stackexchange.com/a/72044/50076))
@@ -47,11 +66,12 @@ Synchrosqueezing is a powerful _reassignment method_ that focuses time-frequency
 </p>
 
 <img src="https://raw.githubusercontent.com/OverLordGoldDragon/ssqueezepy/master/examples/imgs/morlet_5vs20_tf.png">
-<img src="https://raw.githubusercontent.com/OverLordGoldDragon/ssqueezepy/master/examples/imgs/morlet_5vs20_hm.png">
+<img src="https://user-images.githubusercontent.com/16495490/107297978-e6338080-6a8d-11eb-8a11-60bfd6e4137d.png">
 
 <br>
 <hr>
 
+
 ## Minimal example
 
 ```python
@@ -69,27 +89,28 @@ def viz(x, Tx, Wx):
 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]
-x = xo + np.sqrt(2) * np.random.randn(N)
+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()
+plt.plot(x);  plt.show()
 
 #%%# CWT + SSQ CWT ####################################
-Twxo, _, Wxo, *_ = ssq_cwt(xo, 'morlet')
+Twxo, Wxo, *_ = ssq_cwt(xo)
 viz(xo, Twxo, Wxo)
 
-Twx, _, Wx, *_ = ssq_cwt(x, 'morlet')
+Twx, Wx, *_ = ssq_cwt(x)
 viz(x, Twx, Wx)
 
 #%%# STFT + SSQ STFT ##################################
-Tsxo, _, Sxo, *_ = ssq_stft(xo)
+Tsxo, Sxo, *_ = ssq_stft(xo)
 viz(xo, Tsxo, np.flipud(Sxo))
 
-Tsx, _, Sx, *_ = ssq_stft(x)
+Tsx, Sx, *_ = ssq_stft(x)
 viz(x, Tsx, np.flipud(Sx))
 ```
 
+
 ## Learning resources
 
  1. [Continuous Wavelet Transform, & vs STFT](https://ccrma.stanford.edu/~unjung/mylec/WTpart1.html)
@@ -101,6 +122,12 @@ viz(x, Tsx, np.flipud(Sx))
 **DSP fundamentals**: I recommend starting with 3b1b's [Fourier Transform](https://youtu.be/spUNpyF58BY), then proceeding with [DSP Guide](https://www.dspguide.com/CH7.PDF) chapters 7-11.
 The Discrete Fourier Transform lays the foundation of signal processing with real data. Deeper on DFT coefficients [here](https://dsp.stackexchange.com/a/70395/50076), also [3b1b](https://youtu.be/g8RkArhtCc4).
 
+
+## Contributors (noteworthy)
+
+ - [David Bondesson](https://github.com/DavidBondesson): ridge extraction (`ridge_extraction.py`; `examples/`: `extracting_ridges.py`, `README.md`)
+
+
 ## References
 
 `ssqueezepy` was originally ported from MATLAB's [Synchrosqueezing Toolbox](https://github.com/ebrevdo/synchrosqueezing), 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].
@@ -111,6 +138,7 @@ The Discrete Fourier Transform lays the foundation of signal processing with rea
   4. G. Thakur, H.T. Wu. ["Synchrosqueezing-based Recovery of Instantaneous Frequency from Nonuniform Samples"](https://arxiv.org/abs/1006.2533), SIAM Journal on Mathematical Analysis, 43(5):2078-2095, 2011.
   5. Mallat, S. ["Wavelet Tour of Signal Processing 3rd ed"](https://www.di.ens.fr/~mallat/papiers/WaveletTourChap1-2-3.pdf).
 
+
 ## License
 
 ssqueezepy is MIT licensed, as found in the [LICENSE](https://github.com/OverLordGoldDragon/ssqueezepy/blob/master/LICENSE) file. Some source functions may be under other authorship/licenses; see [NOTICE.txt](https://github.com/OverLordGoldDragon/ssqueezepy/blob/master/NOTICE.txt).

examples/cwt_higher_order.py

@@ -0,0 +1,30 @@
+# -*- coding: utf-8 -*-
+"""Show CWT with higher-order Generalized Morse Wavelets on parallel reflect-added
+linear chirps, with and without noise, and show GMW waveforms.
+"""
+import numpy as np
+from ssqueezepy import cwt, TestSignals
+from ssqueezepy.visuals import viz_cwt_higher_order, viz_gmw_orders
+
+#%%# CWT with higher-order GMWs #############################################
+N = 1024
+order = 2
+
+tsigs = TestSignals()
+x, t = tsigs.par_lchirp(N=N)
+x += x[::-1]
+
+for noise in (False, True):
+    if noise:
+        x += np.random.randn(len(x))
+    Wx_k, scales = cwt(x, 'gmw', order=range(order + 1), average=False)
+
+    viz_cwt_higher_order(Wx_k, scales, 'gmw')
+    print("=" * 80)
+
+#%%# Higher-order GMWs #######################################################
+gamma, beta, norm = 3, 60, 'bandpass'
+n_orders = 3
+scale = 5
+
+viz_gmw_orders(N, n_orders, scale, gamma, beta, norm)

examples/scales_selection.py

@@ -0,0 +1,73 @@
+# -*- coding: utf-8 -*-
+"""Shows methods to use for CWT scales selection; also see their docstrings."""
+import numpy as np
+from ssqueezepy import ssq_cwt, Wavelet
+from ssqueezepy.visuals import imshow, plot
+from ssqueezepy.utils import cwt_scalebounds, make_scales, p2up
+from ssqueezepy.utils import logscale_transition_idx
+
+#%%# Helper visual method ####################################################
+def viz(wavelet, scales, scaletype, show_last, nv):
+    plot(scales, show=1, title="scales | scaletype=%s, nv=%s" % (scaletype, nv))
+    if scaletype == 'log-piecewise':
+        extra = ", logscale_transition_idx=%s" % logscale_transition_idx(scales)
+    else:
+        extra = ""
+    print("n_scales={}, max(scales)={:.1f}{}".format(
+        len(scales), scales.max(), extra))
+
+    psih = wavelet(scale=scales)
+    last_psihs = psih[-show_last:]
+
+    # find xmax of plot
+    least_large = last_psihs[0]
+    mx_idx = np.argmax(least_large)
+    last_nonzero_idx = np.where(least_large[mx_idx:] < least_large.max()*.1)[0][0]
+    last_nonzero_idx += mx_idx + 2
+
+    plot(last_psihs.T[:last_nonzero_idx], color='tab:blue', show=1,
+         title="Last %s largest scales" % show_last)
+
+#%%# EDIT HERE ###############################################################
+# signal length
+N = 2048
+# your signal here
+t = np.linspace(0, 1, N, endpoint=False)
+x = np.cos(2*np.pi * 16 * t) + np.sin(2*np.pi * 64 * t)
+
+# choose wavelet
+wavelet = 'gmw'
+# choose padding scheme for CWT (doesn't affect scales selection)
+padtype = 'reflect'
+
+# one of: 'log', 'log-piecewise', 'linear'
+# 'log-piecewise' lowers low-frequency redundancy; see
+# https://github.com/OverLordGoldDragon/ssqueezepy/issues/29#issuecomment-778526900
+scaletype = 'log-piecewise'
+# one of: 'minimal', 'maximal', 'naive' (not recommended)
+preset = 'maximal'
+# number of voices (wavelets per octave); more = more scales
+nv = 32
+# downsampling factor for higher scales (used only if `scaletype='log-piecewise'`)
+downsample = 4
+# show this many of lowest-frequency wavelets
+show_last = 20
+
+#%%## Make scales ############################################################
+# `cwt` uses `p2up`'d N internally
+M = p2up(N)[0]
+wavelet = Wavelet(wavelet, N=M)
+
+min_scale, max_scale = cwt_scalebounds(wavelet, N=len(x), preset=preset)
+scales = make_scales(N, min_scale, max_scale, nv=nv, scaletype=scaletype,
+                     wavelet=wavelet, downsample=downsample)
+
+#%%# Visualize scales ########################################################
+viz(wavelet, scales, scaletype, show_last, nv)
+wavelet.viz('filterbank', scales=scales)
+
+#%%# Show applied ############################################################
+Tx, Wx, ssq_freqs, scales, *_ = ssq_cwt(x, wavelet, scales=scales,
+                                        padtype=padtype)
+imshow(Wx, abs=1, title="abs(CWT)")
+imshow(Tx, abs=1, title="abs(SSQ_CWT)")

examples/se_ans0.py

@@ -33,15 +33,15 @@ def cos_f(freqs, N=128, phi=0):
 plt.show()
 
 #%%## Animate rows ###########################################################
-def row_anim(Wxz, idxs, scales, superimposed=False):
+def row_anim(Wxz, idxs, scales, superposed=False):
     mx = np.max(np.abs(Wxz))
     for scale, row in zip(scales[idxs], Wxz):
         if row.max() == Wxz.max():
             plt.plot(row.real, color='r')
         else:
             plt.plot(row.real, color='tab:blue')
         plt.ylim(-1.05*mx, 1.05*mx)
-        if not superimposed:
+        if not superposed:
             plt.annotate("scale=%.1f" % scale, weight='bold', fontsize=14,
                          xy=(.85, .93), xycoords='axes fraction')
             plt.show()
@@ -51,11 +51,11 @@ def row_anim(Wxz, idxs, scales, superimposed=False):
 #%%
 row_anim(Wxz, idxs, scales)
 #%%## Superimpose ####
-row_anim(Wxz, idxs, scales, superimposed=True)
+row_anim(Wxz, idxs, scales, superposed=True)
 #%%## Synchrosqueeze
-Tx, fs, *_ = ssq_cwt(x, wavelet, t=_t(0, 1, N))
+Tx, _, ssq_freqs, *_ = ssq_cwt(x, wavelet, t=_t(0, 1, N))
 #%%
-imshow(Tx, abs=1, title="abs(SSWT)", yticks=fs, show=1)
+imshow(Tx, abs=1, title="abs(SSWT)", yticks=ssq_freqs, show=1)
 
 #%%# Damped pendulum example ################################################
 N, w0 = 4096, 25
@@ -71,7 +71,7 @@ def row_anim(Wxz, idxs, scales, superimposed=False):
 
 #%%# Now SSWT ##
 wavelet = ('morlet', {'mu': 5})
-Tx, fs, *_ = ssq_cwt(s, wavelet, t=t)
+Tx, *_ = ssq_cwt(s, wavelet, t=t)
 #%%# 'cheat' a little; could use boundary wavelets instead (not implemented)
 aTxz = np.abs(Tx)[:, len(t) // 8:]
 imshow(aTxz, abs=1, title="abs(SSWT(s(t)))", show=1)

examples/test_transforms.py

@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+import scipy.signal as sig
+from ssqueezepy import Wavelet, TestSignals
+from ssqueezepy.utils import window_resolution
+
+tsigs = TestSignals(N=512)
+#%%# Viz signals #############################################################
+# set `dft` to 'rows' or 'cols' to also plot signals' DFT, along rows or columns
+dft = (None, 'rows', 'cols')[0]
+tsigs.demo(dft=dft)
+
+#%%# How to specify `signals` ################################################
+signals = [
+    'am-cosine',
+    ('hchirp', dict(fmin=.2)),
+    ('sine:am-cosine', (dict(f=32, phi0=1), dict(amin=.3))),
+]
+tsigs.demo(signals, N=512)
+
+#%%# With `dft` ##################
+tsigs.demo(signals, dft='rows')
+tsigs.demo(signals, dft='cols')
+
+#%%# Viz CWT & SSQ_CWT with different wavelets ###############################
+tsigs = TestSignals(N=512)
+wavelets = [Wavelet(('gmw', {'beta': 5})),
+            Wavelet(('gmw', {'beta': 22}))]
+tsigs.wavcomp(wavelets)
+
+#%%#
+tsigs.wavcomp(wavelets, signals=[('#echirp', dict(fmin=.1))], N=512)
+
+#%%# Viz CWT vs STFT (& SSQ'd) ###############################################
+# (N, beta, NW): (512, 42.5, 255); (256, 21.5, 255)
+N = 512
+n_fft = N
+win_len = n_fft // 2
+tsigs = TestSignals(N=N)
+wavelet = Wavelet(('GMW', {'beta': 21.5}))
+
+NW = win_len//2 - 1
+window = np.abs(sig.windows.dpss(win_len, NW))
+window = np.pad(window, win_len//2)
+window_name = 'DPSS'
+config_str = '\nNW=%s, win_pad_len=%s' % (NW, len(window) - win_len)
+
+# ensure `wavelet` and `window` have ~same time & frequency resolutions
+print("std_w, std_t, harea\nwavelet: {:.4f}, {:.4f}, {:.8f}"
+      "\nwindow:  {:.4f}, {:.4f}, {:.8f}".format(
+          wavelet.std_w, wavelet.std_t, wavelet.harea,
+          *window_resolution(window)))
+#%%
+tsigs.cwt_vs_stft(wavelet, window, signals='all', N=N, win_len=win_len,
+                  n_fft=n_fft, window_name=window_name, config_str=config_str)

requirements-dev.txt

@@ -4,4 +4,5 @@ coverage
 coveralls
 pytest
 pytest-cov
-pycode
+pycode
+librosa

setup.py

@@ -1,3 +1,21 @@
+# -*- coding: utf-8 -*-
+#
+# Copyright © OverLordGoldDragon
+# Licensed under the terms of the MIT License
+# (see ssqueezepy/__init__.py for details)
+
+"""
+ssqueezepy
+======
+
+Synchrosqueezing, wavelet transforms, and time-frequency analysis in Python
+
+ssqueezepy features time-frequency analysis written for performance, flexibility,
+and clarity. Included are Continuous Wavelet Transform (CWT), Short-Time Fourier
+Transform (STFT), CWT & STFT synchrosqueezing, Generalized Morse Wavelets,
+visualizations, a signal testing suite, and automatic ridge extraction.
+"""
+
 import os
 import re
 from setuptools import setup, find_packages
@@ -33,12 +51,12 @@ def find_version(*file_paths):
     author="OverLordGoldDragon",
     author_email="16495490+OverLordGoldDragon@users.noreply.github.com",
     description=("Synchrosqueezing, wavelet transforms, and "
-                 "time-frequency analysis in Python "),
+                 "time-frequency analysis in Python"),
     long_description=read_file('README.md'),
     long_description_content_type="text/markdown",
     keywords=(
-        "signal-processing python synchrosqueezing wavelet-transform cwt "
-        "time-frequency time-frequency-analysis"
+        "signal-processing python synchrosqueezing wavelet-transform cwt stft "
+        "morse-wavelet ridge-extraction time-frequency time-frequency-analysis"
     ),
     install_requires=get_requirements('requirements.txt'),
     tests_require=["pytest>=4.0", "pytest-cov"],