overlordgolddragon / ssqueezepy Goto Github PK

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

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.

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 😉

I run the minimum example code with and without multi-threaded executions.

• os.environ['SSQ_PARALLEL'] = '0'
• os.environ['SSQ_PARALLEL'] = '1'

Results

1. It works and with multi-threaded is faster.
2. However, I saw small results with and without GPU. Attached is the image. Is this normal?
3. Also, I repeated the process and I have observed changes both in computation time (although the GPU is faster) and different plots.

I have tried to run the documented example of ridge extraction.

However, I encountered an index error everytime I use extract_ridges.

Hi, it's me again.
Seems like in your 0.6.3 release, you switched by anadvertance, while rewriting that function, ridge_f and ridge_e, lines 138 and 139.

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.

Thread for various implementation or theoretical notes / ideas.

Discussed in https://github.com/falseywinchnet/streamcleaner/discussions/16

• drop the padding step on each frame and only process 95 (191) time segments
  I dont know why the padding is done but i assume it has something to do with centering frequency components properly. It is only done on the edge of each frame, but if we switch to processing STFT and appending it to a buffer, we have to change this behavior.
  @OverLordGoldDragon why is the padding done on the stft?

You may be aware that there has been recent research effort to combine synchrosqueezing with Empirical mode decomposition.
However, EMD has some shortcomings and it would seem there is a 21st century decomposition technique called intrinsic time-scale decomposition " Intrinsic time-scale decomposition: time−frequency−energy analysis and real-time filtering of non-stationary signals" Mark G Frei and Ivan Osorio, doi: 10.1098/rspa.2006.1761

The original publication from 2007 is paywalled, but accessible via sci-hub. It features some insightful representation of the work that makes it evident that ITD is a fundamentally new and powerful weapon, as great as wavelets or FFT, which is missing from our arsenal. I don't know if it's a reversible decomposition method however it seems to be described as such.

I believe this screenshot sums it up pretty nicely.
From page four we have:
"4. The ITD method overcomes essentially all the limitations of EMD described earlier. In particular_, the ITD only involves O(n) computations, operates in realtime,_ and although it can be applied to windows of data, it naturally operates on inter-extrema segments of data, with a time-delay equal to the inter-extrema interval, which is determined by the frequency content of the signal being decomposed. When applied to a window of data, the unavoidable edge effects that stem from uncertainty regarding wave morphology in the intervals before the first and after the last signal extrema do not propagate inward as they can with the EMD, but rather are temporally constrained to margins defined by the first and last two extrema in the window. The ITD does not use splines and thus avoids the ‘difficulties’ they introduce. By construction, the ITD method guarantees proper rotation components2 and precisely preserves temporal information regarding TFE information and all the signal critical points. Finally, and perhaps most importantly, the sifting procedure introduced by Huang et al. (1998) in an attempt to produce proper rotation components with the EMD, but which also causes smearing and smoothing of TFE information, is not needed in the ITD. However, for heuristic purposes, we investigated its application within the ITD framework and determined that, as with the EMD, it has the general effect of smoothing variations in amplitude between adjacent extrema. Exhaustive sifting within a window of data results in a component that is purely frequency modulated with constant amplitude. Unlike the EMD, after each ITD sifting iteration, the extrema points remain identical to those of the input signal."

The best explanation of the core algorithm i could find, and i do not know if this is is correct, i only know that the algorithm isn't presented all at once elsewhere, is found on https://arxiv.org/pdf/1404.3827v1.pdf on page 26.

There is also a matlab implementation found at https://viewer.mathworks.com/?viewer=plain_code&url=https%3A%2F%2Fwww.mathworks.com%2Fmatlabcentral%2Fmlc-downloads%2Fdownloads%2Fbacae6a7-59d2-4ad9-89b8-f4d5dd60f422%2F3d85f570-571f-4e45-8049-b0e10d205469%2Ffiles%2Fitd.m&embed=web

From " The De-noising Algorithm Based on Intrinsic Time-scale Decomposition" doi:10.4028/www.scientific.net/AMR.422.347 jinxia zeng, goufu wang, et al, in 2008, (also paywalled but available on scihub), we find they show evidence supporting that ITD offers similar computational complexity to wavelet transforms but improved noise reduction similar to EMD, while being orders of magnitude easier to compute, as claimed by frei and osorio. Unlike wavelets, the issue of a wavelet base is not encountered making this method superior to both.

