How is the Stockwell Transform frequency resolution determined?

💬 Support & Discussions
1
  • MNE 0.24.0
  • Mac 11.6.1

I’m trying to understand how frequency bins are determined for the stockwell transform of an array. As a simple example, I’m computing it for a 5-second chirp signal with sampling frequency 1000 hz, which starts at 50 hz and ends at 400 hz, increasing in a linear manner.

My code so far:

 import numpy as np
 from scipy.signal import chirp
 from mne.time_frequency import tfr_array_stockwell as stockwell

 fs=1000
 L= 5000 #signal length in ms.
 t=np.arange(0,L/1000,1/fs)
 w=chirp(t, f0=5, f1=400, t1=5, method='linear') # 5000 samples in signal

 w1=np.expand_dims(w,axis=(0,1)) # change to 1 epoch -by -1 channel -by- 5000
 fmin = 0  
 fmax = 25  

 S, itc, freqs=stockwell(w1,fs,fmin=fmin,fmax=fmax)

In this case freqs contains frequencies that start at fmin=0 and end at fmax=25 with 205 bins.

It’s not clear to me how the number of bins is determined using the values chosen above. Is there an easy way to see this without diving into all the source code for the stockwell command?

2

see n_fft parameter. it’s what determines the delta frequency (frequency resolution)

import numpy as np
from scipy.signal import chirp
from mne.time_frequency import tfr_array_stockwell as stockwell

fs = 1000
L = 2**12  # signal length in ms.
n_fft = L
t = np.arange(0, L) / fs
w = chirp(t, f0=5, f1=400, t1=5, method="linear")  # 5000 samples in signal

# w1 = np.expand_dims(w, axis=(0, 1))  # change to 1 epoch -by -1 channel -by- 5000
w1 = w[np.newaxis, np.newaxis, :]
fmin = 0
fmax = 25

S, itc, freqs = stockwell(w1, fs, fmin=fmin, fmax=fmax, n_fft=n_fft)

delta_freq = np.mean(np.diff(freqs))
print(delta_freq)
print(fs / n_fft)

HTH
Alex

1 Like
3

Yes–this is very helpful. Thanks!