fastlib / fCWT

The fast Continuous Wavelet Transform (fCWT) is a library for fast calculation of CWT.
Apache License 2.0
263 stars 53 forks source link

Morlet wavelets with fixed temporal window length #38

Open snazau opened 1 year ago

snazau commented 1 year ago

Hello! Thank you for the implementation!

I'm working with EEG data and want to replace method from mne library with your implementation to speed up my computations

With mne I perform CWT with following code:

freqs = np.arange(1, 40.01, 0.1)
power_spectrum = mne.time_frequency.tfr_array_morlet(
        samples,  # samples.shape = [number_of_samples, number_of_eeg_channels, temporal_dim]
        sfreq=128,
        freqs=freqs,
        n_cycles=freqs,
        output='power',
        n_jobs=-1
)

I tried to implement the exact same behaviour with fCWT but failed:

morl = fcwt.Morlet(2.0)  # The problem is here I suppose

scales = fcwt.Scales(
        morl,
        fcwt.FCWT_LINFREQS,
        fs=int(raw_data.info['sfreq']),
        f0=freqs[0],
        f1=freqs[-1],
       fn=len(freqs),
)

nthreads = 4
use_optimization_plan = False
use_normalization = False
fcwt_obj = fcwt.FCWT(morl, nthreads, use_optimization_plan, use_normalization)

signal = samples[0, 0]  # for example
output = np.zeros((len(freqs), signal.size), dtype=np.complex64)
fcwt_obj.cwt(
        signal,
        scales,
        output,
)
power = np.abs(output) ** 2

Could you give some recommendations on how to choose the right sigma parameter of fcwt.Morlet to create Morlet wavelets with temporal window with fixed length? Is it possible to get the same behaviour as mne with the current implementation of fCWT?

fastlib commented 1 year ago

Hi Snazau,

Thank you for working on an extension of fCWT to MNE. I think this is very usefull! Could you specify a bit more by what you mean with Morlet Wavelets having a temporal window of a fixed length? I guess fCWT is currently not supporting it, but it should only mean creating a custom wavelet to get it to work. However, what should be the constraints of this wavelet? Maybe we can brainstorm together on this topic.