Subtraction of PSD Topomaps

💬 Support & Discussions
, , ,
1

Hi,

I would like to make a plot of two subtracted psd topomaps. Is this possible?
I have epoched my data and would like to subtract different epochs from each other. Kind of like this:
epochs[‘Rec1’].plot_psd_topomap() - epochs[‘Rec2’].plot_psd_topomap()

Thanks!

Ps. I found a "Maling List Archive (Read-only) topic about this. Sorry for the duplicate, but the old one has no answers.

2

Sure, this is possible, but you would have to operate on arrays. First you would use psd_welch function on both you epochs objects (you could also use multitapering if you prefer):

psd1, freq = psd_welch(epochs1, ...)
psd2, freq = psd_welch(epochs2, ...)

(the arguments that would go where the ... are depend on your data and desired frequency resolution etc. see the psd_welch docs)
To subtract the power you can just do:

psd_diff = psd2 - psd1

And then you can plot it using mne.viz.plot_topomap:

# find the index of the frequency you want to plot
frequency = 10
freq_idx = np.argmin(np.abs(freq - frequency))

# plot the desired frequency
mne.viz.plot_topomap(psd_diff[:, freq_idx], epochs.info)
1 Like
3

Great thanks! It works :slight_smile:
Can it also work on frequency bands? I would like to look the alpha and gamma band?

4

Sorry, I can also try myself! :wink:

Is it ok to do like this:

# find the index of the frequency you want to plot
frequency_band = [8, 12]
freq_idx = [np.argmin(np.abs(freq - frequency_band[0])), np.argmax(np.abs(freq - frequency_band[1]))]

# plot the desired frequency
mne.viz.plot_topomap(psd_diff[:, freq_idx].mean(axis = 1), epochs.info)
5

Oh, sorry, I didn’t see your reply, yes it is almost the same solution as mine below.
Be sure to specify a range (i use it with slice) instead of just the two indices.

Sure, instead of one frequency index you could use a range, for example:

freq_limits = [8, 12]
freq_idx = [np.argmin(np.abs(freq - frequency)) for frequency in freq_limits]
freq_range = slice(freq_idx[0], freq_idx[1] + 1)

# average the frequency range:
psd_avg = psd[:, freq_range].mean(axis=1)

# plot
mne.viz.plot_topomap(psd_avg, epochs.info)

(BTW - for gamma you may prefer to use multitapering with some frequency smoothing instead of welch)

2 Likes