I would genuinely like to see it accessible in python and feel this is the best place to announce discussion about it. I have sought help from some engineers from around the world to look into adapting it, perhaps for this repository and directly for inclusion in scipy and numpy libraries. I feel this method will really transform neural network inputs. I have also heard about, but not seen research publications concerning, the application of synchrosqueezing to EMD, so it may even be possible to go further and apply it to ITD, resulting in a truly powerful and transcendent method.

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

A thread to share/discuss example comparisons between implementations.

NOTE: I'm no longer actively developing features for ssqueezepy. Requests can still be submitted, here if short, or as Issues if detailed. If I don't reply within a month, it's probably not happening. I may react with below emoji to confirm.

Contributions as Pull Requests are still welcome, but should mostly complete themselves.

Thread for tracking proposals for extending functionality.

1 checkmark = in progress; 2 = done

Since your synchrosqueezing wavelet transform algorithm takes discrete time sequences as input, why is it called the Continuous Wavelet Tranform and not the Discrete Time Wavelet Transform?

Could you point to some of the equations you used in handling discrete time inputs?

Description

• 'reflect' vs 'zero'
• mention article & relink 1-int inv

Code for reproduction

No response

I am working on binary systems, for this, we are analyzing light curves and some of them are periodic. We have this data through observations with different telescopes. The data we obtain are Time, Flux, and error_flux. And I would like to do an analysis with your ssqueezepy code, to analyze the light curves through a WAVALET transform. how do I define my time, and flux to be able to run the code?

When I run the make_scales function I got an error: 'UnboundLocalError: local variable 'i' referenced before assignment'

the error occurs in find_downsampling_scale(wavelet, scales):

File /usr/local/lib/python3.9/dist-packages/ssqueezepy/utils/cwt_utils.py:580, in find_downsampling_scale(wavelet, scales, span, tol, method, nonzero_th, nonzero_tol, N, viz, viz_last) 575 print(("Failing scale: (idx, scale) = ({}, {:.2f})\n" 576 "out of max: (idx, scale) = ({}, {:.2f})" 577 ).format(i, float(scales[i]), len(scales) - 1, float(scales[-1]))) 578 _viz(psihs, psihs_peaks, joint_peak, psihs_nonzeros, i) --> 580 return i if (i < n_groups - 1) else None

Thus bug occurs when the signal len: 30,000 | low cut: 150 | high_cut: 250 | scaletype: log_piecewise | wavelet: morelet.

#81 (comment)

Most likely you are aware of this. Just in case it should be tracked explicitly.

Out of curiosity I gave synsq_cwt_inv a try, and noticed that currently it doesn't work:

/home/fabian/git/ssqueezepy/ssqueezepy/synsq_cwt.py in synsq_cwt_inv(Tx, fs, opts, Cs, freqband)
    166     for n in range(Cs.shape[1]):
    167         TxMask = np.zeros(Tx.shape)
--> 168         UpperCs = min(max(Cs[:, n] + freqband[:, n], 1), len(fs))
    169         LowerCs = min(max(Cs[:, n] - freqband[:, n], 1), len(fs))
    170 

TypeError: 'int' object is not subscriptable

This makes sense because freqband is initialized as scalar here, so the indexing here will fail.

trying to run with this file :

Channels: 2
Sample width: 2
Frame rate (sample rate): 48000
Frame width: 4
Length (ms): 115136
Frame count: 5526528.0

it's a 22MB wav file

first I load the file with librosa :

#load file with librosa
x, sr = librosa.load(file)

Tx, Wx, ssq_freqs, scales, *_ = ssq_cwt(x, nv=16, wavelet='morlet', astensor=True, cache_wavelet=True)

RuntimeError: CUDA out of memory. Tried to allocate 2.53 GiB (GPU 0; 7.79 GiB total capacity; 5.76 GiB already allocated; 947.44 MiB free; 5.76 GiB reserved in total by PyTorch) If reserved memory is >> allocated memory try setting max_split_size_mb to avoid fragmentation. See documentation for Memory Management and PYTORCH_CUDA_ALLOC_CONF

I checked nvidia-smi and yes it starts taking memory until it reaches the mem limit, then it releases when it crashes

is there any trick to analyze such files without hitting the memory limits ? Thanks

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.

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

Description

Hi. Thanks for the library, is being really handy!

I found a way to crash the icwt function in the current version - 0.6.3

Provided are the time series and code required to reproduce the issue

