Merge 3df7d83 into 5f6e85d

CHANGELOG.md

@@ -32,7 +32,7 @@
  - `_infer_scaletype` -> `infer_scaletype`
  - `_integrate_analytic` -> `integrate_analytic`
  - `find_max_scale` -> `find_max_scale_alt`, but `find_max_scale` is still (but a different) function
- 
+
 #### MISC
  - `phase_cwt`: takes `abs(w)` instead of zeroing negatives
  - `wavelet` in `icwt` and `issq_cwt` now defaults to the default wavelet
@@ -47,13 +47,14 @@
 #### FIXES
  - `visuals.wavelet_heatmap`: string `scales` now functional
  - `visuals`: `w` overreached into negative frequencies for odd `N` in `wavelet_tf`, `wavelet_tf_anim`, & `wavelet_heatmap`
- - `icwt`: `padtype` now functional 
+ - `icwt`: `padtype` now functional
 
 #### FILE CHANGES
  - `ssqueezepy/` added: `_gmw.py`, `_test_signals.py`, `ridge_extraction.py`, `configs.ini`, `README.md`
  - `ssqueezepy/` added `utils/`, split `utils.py` into `common.py`, `cwt_utils.py`, `stft_utils.py`, `__init__.py`, & moved to `utils/`.
  - `tests/` added: `gmw_test.py`, `test_signals_test.py`, `ridge_extraction_test.py`
  - `examples/` added: `extracting_ridges.py`, `scales_selection.py`, `ridge_extract_readme/`: `README.md`, `imgs/*`
+ - Created `MANIFEST.in`
 
 
 ### 0.5.5 (1-14-2021): STFT & Synchrosqueezed STFT

MANIFEST.in

@@ -0,0 +1,5 @@
+include README.md
+include LICENSE
+include requirements.txt
+include ssqueezepy/configs.ini
+include ssqueezepy/README.md

README.md

@@ -46,6 +46,8 @@ Synchrosqueezing is a powerful _reassignment method_ that focuses time-frequency
 
 <img src="https://user-images.githubusercontent.com/16495490/107919540-f4e5d000-6f84-11eb-9f86-dbfd34733084.png">
 
