OverLordGoldDragon / ssqueezepy

Synchrosqueezing, wavelet transforms, and time-frequency analysis in Python
MIT License
618 stars 96 forks source link

Penatly on extract_ridges #51

Closed Bulbille closed 2 years ago

Bulbille commented 2 years ago
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()

Capture

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.

OverLordGoldDragon commented 2 years ago

Ridge extraction was implemented by @DavidBondesson, whose input is welcome.

I suspect the problem is scales='linear', which strictly belongs neither in tranform='cwt' nor 'stft', but 'stft' appears to be a much better fit:

extract_ridges(, transform='stft', n_ridges=2)

I also strongly discourage scales='linear' in the general case; if linear frequency tiling is desired, STFT is expressly designed for it.

DavidBondesson commented 2 years ago

So I am fairly certain that this is indeed due to the linear scaling of your transform and I noticed the behaviour as well during implementation. I did play around a bit with inputting the log of the transforms amplitude (add a small factor to avoid log(0)) which often seemed to generate more expected results. However, I gave that up when I was uncertain if this would actually translate the underlying function when we are trying to account for the distance between frequency bins.

In short I would just agree with @OverLordGoldDragon, but if there is a need for me to adjust implementation, delve deeper into the point or add comments then do let me know. However, due to my current work I will need side projects officially approved by my company which means that I can't just jump in right away.

Bulbille commented 2 years ago

Indeed, my mistake about not checking linear vs cwt/stft. I was toying around because

Twx, Wx, ssq_freqs,scales = ssq_cwt(x,fs=fs,scales='log',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')

gives this weird output for the first and last 20% of the data, so I was playing with other things and forgot I was in linear.

Thanks for your fast reply !