from ssqueezepy import ssq_cwt, icwt
import numpy as np

data = np.array([ 41135817341,  29083562061,  26175547452,  16588370958,
        17264085441,  31947336829,  40770974039,  30242059107,
        21692004719,  29867476527,  27246574439,  34163220274,
        68204556440,  50913575242,  54912007015,  31183975654,
        27132421514,  42009436760,  35329942625,  30818458597,
        28970212744,  28574793478,  26188097173,  24957784918,
        18372538715,  18027170497,  20965695707,  21381535161,
        23552740328,  26267239923,  30767551159,  18100418740,
        16390821947,  21594638208,  26715546990,  24598943708,
        25814972520,  49899834488,  29641127858,  28688807249,
        24150249025,  25810220018,  33042430345,  31158743333,
        25905575359,  24302954056,  22927802083,  39974475562,
        48765202697,  42932549127,  33631012204,  31421555646,
        24021799169,  23565495303,  35574561406,  28624673855,
        31758955233,  40212386158,  35887249746,  28148218301,
        23553591896,  25849159141,  28389250717,  26288169966,
        25120229769,  28881249043,  15978259885,  15886817043,
        28575544847,  23555719219,  32837431722,  37127036580,
        27265804688,  22987346289,  22994133555,  35123501685,
        27753685646,  30931623076,  23747613147,  40509610260,
        27595671000,  23102307723,  31666498758,  31878280659,
        31962253368,  31028679593,  42326789564,  30116729776,
        24366810591,  32637854078,  34483360283,  33225232872,
        30182031010,  29123998928,  23613051457,  25245861652,
        28813460025,  43403978910,  35239757134,  32194477850,
        48469528171,  36913738894,  34493951963,  50212088965,
        51091116622,  37872380889,  36389011503,  30123362273,
        24957448100,  31254779144,  40177002624,  36791346508,
        46363793975,  41135767926,  38896078052,  26149643168,
        23359966112,  44148798321,  58571439619,  53071298734,
        41037843771,  43975248085,  18719537670,  20765955327,
        30484729489,  35887278685,  33223790572,  34711412966,
        29227315390,  16437423167,  16837262532,  27425022774,
        28711532910,  24950173846,  44219840004,  38452356727,
        16192235532,  17988916650,  27472552998,  30580012344,
        22425387184,  24493974420,  32459287866,  16104440957,
        22128794335,  30202235805,  47761524910,  58895950537,
        49625110402,  43994715910,  40369840645,  31486345556,
        45668466815,  39819303159,  55552169483,  43228750179,
        64072727950,  37846047609,  35082693210,  53510852236,
       118992465607, 102905151606,  83202283721,  55871616488,
        29717699419,  27209183682,  49630243054,  36599436183,
        33925512989,  27868914022,  26862218609,  16106223492,
        21313378652,  37429485518,  30726828760,  32958875628,
        26129037414,  18678255976,  18000008764,  20443898509,
        27743025156,  23581685468,  29523576583,  22895392882,
        19539705127,  16217776704,  16824520830,  22209086834,
        19889922369,  19675404389,  20496603770,  20328426366,
        12706781969,  14122486832,  19617581341,  26634741631,
        25534481470,  20964448341,  24031608960,  14463581825,
        10924354698,  17221074814,  22722096615,  14882945045,
        16441573050,  15329265213,   9744636213,  11656379938,
        11886957804,  15748580239,  17005713920,  14472237479,
        15929162910,  11239186456,   9244361700,  12097775227,
        13903079207,  18421743322,  13692758566,  14413662913,
         7714767174,   9768827914,  18624736866,  15808338949,
        18372283782,  34971338710,  29225029694,  38967784639,
        19298407543,  26792494050,  24999983362,  30005625418,
        21152848261,  28799154319,  32442278429,  24746386230,
        26518700512,  26405069715,  30685366709,  26357839322,
        25383335641,  14712928379,  27423687259,  27205595568,
        22837828665,  26683255504,  32066936882,  27083066007,
        15639298538,  19564262605,  23825006542,  27187964471,
        25371367758,  32572572185,  27078406594,  16356226232,
        17821046406,  23918742607,  26792596581,  32483312909,
        39316664596,  41358451255,  19625427158,  25555105670,
        28987376573,  31252098714,  30199996781,  30476264066,
        26811744928,  16100721565,  16644534842,  22660763494,
        20535363434,  24662841200,  20386398516,  26062404610,
        11166012913,  13488941056])


