Added folder/files:

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

Regarding:

Licensing; unsure how to proceed here; original's says to "Redistributions in binary form must reproduce the above copyright notice" - but I'm not "redistributing" it, I'm distributing my rewriting of it

The original code is published under a 2-clause BSD, so basically you are free to chose your license. It would be great if you can pick one, because code without a license on GitHub is a gray area (impossible to use / contribute) -- ideally something similar like BSD/MIT/Apache :wink:

Hi, it's me again

Wanted to know if there is some example of maprange around so I can do some quick tests and study it from there ?

On my tests, I have observed that if a sudden peak of energy while on a window, is much higher than the already "detected" clusters of energy, the first clusters detected of lower energy are marginalized, in the sense that they become dark ( they lose the color they had )

What is a practical way, if there is any, to tune the scales so this effect does not happen or at least it's less obvious ?

Thank you very much for this wonderful library

I have no idea if this is expected behavior in a SST, but at least to me it was a surprise. I noticed this at first when analyzing raw audio material, which revealed oscillations in overtones much more than what my hear would tell. But it is also possible to reproduce the effect synthetically:

The following example shows a mix of sine waves with frequencies i * f_base. In each plot one more frequency is added; the left side shows the transform of the mix, the right side is the transform of only the last/added harmonic (just as a check how its transform looks in isolation):

import numpy as np
import matplotlib.pyplot as plt

from ssqueezepy import synsq_cwt_fwd

SAMPLE_RATE = 44100


def sine(freq, num_samples, amp=1.0):
    ts = np.arange(num_samples) / SAMPLE_RATE
    wave = amp * np.sin(2.0 * np.pi * freq * ts)
    return wave


def mix(*waves):
    return np.array(waves).sum(axis=0) / len(waves)


def plot(Tx, ax):
    xs = np.arange(Tx.shape[1] + 1)
    ys = np.arange(Tx.shape[0] + 1)
    img = ax.pcolormesh(xs, ys, np.abs(Tx))
    plt.colorbar(img, ax=ax)


def plot_overtone_series(num_overtones, freq=440.0, num_samples=SAMPLE_RATE // 4):
    waves = [
        sine(freq * i, num_samples)
        for i in range(1, num_overtones + 1)
    ]

    Tx_mix, *_ = synsq_cwt_fwd(mix(*waves), fs=1.0, nv=32)
    Tx_add, *_ = synsq_cwt_fwd(waves[-1], fs=1.0, nv=32)

    fig, axes = plt.subplots(1, 2, figsize=(18, 8))
    plot(Tx_mix, axes[0])
    plot(Tx_add, axes[1])
    fig.tight_layout()
    plt.savefig("/tmp/overtones_{}".format(num_overtones))
    plt.show()


if __name__ == "__main__":
    for i in range(1, 9):
        plot_overtone_series(i)

Observations:

• Even though the sine waves have a constant instantaneous frequency, the SST plot shows relatively large oscillations -- each time in the highest harmonic. For instance here is a zoomed in version of the n = 4 case:

  

  If I understand it correctly there are 32 "voices" per octave. The oscillation covers ~10 bins, which in musical terms would correspond to more than 3 semitones -- which seems a lot for a constant frequency signal.

• The highest harmonic seems to cover up the harmonics in between. Perhaps this oscillation in the highest harmonic is actually an reassignment artifact of "stealing" energy from the lower harmonics?

In the Matlab cwt documentation, it gives this useful example (emphasis added by me):

[wt,f] = cwt(___,fs) specifies the sampling frequency, fs, in Hz as a positive scalar. cwt uses fs to determine the scale-to-frequency conversions and returns the frequencies f in Hz. If you do not specify a sampling frequency, cwt returns f in cycles per sample. If the input x is complex, the scale-to-frequency conversions apply to both pages of wt. If x is a timetable, you cannot specify fs. fs is determined from the RowTimes of the timetable.

In this example, if you provide the sampling frequency, then the frequencies are returned. It would be really useful if ssqueezepy could do this too (or at least provide a scale2frequency function to do the required conversion with the specified wavelet).

Hi, for my project https://github.com/falseywinchnet/streamcleaner I have previously been using the librosa stft. However, after i was able to make the ssqueezepy stft behave similar to it by appending 128 samples and then slicing all but the last 127 from the istft output, I switched over to ssqueezepy, believing that for my purposes https://github.com/librosa/librosa/issues/1279 this was an important caveat I should consider.

However, librosa's stft only consumes ~2% CPU, and the ssqueezepy on a similar workload takes up 25%. I'm concerned this is due to my crude attempt at making the two behave similar, when doing this the STFT representation is MSE very close to librosas.

Is there a proper way to make the stft of ssqueezepy generate a 257x376 complex representation from 48,000 samples and return 48,000 samples and perform with similar compute requirements?

I have used ssq_cwt() to get Tx on 0.6.0 version of ssqueezepy. I would like to use speed improvents of 0.6.1 and get the same Tx of 0.6.0 but using 0.6.1 version. Is it possible? How can it be done?

Thanks in advance

EDIT: I have reformulated the question

Currently the default for stft and ssq_stft is to return all scales up to the Nyquist frequency, but it would be helpful if there were a simple way to specify a maximum frequency or, better yet, a specific set of scales that should be queried. https://github.com/OverLordGoldDragon/ssqueezepy/blob/master/ssqueezepy/_ssq_stft.py#L241 Particularly in cases where you don't want to query everything up to the Nyquist frequency for efficiency purposes, it would be helpful to have this option.

Dear Eugene Brevdo and Gaurav Thakur,

First of all, thank you very much for your python library ssqueezepy. it is very useful and is the first library that really implements very well the icwt in python (at least from my concern). I would like to know if It would be possible to improve the ssq_cwt method applying the ConceFT (Daubechies et al., 2016) algorithm. What is your opinion? is it worth it?

If that were the case, what wavelet families could be applied? it could be applied an stack of synchrosqueezed scalograms generated by morlet families at different number of cycles or on the contrary morse wavelet have a clear advantage? For the STFT, what is the advantage of the orthonormal Hermite functions over the Discrete Prolate Spheroidal (Slepian) Sequences?

Thank you very much in advance for your time and for your help

Best Regards, Roberto Cabieces

Hello,

You may be aware of this, and I apologize if this is redundant or if I'm using the code improperly. I upgraded from ssqueezepy 0.6.0 to 0.6.1 today but am running into an issue with ssq_stft input parameter ssq_freqs. I am using a numpy ndarray of frequencies of interest, but in 0.6.1 I get the following error:

File "/home/user/ssqtest.py", line 90, in SSQ_Tvst, Sxo, ssq_freqs, scaleS = ssq.ssq_stft(xo, ssq_freqs=fcm, fs=10, t=time)

File "/home/user/anaconda3/envs/base/lib/python3.7/site-packages/ssqueezepy/_ssq_stft.py", line 112, in ssq_stft maprange='maximal', transform='stft')