+[More](https://github.com/OverLordGoldDragon/ssqueezepy/tree/master/examples/ridge_extraction)
+
 ### 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">
@@ -68,7 +70,6 @@ Synchrosqueezing is a powerful _reassignment method_ that focuses time-frequency
 <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">
 
-<br>
 <hr>
 
 
@@ -110,6 +111,8 @@ Tsx, Sx, *_ = ssq_stft(x)
 viz(x, Tsx, np.flipud(Sx))
 ```
 
+Also see ridge extraction [README](https://github.com/OverLordGoldDragon/ssqueezepy/tree/master/examples/ridge_extraction).
+
 
 ## Learning resources
 

examples/ridge_extract_readme/README.md → examples/ridge_extraction/README.md

@@ -4,7 +4,7 @@ Extracts the `n_ridges` (user selected integer) most prominent frequency ridges
 
 The method is based on a forward-backward greedy path optimization algorithm that penalises frequency jumps similar to the MATLAB function 'tfridge' (https://de.mathworks.com/help/signal/ref/tfridge.html). 
 
-Further information about the particular algorithm (version of eq. III.4 in publication) as well as limitations and comparisson to other ridge extraction schemes can be found in publication [2] (below):
+Further information about the particular algorithm (version of eq. III.4 ref [1]) as well as limitations and comparisson to other ridge extraction schemes can be found in publication [1] (below).
 
 
 
@@ -34,9 +34,9 @@ ridge_idxs = extract_ridges(Tx, scales, penalty, n_ridges=2, bw=4)
 viz(x, Tx, ridge_idxs, ssq_freqs, ssq=True)
 ```
 
-![signal_1](/examples/ridge_extract_readme/imgs/signal_1.png)
-![cwt_ridge_signal_1](/examples/ridge_extract_readme/imgs/cwt_signal_1_ridge.png)
-![ssq_ridge_signal_1](/examples/ridge_extract_readme/imgs/ssq_signal_1_ridge.png)
+![signal_1](/examples/ridge_extraction/imgs/signal_1.png)
+![cwt_ridge_signal_1](/examples/ridge_extraction/imgs/cwt_signal_1_ridge.png)
+![ssq_ridge_signal_1](/examples/ridge_extraction/imgs/ssq_signal_1_ridge.png)
 
 ## 2: Two chirp signals with linear and quadratic frequency variation
 
@@ -62,9 +62,9 @@ ridge_idxs = extract_ridges(Tx, scales, penalty, n_ridges=2, bw=2)
 viz(x, Tx, ridge_idxs, ssq=True)
 ```
 
-![signal_2](/examples/ridge_extract_readme/imgs/signal_2.png)
-![cwt_ridge_signal_2](/examples/ridge_extract_readme/imgs/cwt_signal_2_ridge.png)
-![ssq_ridge_signal_2](/examples/ridge_extract_readme/imgs/ssq_signal_2_ridge.png)
+![signal_2](/examples/ridge_extraction/imgs/signal_2.png)
+![cwt_ridge_signal_2](/examples/ridge_extraction/imgs/cwt_signal_2_ridge.png)
+![ssq_ridge_signal_2](/examples/ridge_extraction/imgs/ssq_signal_2_ridge.png)
 
 ## 3: One sweep signal and one constant frequency signal
 
@@ -91,9 +91,9 @@ ridge_idxs = extract_ridges(Tx, scales, penalty, n_ridges=2, bw=2)
 viz(x, Tx, ridge_idxs, ssq=True)
 ```
 
-![signal_3](/examples/ridge_extract_readme/imgs/signal_3.png)
-![cwt_ridge_signal_3](/examples/ridge_extract_readme/imgs/cwt_signal_3_ridge.png)
-![ssq_ridge_signal_3](/examples/ridge_extract_readme/imgs/ssq_signal_3_ridge.png)
+![signal_3](/examples/ridge_extraction/imgs/signal_3.png)
+![cwt_ridge_signal_3](/examples/ridge_extraction/imgs/cwt_signal_3_ridge.png)
+![ssq_ridge_signal_3](/examples/ridge_extraction/imgs/ssq_signal_3_ridge.png)
 
 ## 4: Example of insufficient penalty term 
 
@@ -119,26 +119,27 @@ ridge_idxs = extract_ridges(Wx, scales, penalty=0.5, n_ridges=2, bw=25)
 viz(x, Wx, ridge_idxs)
 ```
 
-![signal_failed_wsst](/examples/ridge_extract_readme/imgs/signal_failed_wsst.png)
-![cwt_ridge_signal_failed_wsst_noPen](/examples/ridge_extract_readme/imgs/cwt_signal_failed_wsst_ridge_pen00.png)
-![cwt_ridge_signal_failed_wsst_penalty05](/examples/ridge_extract_readme/imgs/cwt_signal_failed_wsst_ridge_pen05.png)
+![signal_failed_wsst](/examples/ridge_extraction/imgs/signal_failed_wsst.png)
+![cwt_ridge_signal_failed_wsst_noPen](/examples/ridge_extraction/imgs/cwt_signal_failed_wsst_ridge_pen00.png)
+![cwt_ridge_signal_failed_wsst_penalty05](/examples/ridge_extraction/imgs/cwt_signal_failed_wsst_ridge_pen05.png)
 
 
 ### Practical considerations
 
 ## 1. Ridge extraction on ssq_cwt
 
-The idea of a syncrosqueezed time-frequency representation is to improve upon the localization process and more accurately recover the frequency components compared to inverting the CWT over the entire time-scale plan [1]. However, whether ridge extraction is improved on a syncrosqueezed transform is still uncertain and may require further considerations [2] and expertise from the users on appropriate parameter tuning for stable and improved results.  
+The idea of a syncrosqueezed time-frequency representation is to improve upon the localization process and more accurately recover the frequency components compared to inverting the CWT over the entire time-scale plane [2]. However, whether ridge extraction is improved on a syncrosqueezed transform is still uncertain and may require further considerations [1] and expertise from the users on appropriate parameter tuning for stable and improved results.  
+
+[1] On the extraction of instantaneous frequencies from ridges in time-frequency representations of signals. D. Iatsenko, P. V. E. McClintock, A. Stefanovska. https://arxiv.org/pdf/1310.7276.pdf <br>
+[2] The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications. Gaurav Thakur, Eugene Brevdo, Neven S. Fuˇckar, and Hau-Tieng Wu. https://arxiv.org/pdf/1105.0010.pdf
 
-[1] The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications. Gaurav Thakur, Eugene Brevdo, Neven S. Fuˇckar, and Hau-Tieng Wu. https://arxiv.org/pdf/1105.0010.pdf <br>
-[2] On the extraction of instantaneous frequencies from ridges in time-frequency representations of signals. D. Iatsenko, P. V. E. McClintock, A. Stefanovska. https://arxiv.org/pdf/1310.7276.pdf
 
 Example 1 particularly highlights how edge-effects may make ridge extraction on the syncrosqueezed transform more unstable. 
 
 ## 2. Signal padding
 
 Particular consideration to choice of wavelet type, signal padding and other input parameters are advised before using the time-frequency ridge extraction. 
-Different padding schems suited to your time signals may be nessecary to handle edge effects in your time frequency maps. Example 2 & 3 highlights an appropriate padding along with usage of wavelet.
+Different padding schems suited to your time signals may be nessecary to handle edge effects in your time frequency maps. Examples 2 & 3 highlight an appropriate padding along with usage of wavelet.
 
 
 ## 3. Penalty term

examples/test_transforms.py

@@ -4,7 +4,7 @@
 from ssqueezepy import Wavelet, TestSignals
 from ssqueezepy.utils import window_resolution
 
-tsigs = TestSignals(N=512)
+tsigs = TestSignals(N=2048)
 #%%# Viz signals #############################################################
 # set `dft` to 'rows' or 'cols' to also plot signals' DFT, along rows or columns
 dft = (None, 'rows', 'cols')[0]
@@ -16,40 +16,67 @@
     ('hchirp', dict(fmin=.2)),
     ('sine:am-cosine', (dict(f=32, phi0=1), dict(amin=.3))),
 ]
-tsigs.demo(signals, N=512)
+tsigs.demo(signals, N=2048)
 
 #%%# 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 = TestSignals(N=2048)
+wavelets = [
+    Wavelet(('gmw', {'beta': 60})),
+    Wavelet(('gmw', {'beta': 5})),
+]
+tsigs.wavcomp(wavelets, signals='all')
 
 #%%#
-tsigs.wavcomp(wavelets, signals=[('#echirp', dict(fmin=.1))], N=512)
+tsigs.wavcomp(wavelets, signals=[('#echirp', dict(fmin=.1))], N=2048)
 
 #%%# Viz CWT vs STFT (& SSQ'd) ###############################################
 # (N, beta, NW): (512, 42.5, 255); (256, 21.5, 255)
-N = 512
+N = 2048
+signals = 'all'
+
 n_fft = N
-win_len = n_fft // 2
+win_len = 720
 tsigs = TestSignals(N=N)
-wavelet = Wavelet(('GMW', {'beta': 21.5}))
+wavelet = Wavelet(('GMW', {'beta': 60}))
 
 NW = win_len//2 - 1
 window = np.abs(sig.windows.dpss(win_len, NW))
-window = np.pad(window, win_len//2)
+window = np.pad(window, (N - len(window))//2)
+assert len(window) == N
 window_name = 'DPSS'
-config_str = '\nNW=%s, win_pad_len=%s' % (NW, len(window) - win_len)
+config_str = 'NW=%s, win_len=%s, win_pad_len=%s' % (
+    NW, win_len, 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=None,
-                  n_fft=n_fft, window_name=window_name, config_str=config_str)
+kw = dict(wavelet=wavelet, window=window, win_len=None, n_fft=n_fft,
+          window_name=window_name, config_str=config_str)
+tsigs.cwt_vs_stft(N=N, signals=signals, **kw)
+
+#%%# Noisy example ###########################################################
+N = 2048
+snr = -2  # in dB
+signals = 'packed-poly'
+
+tsigs = TestSignals(N=N, snr=snr)
+tsigs.cwt_vs_stft(N=N, signals=signals, **kw)
+
+#%%# Ridge extraction ########################################################
+N = 512
+signals = 'poly-cubic'
+snr = None
+n_ridges = 3
+penalty = 25
+
+tsigs = TestSignals(N=N, snr=snr)
+kw = dict(N=N, signals=signals, n_ridges=n_ridges, penalty=penalty)
+tsigs.ridgecomp(transform='cwt',  **kw)
+tsigs.ridgecomp(transform='stft', **kw)

ssqueezepy/README.md

@@ -1,7 +1,12 @@
 # Usage guide
 
 See full README and `examples/` at https://github.com/OverLordGoldDragon/ssqueezepy/ . Also 
-see method/object docstrings via e.g. `help(cwt)`, `help(wavelet)`.
+see method/object docstrings via e.g. `help(cwt)`, `help(wavelet)`. Other READMEs/references:
+
+ - Ridge extraction: https://github.com/OverLordGoldDragon/ssqueezepy/tree/master/examples/ridge_extraction
+ - Generalized Morse Wavelets: https://overlordgolddragon.github.io/generalized-morse-wavelets/
+ - Testing suite: https://overlordgolddragon.github.io/test-signals/
+
 
 ## Selecting `scales` (CWT)
 

ssqueezepy/_cwt.py

@@ -20,7 +20,8 @@ def cwt(x, wavelet='gmw', scales='log-piecewise', fs=None, t=None, nv=32,
 
         wavelet: str / tuple[str, dict] / `wavelets.Wavelet`
             Wavelet sampled in Fourier frequency domain.
-                - str: name of builtin wavelet. See `ssqueezepy.wavs()`/
+                - str: name of builtin wavelet. See `ssqueezepy.wavs()`
+                  or `Wavelet.SUPPORTED`.
                 - tuple: name of builtin wavelet and its configs.
                   E.g. `('morlet', {'mu': 5})`.
                 - `wavelets.Wavelet` instance. Can use for custom wavelet.
@@ -391,13 +392,14 @@ def _process_gmw_wavelet(wavelet, l1_norm):
 def cwt_higher_order(x, wavelet='gmw', order=1, average=None, **kw):
     """Compute `cwt` with GMW wavelets of order 0 to `order`. See `help(cwt)`.
 
-    Yields lower variance and more noise robust representation. See ref[1].
+    Yields lower variance and more noise robust representation. See VI in ref[1].
 
     # Arguments:
         x: np.ndarray
             Input, 1D.
 
         wavelet: str / wavelets.Wavelet
+            CWT wavelet.
 
         order: int / tuple[int] / range
             Order of GMW to use for CWT. If tuple, will compute for each

ssqueezepy/_gmw.py

@@ -2,8 +2,7 @@
 """Generalized Morse Wavelets.
 
 For complete functionality, utility functions have been ported from jLab, and
-largely validated to match jLab's behavior. These can be used to compute
-higher-order wavelets, or pertinent properties thereof. jLab tests not ported.
+largely validated to match jLab's behavior. jLab tests not ported.
 """
 import numpy as np
 from numpy.fft import ifft

ssqueezepy/_stft.py

@@ -43,7 +43,7 @@ def stft(x, window=None, n_fft=None, win_len=None, hop_len=1, fs=None, t=None,
             windowings. Relates to 'overlap' as `overlap = n_fft - hop_len`.
             Must be 1 for invertible synchrosqueezed STFT.
 
-        fs: float
+        fs: float / None
             Sampling frequency of `x`. Defaults to 1, which makes ssq
             frequencies range from 0 to 0.5*fs, i.e. as fraction of reference
             sampling rate up to Nyquist limit. Used to compute `dSx` and