wavelet = 'morlet'

_, cwt, _, scales = ssq_cwt(data, wavelet=wavelet)

icwt(cwt, wavelet=wavelet, scales=scales)

I am able to get the examples shown in the front page to work fine, however, all attempts towards ridge extraction fail on my system. The frustrating part is that I don't even get a traceback, only total kernel failure. I've tried through the command line and through a jupyter notebook (I realize they are linking to the same kernel). Any ideas on how I might rectify this situation?

I've attempted to run the examples on the readme for ridge extraciton and observed the same result when running the core example. On the latter I can't proceed past line 63 (https://raw.githubusercontent.com/OverLordGoldDragon/ssqueezepy/master/examples/extracting_ridges.py).

Synchrosqueezing is an intriguing concept and am interested to try it out on our data.

System Profile:
OSX 10.13.6 (High Sierra)
Python 3.8 (Anaconda, Fresh Install)

Install/Dependency Trace:
Downloading ssqueezepy-0.6.1-py3-none-any.whl (121 kB)
|████████████████████████████████| 121 kB 1.2 MB/s
Requirement already satisfied: numba in /opt/anaconda3/lib/python3.8/site-packages (from ssqueezepy) (0.53.1)
Requirement already satisfied: scipy in /opt/anaconda3/lib/python3.8/site-packages (from ssqueezepy) (1.6.2)
Requirement already satisfied: matplotlib in /opt/anaconda3/lib/python3.8/site-packages (from ssqueezepy) (3.3.4)
Requirement already satisfied: numpy in /opt/anaconda3/lib/python3.8/site-packages (from ssqueezepy) (1.20.1)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /opt/anaconda3/lib/python3.8/site-packages (from matplotlib->ssqueezepy) (2.4.7)
Requirement already satisfied: pillow>=6.2.0 in /opt/anaconda3/lib/python3.8/site-packages (from matplotlib->ssqueezepy) (8.2.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /opt/anaconda3/lib/python3.8/site-packages (from matplotlib->ssqueezepy) (1.3.1)
Requirement already satisfied: cycler>=0.10 in /opt/anaconda3/lib/python3.8/site-packages (from matplotlib->ssqueezepy) (0.10.0)
Requirement already satisfied: python-dateutil>=2.1 in /opt/anaconda3/lib/python3.8/site-packages (from matplotlib->ssqueezepy) (2.8.1)
Requirement already satisfied: six in /opt/anaconda3/lib/python3.8/site-packages (from cycler>=0.10->matplotlib->ssqueezepy) (1.15.0)
Requirement already satisfied: llvmlite<0.37,>=0.36.0rc1 in /opt/anaconda3/lib/python3.8/site-packages (from numba->ssqueezepy) (0.36.0)
Requirement already satisfied: setuptools in /opt/anaconda3/lib/python3.8/site-packages (from numba->ssqueezepy) (52.0.0.post20210125)
Installing collected packages: ssqueezepy
Successfully installed ssqueezepy-0.6.1

_cwt freq is of shape (N,1)
_ssq_cwt freq is of shape (N)
_ssq_cwt scales is of shape (N,1)

Everything should be of the same shape

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
)

I read a recent publication on fastCWT which may help speed up the computation. I posted a question and link here: https://dsp.stackexchange.com/q/83469/62730

Hi,

Is there any straightforward way to extract the averaged wavelet power of a given frequency band using ssqueezepy ? For instance :

I have a 0.1s EEG recording (sf = 3500), and I want the mean power of the entire recording regarding alpha/beta (10-30Hz).

Thanks a lot,

Best,

Hi @OverLordGoldDragon

I haven't gotten a chance to look at this in much detail yet, but here are some initial overall comments regarding the code.

General comments

• 1.) This toolbox reimplements its own CWT. Is it possible to reuse the one already in PyWavelets?
  I know there are some numerical limitations to the CWT as currently implemented in PyWavelets (e.g. PyWavelets/pywt#535, PyWavelets/pywt#531). I don't yet know if these are specifically a problem for this application.
  Alternatively, would the SciPy implementation work?

• 2.) PyWavelets uses the Numpydoc for docstrings as described at:https://docs.scipy.org/doc/numpy-1.15.0/docs/howto_document.html

I wouldn't make refactoring the docstrings a first priority though (The same applies for PEP8 issues)

