overlordgolddragon / ssqueezepy Goto Github PK
View Code? Open in Web Editor NEWSynchrosqueezing, wavelet transforms, and time-frequency analysis in Python
License: MIT License
Synchrosqueezing, wavelet transforms, and time-frequency analysis in Python
License: MIT License
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)
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.
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)
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
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?
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
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.
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.
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).
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
)
A thread to share/discuss example comparisons between implementations.
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.
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".
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
Status | Description | Code / Resources |
1. [] [] |
padsignal : more methods in PyWavelets; pull?
|
|
2. [] [] |
padsignal : Interpolation, Fourier series extension
|
|
3. [] [] |
padsignal : Boundary wavelets are described as superior to zero-padding, periodic, and symmetric extensions
|
|
4. [] [] | Generalized Morse wavelets "capture the essential idea of the Morlet, while avoiding its deficiencies ... may be recommended as an ideal starting point for general purpose use". | |
5. [x][x] | Fast wavelets -- Implement wavelet recomputation in C or wrap `@numba.jit` | |
6. [] [] | Cachelets -- Allow wavelets to be cached (on disk) for reuse; provide database of default-parameter wavelets computed for several common powers of 2 |
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
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)
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?
Thread for tracking ideas in need of explaining, and found explanations.
1 checkmark = explanation found; 2 = documented
Status | Description | Findings / Resources |
1. [ ] [x] |
cwt_fwd : comparison against PyWavelets, Scipy
|
|
2. [ ] [x] |
padsignal , boundary effects
|
|
3. [ ] [ ] |
cwt_fwd coeff computation method
|
|
4. [ ] [ ] |
synsq_cwt_fwd Oscillations in overtones
|
I run the minimum example code with and without multi-threaded executions.
os.environ['SSQ_PARALLEL'] = '0'
os.environ['SSQ_PARALLEL'] = '1'
Results
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.
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,
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,
Is there a way / built in function to easily convert between scales and frequencies, given sampling rate?
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?
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?
'reflect'
vs 'zero'
No response
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.
_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
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.
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
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!
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
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!
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?
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
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 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?
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 !
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 ๐
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.
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.
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.
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.
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
Error message:
ValueError Traceback (most recent call last)
Input In [2], in <cell line: 80>()
76 wavelet = 'morlet'
78 _, cwt, _, scales = ssq_cwt(data, wavelet=wavelet)
---> 80 icwt(cwt, wavelet=wavelet, scales=scales)
File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/_cwt.py:391, in icwt(Wx, wavelet, scales, nv, one_int, x_len, x_mean, padtype, rpadded, l1_norm)
389 idx = logscale_transition_idx(scales)
390 x = icwt(Wx[:idx], scales=scales[:idx], **kw)
--> 391 x += icwt(Wx[idx:], scales=scales[idx:], **kw)
392 return x
393 ##########################################################################
394
395 #### Invert ##############################################################
File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/_cwt.py:378, in icwt(Wx, wavelet, scales, nv, one_int, x_len, x_mean, padtype, rpadded, l1_norm)
376 wavelet = Wavelet._init_if_not_isinstance(wavelet)
377 # will override nv
to match scales
's
--> 378 scales, scaletype, _, nv = process_scales(scales, x_len, wavelet, nv=nv,
379 get_params=True)
380 assert (len(scales) == na), "%s != %s" % (len(scales), na)
382 #### Handle piecewise scales case ########################################
383 # nv
must be left unspecified so it's inferred automatically from scales
384 # in process_scales
for each piecewise case
File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:247, in process_scales(scales, N, wavelet, nv, get_params, use_padded_N)
244 nv = int(nv)
245 return scaletype, nv, preset
--> 247 scaletype, nv, preset = _process_args(scales, nv, wavelet)
248 if isinstance(scales, (np.ndarray, torch.Tensor)):
249 scales = scales.reshape(-1, 1)
File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:228, in process_scales.._process_args(scales, nv, wavelet)
225 if scales.squeeze().ndim != 1:
226 raise ValueError("scales
, if array, must be 1D "
227 "(got shape %s)" % str(scales.shape))
--> 228 scaletype, _nv = infer_scaletype(scales)
229 if scaletype == 'log':
230 if nv is not None and _nv != nv:
File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:289, in infer_scaletype(scales)
286 # ceil to avoid faulty float-int roundoffs
287 nv = int(np.round(1 / np.diff(np.log2(scales), axis=0)[0]))
--> 289 elif logscale_transition_idx(scales) is None:
290 raise ValueError("could not infer scaletype
from scales
; "
291 "scales
array must be linear or exponential. "
292 "(got diff(scales)=%s..." % np.diff(scales, axis=0)[:4])
294 else:
File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:380, in logscale_transition_idx(scales)
378 scales = asnumpy(scales)
379 scales_diff2 = np.abs(np.diff(np.log(scales), 2, axis=0))
--> 380 idx = np.argmax(scales_diff2) + 2
381 diff2_max = scales_diff2.max()
382 # every other value must be zero, assert it is so
File <array_function internals>:180, in argmax(*args, **kwargs)
File /opt/conda/envs/gst/lib/python3.8/site-packages/numpy/core/fromnumeric.py:1216, in argmax(a, axis, out, keepdims)
1129 """
1130 Returns the indices of the maximum values along an axis.
1131
(...)
1213 (2, 1, 4)
1214 """
1215 kwds = {'keepdims': keepdims} if keepdims is not np._NoValue else {}
-> 1216 return _wrapfunc(a, 'argmax', axis=axis, out=out, **kwds)
File /opt/conda/envs/gst/lib/python3.8/site-packages/numpy/core/fromnumeric.py:57, in _wrapfunc(obj, method, *args, **kwds)
54 return _wrapit(obj, method, *args, **kwds)
56 try:
---> 57 return bound(*args, **kwds)
58 except TypeError:
59 # A TypeError occurs if the object does have such a method in its
60 # class, but its signature is not identical to that of NumPy's. This
(...)
64 # Call _wrapit from within the except clause to ensure a potential
65 # exception has a traceback chain.
66 return _wrapit(obj, method, *args, **kwds)
ValueError: attempt to get argmax of an empty sequence
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)
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.
Is the code for EEG use case available?
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
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!
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.