diff --git a/.github/workflows/conda.yml b/.github/workflows/conda.yml index 3db1c1c..d426ab2 100644 --- a/.github/workflows/conda.yml +++ b/.github/workflows/conda.yml @@ -46,9 +46,9 @@ jobs: conda run python3 -m pip install --upgrade pip conda run python3 -m pip install -r requirements.txt conda run python3 -m pip install -r requirements-dev.txt - + - name: Set up ssqueezepy - run: conda run python3 setup.py develop + run: conda run pip install -e . - name: Test shell: bash -l {0} diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ef5cee..876b52e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -31,7 +31,10 @@ - `nan_checks` now defaults to `True` only for NumPy inputs - CI changes: conda -> micromamba - docstring references improved + - docs: added missing references (w.r.t. original toolbox) - README section added + - README examplesp improved + - updated github workflows code ### 0.6.3 (1-23-2022): QoL, cleanups, fixes diff --git a/README.md b/README.md index 615f6a7..901ab4e 100644 --- a/README.md +++ b/README.md @@ -66,9 +66,9 @@ See [Performance guide](https://github.com/OverLordGoldDragon/ssqueezepy/blob/ma -### 3. Testing suite: CWT vs STFT, reflect-added parallel linear chirp +### 3. Testing suite: CWT vs STFT, reflect-added parallel A.M. linear chirp - + ### 4. Ridge extraction: cubic polynom. F.M. + pure tone; noiseless & 1.69dB SNR @@ -78,7 +78,7 @@ See [Performance guide](https://github.com/OverLordGoldDragon/ssqueezepy/blob/ma ### 5. Testing suite: GMW vs Morlet, reflect-added hyperbolic chirp (extreme time-loc.) - + ### 6. Higher-order GMW CWT, reflect-added parallel linear chirp, 3.06dB SNR diff --git a/ssqueezepy/_cwt.py b/ssqueezepy/_cwt.py index 5f9414e..eed8b05 100644 --- a/ssqueezepy/_cwt.py +++ b/ssqueezepy/_cwt.py @@ -13,8 +13,7 @@ def cwt(x, wavelet='gmw', scales='log-piecewise', fs=None, t=None, nv=32, l1_norm=True, derivative=False, padtype='reflect', rpadded=False, vectorized=True, astensor=True, cache_wavelet=None, order=0, average=None, nan_checks=None, patience=0): - """Continuous Wavelet Transform, discretized, as described in - Sec. 4.3.3 of [1] and Sec. IIIA of [2]. Uses FFT convolution via frequency- + """Continuous Wavelet Transform. Uses FFT convolution via frequency- domain wavelets matching (padded) input's length. Uses `Wavelet.dtype` precision. @@ -72,8 +71,8 @@ def cwt(x, wavelet='gmw', scales='log-piecewise', fs=None, t=None, nv=32, l1_norm: bool (default True) Whether to L1-normalize the CWT, which yields a more representative - distribution of energies and component amplitudes than L2 (see [3], - [6]). If False (default True), uses L2 norm. + distribution of energies and component amplitudes than L2 (see [3]). + If False (default True), uses L2 norm. derivative: bool (default False) Whether to compute and return `dWx`. Requires `fs` or `t`. @@ -107,6 +106,9 @@ def cwt(x, wavelet='gmw', scales='log-piecewise', fs=None, t=None, nv=32, > 0 computes `cwt` with higher-order GMWs. If tuple, computes `cwt` at each specified order. See `help(_cwt.cwt_higher_order)`. + NOTE: implementation may be not entirely correct. Specifically, + alignment by center frequency rather than scales may be optimal. + average: bool / None Only used for tuple `order`; see `help(_cwt.cwt_higher_order)`. @@ -129,29 +131,36 @@ def cwt(x, wavelet='gmw', scales='log-piecewise', fs=None, t=None, nv=32, differentiation (effectively, derivative of trigonometric interpolation; see [4]). Implements as described in Sec IIIB of [2]. - # References: - 1. Wavelet Tour of Signal Processing, 3rd ed. S. Mallat. - https://www.di.ens.fr/~mallat/papiers/WaveletTourChap1-2-3.pdf + # Note: + CWT is cross-correlation of wavelets with input. For zero-phase wavelets + (real-valued in Fourier), this is equivalent to convolution. All + ssqueezepy wavelets are zero-phase. If a custom general wavelet is + used, it must be conjugated in frequency, and it should *not* be used + with synchrosqueezing (see one-integral inverse References in `icwt`). - 2. The Synchrosqueezing algorithm for time-varying spectral analysis: + # References: + 1. The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications. G. Thakur, E. Brevdo, N.-S. Fučkar, and H.-T. Wu. https://arxiv.org/abs/1105.0010 - 3. How to validate a wavelet filterbank (CWT)? John Muradeli. + 2. How to validate a wavelet filterbank (CWT)? John Muradeli. https://dsp.stackexchange.com/a/86069/50076 + 3. Wavelet "center frequency" explanation? Relation to CWT scales? + John Muradeli. + https://dsp.stackexchange.com/a/76371/50076 + 4. The Exponential Accuracy of Fourier and Chebyshev Differencing Methods. E. Tadmor. http://webhome.auburn.edu/~jzl0097/teaching/math_8970/Tadmor_86.pdf - 5. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. + 5. Wavelet Tour of Signal Processing, 3rd ed. S. Mallat. + https://www.di.ens.fr/~mallat/papiers/WaveletTourChap1-2-3.pdf + + 6. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ cwt_fw.m - - 6. Rectification of the Bias in the Wavelet Power Spectrum. - Y. Liu, X. S. Liang, R. H. Weisberg. - http://ocg6.marine.usf.edu/~liu/Papers/Liu_etal_2007_JAOT_wavelet.pdf """ def _vectorized(xh, scales, wavelet, derivative, cache_wavelet): if cache_wavelet: @@ -382,7 +391,7 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True, 7. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ - synsq_cwt_fw.m + synsq_cwt_iw.m """ #### Prepare for inversion ############################################### na, n = Wx.shape diff --git a/ssqueezepy/_ssq_cwt.py b/ssqueezepy/_ssq_cwt.py index e0964ef..a1bfd14 100644 --- a/ssqueezepy/_ssq_cwt.py +++ b/ssqueezepy/_ssq_cwt.py @@ -180,9 +180,8 @@ def ssq_cwt(x, wavelet='gmw', scales='log-piecewise', nv=None, fs=None, t=None, Decomposition. I. Daubechies, J. Lu, H.T. Wu. https://arxiv.org/pdf/0912.2437.pdf - 4. Synchrosqueezing-based Recovery of Instantaneous Frequency from - Nonuniform Samples. G. Thakur and H.-T. Wu. - https://arxiv.org/abs/1006.2533 + 4. Synchrosqueezed Wavelet Transform Explanation. John Muradeli. + https://dsp.stackexchange.com/a/71399/50076 5. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ diff --git a/ssqueezepy/_ssq_stft.py b/ssqueezepy/_ssq_stft.py index bda9aca..05b6fc9 100644 --- a/ssqueezepy/_ssq_stft.py +++ b/ssqueezepy/_ssq_stft.py @@ -70,6 +70,10 @@ def ssq_stft(x, window=None, n_fft=None, win_len=None, hop_len=1, fs=None, t=Non 1. Synchrosqueezing-based Recovery of Instantaneous Frequency from Nonuniform Samples. G. Thakur and H.-T. Wu. https://arxiv.org/abs/1006.2533 + + 2. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. + https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ + synsq_stft_fw.m """ if x.ndim == 2 and get_w: raise NotImplementedError("`get_w=True` unsupported with batched input.") @@ -157,6 +161,10 @@ def issq_stft(Tx, window=None, cc=None, cw=None, n_fft=None, win_len=None, 2. Fourier synchrosqueezed transform MATLAB docs. https://www.mathworks.com/help/signal/ref/fsst.html + + 3. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. + https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ + synsq_stft_iw.m """ def _process_args(Tx, window, cc, cw, win_len, hop_len, n_fft, modulated): if not modulated: @@ -224,6 +232,10 @@ def phase_stft(Sx, dSx, Sfs, gamma=None, parallel=None): robustness properties and new paleoclimate applications. G. Thakur, E. Brevdo, N.-S. Fučkar, and H.-T. Wu. https://arxiv.org/abs/1105.0010 + + 3. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. + https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ + phase_stft.m """ S.warn_if_tensor_and_par(Sx, parallel) if gamma is None: diff --git a/ssqueezepy/_stft.py b/ssqueezepy/_stft.py index 2f040fa..63728a0 100644 --- a/ssqueezepy/_stft.py +++ b/ssqueezepy/_stft.py @@ -65,7 +65,7 @@ def stft(x, window=None, n_fft=None, win_len=None, hop_len=1, fs=None, t=None, Whether to use "modified" variant as in [1], which centers DFT cisoids at the window for each shift `u`. `False` will not invert once synchrosqueezed. - Recommended to use `True`; see "Modulation" below. + Recommended `True`. See "Modulation" and [2] below. derivative: bool (default False) Whether to compute and return `dSx`. Uses `fs`. @@ -93,6 +93,8 @@ def stft(x, window=None, n_fft=None, win_len=None, hop_len=1, fs=None, t=None, of high frequency, making inversion and synchrosqueezing very unstable. Details & visuals: https://dsp.stackexchange.com/a/72590/50076 + Better explanation in ref [2]. + # Returns: Sx: [(n_fft//2 + 1) x n_hops] np.ndarray STFT of `x`. Positive frequencies only (+dc), via `rfft`. @@ -110,6 +112,17 @@ def stft(x, window=None, n_fft=None, win_len=None, hop_len=1, fs=None, t=None, 1. Synchrosqueezing-based Recovery of Instantaneous Frequency from Nonuniform Samples. G. Thakur and H.-T. Wu. https://arxiv.org/abs/1006.2533 + + 2. Equivalence between "windowed Fourier transform" and STFT as + convolutions/filtering. John Muradeli. + https://dsp.stackexchange.com/a/86938/50076 + + 3. STFT: why overlapping the window? John Muradeli. + https://dsp.stackexchange.com/a/88124/50076 + + 4. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. + https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ + stft_fw.m """ def _stft(xp, window, diff_window, n_fft, hop_len, fs, modulated, derivative): Sx = buffer(xp, n_fft, n_fft - hop_len, modulated) @@ -207,6 +220,10 @@ def istft(Sx, window=None, n_fft=None, win_len=None, hop_len=1, N=None, 2. Invertibility of overlap-add processing. B. Sharpe. https://gauss256.github.io/blog/cola.html + + 3. Synchrosqueezing Toolbox, (C) 2014--present. E. Brevdo, G. Thakur. + https://github.com/ebrevdo/synchrosqueezing/blob/master/synchrosqueezing/ + stft_iw.m """ ### process args ##################################### n_fft = n_fft or (Sx.shape[0] - 1) * 2