• 3.) It would be good to make the argument handling more Pythonic.
  For example, rather than specify opts=None or opts={} and then parse it with a _get_opts helper, I would explicitly make the signature of synsq_cwt_fwd

def synsq_cwt_fwd(x, t=None, fs=None, nv=32, type='morlet', difftype='direct', gamma=None)```

so that it is easy for the user to determine the names and defaults of all options.

The same applies to the various other functions such as synsq_squeeze that are called within. It is better to explicitly list the arguments by name than hiding them within an opts dictionary.

I'm new new CWT, so this may be a silly question, but I'm unsure why issq_cwt can simply sum the reals of the output of ssq_cwt, without considering the frequencies. In this Matlab page, the formula given at the bottom is nearly identical to the formula in issq_cwt, except it includes division by a (which I believe is the frequency).

Is this a bug in issq_cwt? If not, how should I reconstruct my original signal if I have modified the frequency bins? If it helps, I used the following parameters when calling ssq_cwt: maprange=(16.35, 22100), ssq_freqs="log", fs=44100.

Is there an environ setting to enable which GPU to use and the ability to use more than one for CUPY and Torch?

I am on a machine with 4 GPUs but it only uses GPU #1 which is in use by other processes.

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.

Setting parameters (vectorized=False, cache_wavelet=True) to ssq_cwt gives "WARNING: cache_wavelet=True requires vectorized=True; setting to False."

It looks like in the _cwt.py _process_args vectorized is undefined?

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

1. I'm using ssqueezepy to calculate signal's instantaneous frequency, The method like below:

Tx, Wx, ssq_freqs, scales, *_ = ssq_cwt(data, ('gmw', {'beta': 2, 'gamma': 3}), order=0, padtype='wrap', fs=fs)
ridge_idxs = extract_ridges(Tx, scales, penalty=0.5)
frequencies = ssq_freqs[::-1][ridge_idxs[:, 0]]  # this is the result

2. The result like below:
   

3. I think this method can get the correct result, but the result have some steps which not smooth.So what cause this?Which parameters should change?

Thanks!

tmax,N, f1, f2 = 1, 1000, 20,50
fs = N/tmax
t = np.linspace(0, tmax, N, endpoint=True)
x1,x2 = np.sin(2*np.pi *f1 * t), np.cos(2*np.pi * f2 * t)
x = x1 + x2 

Twx, Wx, ssq_freqs,scales = ssq_cwt(x,fs=fs,scales='linear',maprange='peak',nv=100,wavelet='morlet')
ridge_idxs,fridge,max_energy = extract_ridges(Twx, scales, penalty=20, n_ridges=2, bw=2,get_params=True,transform='cwt')
plt.plot(ssq_freqs[::-1][ridge_idxs])
plt.show()

We have two signals, with frequency 20 and 50, no noise. extract_ridges manages to get the right frequencies, but the jumps between both of them are very surprising. Going up to penalty=400 still keeps those jums.

Is the code for EEG use case available?

Great work with this library!

I am attempting to generate scaleogram images as shown by imshow() , in order to classify them with a convolutional neural network. Of course, I am able to generate the 2D array that is passed to imshow() with cwt(data, 'morlet') , where data is my 1D signal. However, I am hoping to derive the 3D array with the colormap applied as it is in your imshow() function, so that the resulting 3D array can be passed into a neural network.

Any suggestions? Thanks!

My understanding is that the ssq_cwt function optionally accepts a set of wavelet scales rather than corresponding frequencies. While I understand that there are a few different interpretations of "corresponding frequencies" for wavelet scales, I've noted that there is a simple, analytical relationship (see Equation A3 of "An Introduction to Wavelet Analysis in Oceanography and Meteorology", also referenced in "A Practical Guide to Wavelet Analysis") between wavelet scale and frequency for certain wavelet functions. Would it be difficult to add functionality for the user to pass in frequences rather than wavelet scales and have those converted to wavelet scales on the backend?

Amazing package! Many thanks!

As a side note: the SE explanation that this uses a prior to improve on the traditional time / frequency uncertainty principle is really what meant most to me; I am used to signal processing and Kalman filters, so I know how priors help, and I was aware of the uncertainty principle, so when I read this SE explanation I understood at once that i) this is amazing ii) how this works "deep down", iii) this is exactly what I need. I think I am not alone being in this situation :) .

Wondering if this could be stated at the start of the readme, in big font, something in the kind of the extension to your catchline: "Synchrosqueezing is a powerful reassignment method that focuses time-frequency representations, and allows extraction of instantaneous amplitudes and frequencies. It uses a prior (that the signal follows the modulation approximation) to improve on the frequency / time uncertainty principle".

Installed using the following command:
pip install git+https://github.com/OverLordGoldDragon/ssqueezepy
and
from ssqueezepy.experimental import phase_ssqueeze
gives me this error:

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
/tmp/ipykernel_25228/4269281225.py in <module>
     18 from ssqueezepy import ssq_cwt, ssq_stft, TestSignals, Wavelet
     19 from ssqueezepy.visuals import imshow
---> 20 from ssqueezepy.experimental import phase_ssqueeze

ImportError: cannot import name 'phase_ssqueeze' from 'ssqueezepy.experimental' (/home/lucy/Downloads/general/lib/python3.9/site-packages/ssqueezepy/experimental.py)

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 librosa/librosa#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?

Thread for tracking ideas in need of explaining, and found explanations.

1 checkmark = explanation found; 2 = documented

Hi, I wanted to ask if the CWT can be used on a discrete time-sampled digital signal

This simulator samples incoming messages, realtime every 1ms. The result is an array of integers ( amount of messages per ms ) that ranges from 0 to N like this [0,12,4,26,32,2,5] and so on.

I wanted to know if I need to do something to the data first or I can just use it like this with the CWT

Thanks !

The doc for extract_ridges

# Returns
        ridge_idxs: np.ndarray
            Indices for maximum frequency ridge(s).
        fridge: np.ndarray
            Frequencies tracking maximum frequency ridge(s).
        max_energy: np.ndarray [n_timeshifts x n_ridges]
            Energy maxima vectors along time axis.

Though in the code, fridge are the scales. (line 137). Please clarify the doc (or change the doc), since fridge has the same dimension as scales. It feels like one of them should change.

Also it feels weird that the following code does not give the expected scales

Twx, Wx, ssq_freqs,scales = ssq_cwt(x,fs=fs,scales='log',maprange='peak',nv=32,wavelet='morlet')
ridge_idxs,fridge,max_energy = extract_ridges(Wx, scales, penalty=10, n_ridges=2, bw=5,get_params=True,transform='cwt')
np.allclose(np.exp(fridge),np.squeeze(scales[ridge_idxs])) # Gives True, shall not be np.exp

In the first line, scales is logarithmic. Using transform='cwt', the outputed scales from extract_ridges are the log of ssq_cwt scales, which are already in log. There's some confusion there.

Third thing :
It feels like scales and ssq_freqs arrays are reversed, and to retrieve the right frequencies, I should use the following snipet

Twx, Wx, ssq_freqs,scales = ssq_cwt(x,fs=fs,scales='log',maprange='peak',nv=64,wavelet='morlet')
ridge_idxs,fridge,max_energy = extract_ridges(Twx,ssq_freqs[::-1], penalty=10, n_ridges=2, bw=8,get_params=True,transform='linear')

With the use of the [::-1] to reverse ssq_freqs array., which originates from lines 229 and 230 of ssqueezing.py. Tell me if that is the intended functionnality.

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).

Is there a way / built in function to easily convert between scales and frequencies, given sampling rate?

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?

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!

As promised (in a StackExchange thread) ... some examples for you to use for testing (and comparison) and as examples and to add to your test suite or examples archive. Just grab the sounds from the videos.

15 second segment of a drum solo at 1/4 speed
https://www.youtube.com/watch?v=6orozX1GD1w

(House music with preacher talking)
https://www.youtube.com/watch?v=nd2J4xTrSHQ

(The drum at full and 1/4 speed, a segment of the 2001 theme and the House segment with preacher at full and 1/4 speed)
https://www.youtube.com/watch?v=itUSUau6DJM

((The Vampire Song)
https://www.youtube.com/watch?v=OugT7uGGtNg

You may want to consider expanding your software to include a fuller range of transforms (or to use for comparison). The time-scale transform, which the videos show an early implementation of (in the process of being upgraded), is described in more detail in the description page of the Vampire Song. Its main property is that it is invariant with respect to arbitrary frequency reassignment (e.g. even a random scrambling of frequencies): an arbitrarily-reassigned scalogram will yield, upon inversion, the same result as the original scalogram - both reproducing the original signal.

Originally posted by @RockBrentwood in #9 (comment)