nelpy / ghost

Spectral analysis tools for neuroscience data
MIT License
10 stars 1 forks source link

Adaptive sampling of wavelets #4

Open jchutrue opened 5 years ago

jchutrue commented 5 years ago

At the lower frequencies, we need more samples to get an accurate representation of the wavelet, but at higher frequencies, we need fewer. I need to figure out a better way of sampling the wavelets rather than the default of 16000 samples or so.

ckemere commented 5 years ago

I'm confused. I always define the root wavelet (i.e., scale = 1), which corresponds to something close to the highest frequency being sampled. Lower frequencies are just stretched versions of this. The duration of the wavelet should scale with increasing scale (decreasing frequency)?

jchutrue commented 5 years ago

Hmm I was doing it backward, where I computed the number of points that the lowest frequency wavelet needed to be sampled sufficiently. I set this level to be 1. Then for example, if the frequency doubled, the number of samples should be halved.

For the example I was testing, the lowest frequency was around 0.15 Hz and needed 25000 points. At 150 Hz, my strategy said I only need 25, but when I checked the wavelet's representation, it was definitely undersampled

ckemere commented 5 years ago

This is the benefit of the "scales"-focused approach. In your example, you might have been undersampling 0.15Hz without realizing.

jchutrue commented 5 years ago

You are correct. The wavelet implementation I adapted my code around dealt only with frequencies, not scales. I'll take another look at the Morse papers and see what they did.

jchutrue commented 5 years ago

There's a formula here that is analogous to the temporal bandwidth: https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2016.0776

See equation 2.8. From preliminary tests, about 4 times this value seemed to give sufficient sampling. I'll go verify on the lowest frequencies too.

ckemere commented 5 years ago

The whole point of the wavelet is that you convolve with h( (t -a)/b) ). b is the scale. So if the wavelet h(t) is time limited, it should be time limited the same amount regardless of scale. Equation 2.8 is just saying this. I don't understand why you don't just use scale internally???

jchutrue commented 5 years ago

From a maintenance perspective, I find frequencies easier to work with, and this is how Lilly's own Morse wavelet package is implemented. In fact, all of the functions I've adapted use normalized radians, not scales.

Internally for the transform object, I do scale the number of wavelet samples by the ratio of the frequencies to the mother wavelet peak frequencies, per equation 2.8.

This adaptive wavelet sampling based on the scale/frequency is now up on the develop and unstable branches, and is working. You can also specify a voices_per_octave argument like MATLAB's cwt function.