Added time-frequency ridge extraction

[DavidBondesson]

Added folder/files:

ssqueezepy/ridge_extraction.py test/ridge_extract_test.py test/ridge_extract_readme -->imgs test/ridge_extract_readme -->Readme

[OverLordGoldDragon]

Thank you for this significant contribution. Some comments:

1. Wx /= sqrt(scales) no longer needed to L1-norm since v0.5.5+; I've edited accordingly
2. extract_fridges (...great name, but edited to what's more expected) docstring reads "integer value" for penalty, but 0.5 is used in test; I presumed this a typo and removed "integer"
3. ssq_cwt return order will change in v0.6.0, I'll edit later to accommodate
4. Reformatted for Pythonic-ness

I've briefly reviewed reference 2, and confirmed functionality for example signals. For completeness, please add an example where the method fails or works poorly on ssq_cwt (and optionally on cwt). Also let me know if you disagree with some of my edits.

[OverLordGoldDragon]

It'd also help to add as inline comments, where applicable, references to formulae used (e.g. equation numbers from ref 2).

[Rhyst223]

This is a question rather than comment:

I am trying to use SSQ for automatic wave detection in time series, where waves do not necessarily appear for the entirety of the signal. Looking at:

On the extraction of instantaneous frequencies fromridges in time-frequency representations of signals. D. Iatsenko, P. V. E. McClintock, A. Stefanovska. https://arxiv.org/pdf/1310.7276.pdf

I know the signal not appearing for the entirety of the signal poses a problem. However, I have been trying to look at a workaround for this by

• Identifying connected frequency bands in the signal surrounding local maxima in the SSQ amplitude (to constitute a singular event, connectivity determined by scipy.ndimage.label)
• Clip the signal to each of these events
• Perform ridge detection and signal reconstruction

While this originally appeared to be a nice solution initially, the magnitudes of SSQ coefficients can pose some issues. For example, consider this dummy signal containing 4 different frequencies which occur for varying times (can provide code for this if necessary).

With linear colour scaling after SSQ it appears obvious what areas to extract as maxima, however, with log scaling you can see the smaller values are not zero but resemble noise. Within the noise are further "noise local maxima". There is also some frequency leakage for shorter events. Pink crosses in the image below are the true frequencies and the times they occur.

Therefore I need to find a way to mask the noise in the SSQ, but haven't been successful so far. I wondered if you might have any techniques for this? I can isolate the maximal ridges easily enough, but with the pitfall of masking the bands surrounding instantaneous frequencies. I also thought to mask by amplitudes I expect in the wave, but I need the full reconstruction band to compare to this.

An alternative approach would be to isolate events in the abs(CWT) which is smoother, and then map to the SSQ-plane, but since scale -> frequency is not a one-to-one mapping this doesn't seem a mathematically sound approach to me?

If I have posted this in the wrong place or my question is not clear enough please let me know. I can also provide the methods which I have tried so far which fell short if that helps.

[OverLordGoldDragon]

@Rhyst223 Fair question. Please do provide reproducible code, and open a separate Issue as I rather focus this thread on PR progress. Can pastebin, or

<details>
  <summary>Code</summary>

```python
# code here; ignore the "s" below. Keep newline between </summary> & ```
```s

</details>

Code

```python # code here; ignore the "s" below. Keep newline between & ``` ```

[OverLordGoldDragon]

If optimal BW is frequency resolution dependent, we just need to know what the dependence is, and can use wavelets.freq_resolution or utils.window_resolution. I've confirmed 15 to work better than 25 on your test examples; this may be per different wavelet defaults here vs MATLAB's.

Comment on test cases: in test_chirp_lq I found stft_bw=4 to work best, which is quite low; despite 15 working fine for pure sines, it fails with more complicated FM despite it having even greater separability. Mind I have different default configs for wavelet and window, and the two don't necessarily share similar time-frequency localizations. A problem with making BW dependent on frequency resolution is... there's no single "resolution" in CWT, it's literally "multi-resolution", varying with scale. The closest we have to an absolute is a "nondimensional" measure, i.e. adjusted for center frequency.

---

My edits summary

1. `extract_ridges`: support for `transform='cwt'` or `'stft'` 2. `extract_ridges`: support for `get_params` which decides whether `fridges` and `max_energy` are computed and returned (default `False` as I figure only `ridge_idxs` will be of interest most of the time) 3. `extract_ridges`: changed default `BW=25` -> `15` as latter seems to work better with lib's defaults 4. `extract_ridges_test`: added STFT tests, renamed `sig_len, t_vec` -> `N, t` for consistency with other tests & `ssqueezepy` 5. `extract_ridges_test`: added unified function for testing and visualizing `cwt`, `ssq_cwt`, `stft`, `ssq_stft` on a signal 6. An "ideal" sine for CWT purposes is generated with `endpoint=True` _and_ `N=` power of 2 + 1; adjusted accordingly. 7. Docstrings reformat; methods within methods don't need to adhere to full format, better to save lines 8. Added file authorship

Some tests may not run because I've not accounted for all branch differences, for example only mine has t support for stft. Feel free to edit locally or even push, I'll adjust upstream with GMW merger.

[DavidBondesson]

In regards of appropriate BW selection: it makes sense to scale it with the effect of the frequency resolution and the window but at some point I think that a user will still have to select it themselves as it will of course depend heavily on their signal to be evaluated. Thus, I find it sufficient (for the moment at least) to give a suggestion as we have done now through the examples and if anybody has a good idea for how to automate BW- (or other parameter-) selection they can feel free to suggest it.

I adjusted the local settings for my PR but left your edits as is, since I accept them all.

[OverLordGoldDragon]

Have you any info on the approximate relationship? e.g. is BW to grow proportionally with frequency width, or inverse? And we might detach the signal factor in sense that, if the wavelet/window works well on it, so will the BW. Problem is, unsure this holds, as again test_chirp_lq signals had greater separation than pure tones but needed much lesser BW for STFT (and my impression is lesser BW is to handle tighter separations).

[DavidBondesson]

So the concept of BW in this sense i would say can be simple generalized to 'what width does your frequency component have in this time-frequency representation?' Then you simply zero this component after the ridge has been extracted and continue with another run to find the 2nd, 3rd.. strongest component.

To automatically select BW, we must have a way to estimate how wide your signal will appear in your TF representation. If it is wider BW must be increased.
to use your example imagining two signal components with very tight separation, here this algorithm breaks down as the components spill on to each other and there may be cases where you can not separate them appropriately at all. With the case of the test_chirp_lq you would try to find the optimal BW that removes all spill from the 1st components that would be larger (and stable due to the forward backward tracking and penalty factor to frequency jumps) than the 2nd strongest component, while trying not to remove the later partially or all together.

To answer your first question though. I have so far not found any formula for this and my guess is that since Mathworks is also not offering more than letting the users set their required bins to be removed around your main component, it may be a bit trickier than expected.

[OverLordGoldDragon]

So just like in _invert_components we zero the parts we "already visited", discounting them from the operation (here ridge extraction)? That's very helpful, and I figure code is here. So ideally we drop just the right amount from a component so to leave other components intact.

Indeed this cannot be accounted from the wavelet/window alone, at least not entirely. While exact amount to drop varies with "spectral leakage", a wavelet's interaction with said leakage is well-predictable, and I can see tying BW to frequency width working lot better as a default. Unsure I have time to pursue this, but I'll at the least show what to look for; for now I'll just adjust docstrings to better guide on BW.

Thanks for clarifying.

[OverLordGoldDragon]

Here's my take on automating bw selection; it depends on:

1. wavelet/window frequency width (use wavelets.freq_resolution(nondim=False) or utils.window_resolution())
2. number of total rows (cwt scales or stft n_fft) - or more precisely, spacing of rows' center frequencies
3. amount of expected spectral leakage in input. One can manually set an "upper bound" here, but some "measure" is needed (e.g. user specifies tolerance of 0.8 out of 1, we need to know what "1" is)

Illustrated as follows for CWT of sine + cosine:

• Each CWT row is generated via ifft() of product of frequency-domain wavelet (blue, orange) and fft(x_padded) (black)
• The wiggles of black about peaks are the "leakage" that e.g. blue can capture from f1
• Reducing nv increases spacing of scales and thus wavelets, thus fewer CWT rows in total will strongly capture f2
• Reducing frequency width has the same effect, independent of scales spacing
• Example of point 3: brute force generating a variety of cos + sin at f1 & f2 and taking worst case scenario

With STFT the matter's simpler as frequency width is fixed - with CWT, width increases with increasing center frequency (decreasing scale). Which very much asks... why is bw a constant? Makes more sense to vary it linearly from a min to max, like np.linspace(min_bw, max_bw, len(scales)).

[Rhyst223]

I might be able to help somewhat regarding spectral leakage. This week I've been working on some code to identify independent events in the SSQ output, finding the ridges in each event using your function (ran on a slice of abs(SSQ) corresponding to the event), and finally reconstructing the signal with BW estimated from the edges in the SSQ output.

I find the ridges and edges using a Hybrid Hessian filter:

Ng, C. C., Yap, M. H., Costen, N., & Li, B. (2014,). Automatic wrinkle detection using hybrid Hessian filter. In Asian Conference on Computer Vision (pp. 609-622). Springer International Publishing. DOI:10.1007/978-3-319-16811-1_40

Once the ridge extraction function is ran over the period of interest I locate the edge boundaries (lb,ub) above and below the ridge. BW is then estimated as int(max(ub)-min(lb)).

It is not necessarily smooth per se, but it does seem to work. For example, here is a ~5.3e-3Hz wave (Hanning windowed), and lasting for 4 periods:

There is clearly something weird going on beyond the boundaries of the event, but I do reconstruct the signal quite well (the outside boundary bits just have zero amplitude):

Another example for the same frequency wave lasting for 8 periods, for illustration:

Finally another example when there are multiple events going on:

It is far from perfect but thought it might help might regarding BW calculation because in this approach it is fully automated. I hope this is useful, and I wanted to say thank you as your SSQ code has been very helpful in my current PhD project!

I can provide code if necessary for reproducibility.

[OverLordGoldDragon]

@Rhyst223 Glad my work's been helpful, and thanks for your input.

So the idea is, we draw bands about ridges? If we can do that, why not just take the midpoint and call it a ridge? Seems this method should have its own bw-like parameter, else it could act a stabler alternative for ridge extraction in some cases.

Whatever the case I see this as a worthwhile contribution, so you are welcome to PR if you get the chance. Unfortunately I lack time to dig into this but will minimally validate the algorithm and assist in library integration (Pythonicness etc).

[DavidBondesson]

Since I and colleagues are certainly interesting in frequency tracking I would welcome further contribution from @Rhyst223 's suggested method as well! At the very least this may yield a nice way to offer automated BW selection for the current algorithm. Particular thing to consider for ridge extraction for our purpose would be to track signal components over the whole time which is why I would be very interested to see how the hybrid hessian filter algorithm performs on our standard examples.

@OverLordGoldDragon I added one more reference to specific equation in 'ridge_extraction.py' and would say that I feel sufficiently satisfied that this can be merged if you have no other comments.

[OverLordGoldDragon]

Thanks again for your contribution, hope the library serves well your endeavors.

[OverLordGoldDragon]

An idea on auto bw: for a given wavelet and nv, the number of rows for a pure tone is fixed regardless of frequency (I was wrong earlier; while wavelet resolution changes, log scaling 'offsets' it). One can thus take CWT of a pure tone, and based on a set cutoff threshold, find distance from max to threshold, and difference in number of rows:

Boundary effects naturally complicate things. Neglecting them and assuming no leakage, one can compute bw for a given wavelet and nv analytically or with a simple black box algorithm. More will be known on handling leakage and boundaries with these done.

One could repeat this cutoff logic to an expected F.M., like linear chirp, or take average of that and a pure tone. Many possibilities.

[Rhyst223]

I like this idea as a starting point, and it looks like it should be easy to get bw as a function of wave duration for however long your signal is.

I have a few questions: 1) What are the expected boundary effects (possibly analytically), and does it change with the frequency of the wave?

2) You have used 20% of the peak height to threshold and get the automated bw, was this just an arbitrary cvalue that seemed to do the trick?

and a comment: I'd love to talk about/figure out how to get bw from (wavelet, nv, spectral leakage) analytically!

I have almost finished writing up my code using the Hybrid-Hessian filter so I should be able to contribute that soon. The main issue I'm finding with my approach is that is seems to be overestimating when an event occurs. But, I plan to add a noise term so that if there is something you know about the waves you seek (ie above a certain amplitude), you can remove times you don't want.

[OverLordGoldDragon]

1. I'm not well-informed on ascertaining boundary effects analytically but have come across publications that do (no longer recall where though). What I do know, the lower the frequency, the worse the boundary effects; think in time domain convolution, where a large scale wavelet is nonzero over a larger portion of input, whereas for highest frequencies it's just a few samples. Without providing exact math, I could answer more specific questions, particularly on remedies with padding.
2. Totally arbitrary.

Look forward to the contribution, though am currently occupied with adding GPU & CPU acceleration to ssqueezepy, and may take up on Time-frequency scattering (I recommend looking into for machine learning classification applications).