Merge c82be23 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,18 @@ 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. 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">
+
+[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 +59,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 +82,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]  # superimpose 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 +115,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`)
+
+
 ## 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 +131,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/extracting_ridges.py

@@ -0,0 +1,122 @@
+# -*- coding: utf-8 -*-
+"""Authors: David Bondesson, OverLordGoldDragon
+
+Ridge extraction on signals with varying time-frequency characteristics.
+"""
+import numpy as np
+import scipy.signal as sig
+from ssqueezepy import ssq_cwt, ssq_stft, extract_ridges, TestSignals
+from ssqueezepy.visuals import plot, imshow
+
+import matplotlib
+matplotlib.use('template')  # suppress plots when Spyder unit-testing
+
+#%%## Visual methods #########################################################
+def viz(x, Tf, ridge_idxs, yticks=None, ssq=False, transform='cwt', show_x=True):
+    return
+
+    if show_x:
+        plot(x, title="x(t)", show=1,
+             xlabel="Time [samples]", ylabel="Signal Amplitude [A.U.]")
+
+    if transform == 'cwt' and not ssq:
+        Tf = np.flipud(Tf)
+        ridge_idxs = len(Tf) - ridge_idxs
+
+    ylabel = ("Frequency scales [1/Hz]" if (transform == 'cwt' and not ssq) else
+              "Frequencies [Hz]")
+    title = "abs({}{}) w/ ridge_idxs".format("SSQ_" if ssq else "",
+                                             transform.upper())
+
+    ikw = dict(abs=1, cmap='jet', yticks=yticks, title=title)
+    pkw = dict(linestyle='--', color='k', xlabel="Time [samples]", ylabel=ylabel,
+               xlims=(0, Tf.shape[1]), ylims=(0, len(Tf)))
+
+    imshow(Tf, **ikw, show=0)
+    plot(ridge_idxs, **pkw, show=1)
+
+
+def tf_transforms(x, t, wavelet='morlet', window=None, padtype='wrap',
+                  penalty=.5, n_ridges=2, cwt_bw=15, stft_bw=15,
+                  ssq_cwt_bw=4, ssq_stft_bw=4):
+    kw_cwt  = dict(t=t, padtype=padtype)
+    kw_stft = dict(fs=1/(t[1] - t[0]), padtype=padtype)
+    Twx, Wx, ssq_freqs_c, scales, *_ = ssq_cwt(x,  wavelet, **kw_cwt)
+    Tsx, Sx, ssq_freqs_s, Sfs, *_    = ssq_stft(x, window,  **kw_stft)
+
+    ckw = dict(penalty=penalty, n_ridges=n_ridges, transform='cwt')
+    skw = dict(penalty=penalty, n_ridges=n_ridges, transform='stft')
+    cwt_ridges      = extract_ridges(Wx,  scales,      bw=cwt_bw,      **ckw)
+    ssq_cwt_ridges  = extract_ridges(Twx, ssq_freqs_c, bw=ssq_cwt_bw,  **ckw)
+    stft_ridges     = extract_ridges(Sx,  Sfs,         bw=stft_bw,     **skw)
+    ssq_stft_ridges = extract_ridges(Tsx, ssq_freqs_s, bw=ssq_stft_bw, **skw)
+
+    viz(x, Wx,  cwt_ridges,      scales,      ssq=0, transform='cwt',  show_x=1)
+    viz(x, Twx, ssq_cwt_ridges,  ssq_freqs_c, ssq=1, transform='cwt',  show_x=0)
+    viz(x, Sx,  stft_ridges,     Sfs,         ssq=0, transform='stft', show_x=0)
+    viz(x, Tsx, ssq_stft_ridges, ssq_freqs_s, ssq=1, transform='stft', show_x=0)
+
+#%%# Basic example ###########################################################
+# Example ridge from similar example as can be found at MATLAB:
+# https://www.mathworks.com/help/wavelet/ref/wsstridge.html#bu6we25-penalty
+test_matrix = np.array([[1, 4, 4], [2, 2, 2], [5, 5, 4]])
+fs_test = np.exp([1, 2, 3])
+
+ridge_idxs, *_ = extract_ridges(test_matrix, fs_test, penalty=2.0,
+                                get_params=True)
+print('ridge follows indexes:', ridge_idxs)
+assert np.allclose(ridge_idxs, np.array([[2, 2, 2]]))
+
+#%%# sin + cos ###############################################################
+N, f1, f2 = 257, 5, 20
+padtype = 'wrap'
+penalty = 20
+
+t  = np.linspace(0, 1, N, endpoint=True)
+x1 = np.sin(2*np.pi * f1 * t)
+x2 = np.cos(2*np.pi * f2 * t)
+x = x1 + x2
+
+tf_transforms(x, t, padtype=padtype, penalty=penalty)
+
+#%%# Linear + quadratic chirp ################################################
+N = 257
+penalty = 0.5
+padtype = 'reflect'
+
+t  = np.linspace(0, 10, N, endpoint=True)
+x1 = sig.chirp(t, f0=2,  f1=8, t1=20, method='linear')
+x2 = sig.chirp(t, f0=.4, f1=4, t1=20, method='quadratic')
+x = x1 + x2
+
+tf_transforms(x, t, padtype=padtype, stft_bw=4, penalty=penalty)
+
+#%%# Cubic polynomial frequency variation + pure tone ########################
+N, f = 257, 0.5
+padtype = 'wrap'
+penalty = 2.0
+
+t  = np.linspace(0, 10, N, endpoint=True)
+p1 = np.poly1d([0.025, -0.36, 1.25, 2.0])
+x1 = sig.sweep_poly(t, p1)
+x2 = np.sin(2*np.pi * f * t)
+x = x1 + x2
+
+tf_transforms(x, t, padtype=padtype, stft_bw=4, penalty=penalty)
+
+#%%# Reflect-added linear chirps #############################################
+N = 512
+
+tsigs = TestSignals(N)
+x, t = tsigs.lchirp(N)
+x += x[::-1]
+
+tf_transforms(x, t, cwt_bw=20)
+
+#%%# Parallel F.M. linear chirps  ############################################
+N = 512
+
+tsigs = TestSignals(N)
+x, t = tsigs.par_lchirp(N)
+
+tf_transforms(x, t)

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