Event related Epochs

Support & Discussions
1
  • MNE 0.23.0 - Python version: 3.8.8
  • operating system: Windows 10

Hello,

I’m working with a MEG project. When creating epochs from stim_channel related events, I get an output of over 500 events for the selected channel. I would like to create fixed length epochs for all those events, and from it, plot many time frequency/power spectrum density graphs using Welch’s method. I want to be able to do an average and most importantly, an individual evaluation of each of those epochs graphs, but didn’t find a way to plot it like this. Could someone help me?

Thank you!

2

I don’t think we have implemented any spectro-temporal analyses with Welch’s method. For time-frequency we have stockwell, multitaper and morlet. Any of those functions will work with either Epochs data or Evoked data, so you can do the averaging either before or after the TFR step.

Unfortunately right now EpochsTFR doesn’t have a plot method, but AverageTFR does… so a sort of hackish work-around for the plotting would be to do a for-loop with epochs.iter_evoked(), in which you convert to TFR using one of the functions mentioned above (should yield an AverageTFR object), and then visualize it with its .plot() method.

3

Thank you for the quick reply and great suggestion. I’ll try it !

4

n = 1
for x in epochs.iter_evoked():
if n>= 1 and n<= 10:
mne.time_frequency.tfr_morlet(x, freqs, n_cycles, use_fft=False, return_itc=False, decim=3, n_jobs=1, picks=None, zero_mean=True, average=True, output=‘power’, verbose=None)
n= n + 1

I’m using the lines above to plot a selected number of epochs from the 500 I have. I defined the variables “freqs” and “n_cycles”, however, I’m getting the following message:

"ValueError: At least one of the wavelets is longer than the signal. Use a longer signal or shorter wavelets."

I’ve tried to change the value of all variables but I still get this error. Any suggestions on how to fix it? Thank you.

5

Hi Breno,

as the error says, you are asking for a too long wavelet for the length of your signal. The wavelets method estimates the time-frequency content by sliding this wavelet along the time axis of your data and estimating the fit (roughly speaking). Now, that of course will not work if your wavelet is longer than the data.
This could for example happen if you have 500 ms of data, but ask for 10 cycles of 10 Hz (= 1000 ms).

Hope this helps
Britta

1 Like
6

Unfortunately, still having trouble.

Based on the epochs above, I’m trying to plot a power/frequency graph using an average of a combination of specific(group) MEG channels and a combination of epochs at the same time but I can’t find a way to combine the channels, just the epochs. Here is the code I have been using to do it:

kwargs = dict(fmin=2, fmax=40, n_jobs=1)
psds_welch_mean, freqs_mean = psd_welch(epochs[55:56], average=‘mean’, **kwargs)
psds_welch_median, freqs_median = psd_welch(epochs[55:56], average=‘median’, **kwargs)

psds_welch_mean = 10 * np.log10(psds_welch_mean)
psds_welch_median = 10 * np.log10(psds_welch_median)

channel_indices = mne.pick_types(info, meg=‘mag’, eeg=True)

meg_channels = mne.pick_types(info, meg=True)[:10]
channel_types = [mne.io.pick.channel_type(info, ch) for ch in meg_channels]

ch_name = ‘MEG0121’
ch_idx = epochs.info[‘ch_names’].index(ch_name)
epo_idx = 0

_, ax = plt.subplots()
ax.plot(freqs_mean, psds_welch_mean[epo_idx, ch_idx, :], color=‘k’,
ls=’-’, label=‘mean of segments’)
ax.plot(freqs_median, psds_welch_median[epo_idx, ch_idx, :], color=‘k’,
ls=’–’, label=‘median of segments’)

ax.set(title=‘Welch PSD ({}, Epoch {})’.format(ch_name, epo_idx),
xlabel=‘Frequency (Hz)’, ylabel=‘Power Spectral Density (dB)’)
ax.legend(loc=‘upper right’)
plt.show()

=====================================
Using the code above, I can plot an average of the epochs I defined but can’t find a way to insert more channels into it. This is the code I used to try to insert more channels, but when I plot it, I just get a blank graph.

ch_name = ‘MEG0121’
ch_idx = epochs.info[‘ch_names’].index(ch_name)
epo_idx = 0

b = np.hstack([psds_welch_mean[epo_idx, 14, :], psds_welch_mean[epo_idx, 17, :]])
c = np.hstack([psds_welch_median[epo_idx, 14, :], psds_welch_median[epo_idx, 17, :]])
_, ax = plt.subplots()
ax.plot(freqs_mean, b, color=‘k’,
ls=’-’, label=‘mean of segments’)
ax.plot(freqs_median, c, color=‘k’,
ls=’–’, label=‘median of segments’)

ax.set(title=‘Welch PSD ({}, Epoch {})’.format(ch_name, epo_idx),
xlabel=‘Frequency (Hz)’, ylabel=‘Power Spectral Density (dB)’)
ax.legend(loc=‘upper right’)
plt.show()

14 and 17 are MEG channel indices.
I would like to plot a graph like this one for exemple:
Screenshot_5

Does anyone knows how I can combine those channels + epochs and get a mean?
@alexrockhill

7

Have you seen this tutorial Frequency and time-frequency sensor analysis — MNE 0.23.0 documentation? Specifically the part on psd_welch does pretty much exactly what you are asking.

If you want to take the mean across all epochs you can do

psd_grand_mean = psds_welch_mean.mean(axis=0)

or

psd_grand_median = np.median(psds_welch_median, axis=0)
8

Yes, I got the code above exactly from there, however, I can’t find a way on where/how I would insert more MEG channels to the average. On the line “ch_name =”, if I add more channels it doesn’t work. Do you know how I can do it?

9

You could do

for ch_idx, ch_name in enumerate(epochs.info['ch_names']):
    _, ax = plt.subplots()
    ax.plot(freqs_median, np.median(psds_welch_median, axis=0)[ch_idx], color='k',
            ls='--', label='median of segments')
    ax.set(title='Welch PSD ({})'.format(ch_name),
           xlabel='Frequency (Hz)', ylabel='Power Spectral Density (dB)')
    plt.show()
10

I think you understood me wrong. I need a single plot for the selected channels + epochs. I get an error everytime I select more than 1 channel for the plot

11

If you want to plot more than one channel on the same plot, you can do

_, ax = plt.subplots()
for ch_idx, ch_name in enumerate(epochs.info['ch_names']):
    ax.plot(freqs_median, np.median(psds_welch_median, axis=0)[ch_idx], color='k',
            ls='--', label=ch_name)
    ax.set(title='Welch PSD ({})'.format(ch_name),
           xlabel='Frequency (Hz)', ylabel='Power Spectral Density (dB)')

plt.show()
12

Sorry, but with this code, I will plot all channels in the same plot with many lines (one for each channel), however, what I need to plot is a single line which will be the mean power of all channels I select

14

Solved it. I used “epochs[:].plot_psd” and defined the channels I wanted using picks=[]

Thanks everyone