File "/home/user/anaconda3/envs/base/lib/python3.7/site-packages/ssqueezepy/ssqueezing.py", line 204, in ssqueeze _ssqueeze(Tx, w, Wx, dWx, *args)

File "/home/user/anaconda3/envs/base/lib/python3.7/site-packages/ssqueezepy/ssqueezing.py", line 143, in _ssqueeze gamma, out=Tx, Sfs=Sfs)

File "/home/user/anaconda3/envs/base/lib/python3.7/site-packages/ssqueezepy/algos.py", line 149, in ssqueeze_fast fn(*args, **params)

TypeError: missing argument 'vmin'

I did some poking around in algos.py and to me it looks like **params has a vlmin and dvl, rather than the vmin and dv required for the _ssq_stft function. However, if I simply change the input parameter names to match (aka vmin --> vlmin and dv --> dvl), I get the same error as above with the added errors:

File "/home/user/anaconda3/envs/base/lib/python3.7/site-packages/numba/core/dispatcher.py", line 420, in _compile_for_args error_rewrite(e, 'typing')

File "/home/user/anaconda3/envs/base/lib/python3.7/site-packages/numba/core/dispatcher.py", line 361, in error_rewrite raise e.with_traceback(None)

TypingError: Cannot determine Numba type of <class 'numba.core.ir.UndefinedType'>

I wasn't sure what to do with this error. Thanks!

So far the package isn't really pip installable due to lack of a setup.py. I thought I'd add something basic to get you started, including properly specifying dependencies.

I guess few people in the Python community actually download releases from GitHub. So I guess you can safe you some time, and simply add a pip + git based installation instruction to the README like:

pip install git+https://github.com/OverLordGoldDragon/ssqueezepy

I could also help with publishing to PyPI in the long term.

Hi, Thanks for your wonderful work.

Currently, only NumPy arrays are supported as input. Is there a way to use torch tensors already on the GPU?

If not, if you give me some pointers, I can work on implementing that feature.

Best regards Nabarun

Hi,

I've encountered what i believe to be a bug in a call to ssqueezepy.ssqueeze to synchro-squeeze a continous wavelet transform. I'll try to explain the issue as best i can, because i can't send the exact implementation. This is also my first issue, to bear with me if i get anything wrong in the format or conventions :)

I'm running ssqueezepy v0.6.3

The call is of the form:

Tx, fvec = sq.ssqueeze(
			Mx, w=phase, ssq_freqs="log", wavelet=wavelet, scales="log", fs=sampleRate, maprange="maximal",
			transform="cwt"
		)

where,

wavelet = sq.Wavelet("morlet") is set before calling the method. and is of type(wavelet) = ssqueezepy.wavelets.Wavelet

and Mx, phase are the normal matrices of the transform and phase produced by ssqueezepy.cwt()

In 'squeezing.py/squeeze' after the call to _process_args is completed, the following if statement on line 168 calls the method process_scales if transform == 'cwt'. See below

    if transform == 'cwt':
        scales, cwt_scaletype, _, nv = process_scales(scales, N, get_params=True)

However the call to process_scales does NOT pass the wavelet as an input argument (which it expects and sets to None if not passed). This results in the exception ValueError("must set `wavelet` if `scales` isn't array") being raised on line 220 in cwt_utils.pt/_process_args as called from process_scales.

The argument wavelet passed to ssqueeze() was in fact set correctly.

Hopefully that makes sense...

Thanks, J.

Suppose that we have a highly non-sinusoidal wave that is composed of multiple waves. Assuming that we know the frequency range that we are interested in, how should we define cc and cw in issq_cwt so the reconstructed signal only contains the components in he frequency range of interest. An illustrative example will be helpful.

The following example returns wc = -3.1354 but according to the specs it should be nonnegative, unless I'm misunderstanding something.

from ssqueezepy import wavelets, Wavelet
wavelet = Wavelet(("morlet", {"mu":14.0}))
_ = wavelets.center_frequency(
    wavelet,  
    scale=np.array([1.0]),  
    kind='peak',
    N=1024,
    viz=1
)

Thread for sharing/discussing applications on real-world data. Feel free to open a separate Issue if needed.

Comments occasionally cleared to keep focus on examples.