compute_psd (welch) - Basic questions

đź’¬ Support & Discussions
1
  • MNE version: e.g. 0.24.0
  • operating system: e.g. Windows 10

Hi,
first of all many thanks for your great package and community. Here a common scenario with some basic questions I’m sure that will help many users that are exploring this package and the EEG signal.

I’m running some test on a 10 seconds EEG signal.
I would like to inspect the power in different frequency bands, related to each epoch.
For example, I may compare the theta power of the epoch (from 0 to 1 sec) with the theta power of the epoch (from 1 to 2 sec)

So, thanks to your function “mne.make_fixed_length_events” I created 10 events equally distributed (then I assigned them different ID), in order to create 10 different epochs (the first from 0 to 1 sec, the second from 1 to 2 sec… the last from 9 to 10 s).

The sample frequency is 150Hz.

Here the code block on which I’ need some help:

FREQ_BANDS = {"delta": [0.1, 5.0],
              "delta_standard": [0.5, 5.0],
              "theta": [5.0, 8.0],
              "alpha": [8.0, 13.0],
              "sigma": [13.0, 16.0],
              "beta": [16.0, 25.0],
              "gamma": [30.0, 45.0],
              }

# Optimal Window size expressed in seconds, as indicated in this post (https://mne.discourse.group/t/set-window-size-for-welch-via-compute-psd/5958/4)
n_fft_seconds = 2 / FREQ_BANDS["delta"][0] # if fmin delta is 0.5 --> 4s lengh is necessary
print(n_fft_seconds)

# Optimal Window size expressed in samples, as requested by the welch method
n_fft_sample = int(n_fft_seconds * epochs.info["sfreq"])

# PSD Welch
spectrum = epochs.compute_psd(method="welch",  picks='eeg',
                              fmin=FREQ_BANDS["delta"][0], fmax=FREQ_BANDS["gamma"][1],
                              n_fft = n_fft_sample,
                              # n_overlap = 'TO DEFINE', # The number of points of overlap between segments. Will be adjusted to be <= n_per_seg. The default value is 0.
                              # n_per_seg = 'TO DEFINE', # Length of each Welch segment (windowed with a Hamming window). Defaults to None, which sets n_per_seg equal to n_fft.
                              )

once running the script I receive the following error:

ValueError: If n_per_seg is None n_fft is not allowed to be > n_times. If you want zero-padding, you have to set n_per_seg to relevant length. Got n_fft of 3003 while signal length is 151.

QUESTIONS:

  1. Concerning the error, I’m wondering if integrating “n_overlap” and “n_per_seg” could I solve the issue, considering that I need to keep epochs of 1sec duration, and that the optimal n_fft is equal to 20s or 3003 samples (because the lower frequency of interest is 0.1) . What should I do in this situation?
  2. If the integration of “n_overlap” and “n_per_seg” solves the issue, would you suggest a formula to identify the optimal “n_overlap” and “n_per_seg” based on epoch duration (e.g. 1 sec), sample frequency (e.g. 150Hz), lower and higher frequency of interest (e.g. 0.1 and 45)? This because a formula in this direction could maybe prevent the error.
  3. In some paper, the overlap is reported in percentage (e.g. windows of 30 samples | 50% overlap | 1 Hz of frequency resolution). How to obtain such overlap percentage? In addition the frequency resolution I guess that can be calculated as:

psds_noNorm, freqs = spectrum.get_data(return_freqs=True)
freq_res = freqs[1] - freqs[0] # —> Frequency Resolution

Is it correct?

  1. In your example, you then you generate the final array (n_epochs, n_channels, n_freqs):
    # Normalize the PSDs
    psds /= np.sum(psds, axis=-1, keepdims=True) # in my case is psds = psds_noNorm / np.sum(psds_noNorm, axis=-1, keepdims=True)

    X = []
    for fmin, fmax in FREQ_BANDS.values():
        psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)
        X.append(psds_band.reshape(len(psds), -1))

    return np.concatenate(X, axis=1)

Just to make sure, could you please let me know, concerning the returning array (n_epochs, n_channels, n_freqs), what is the unit measure of each value? Both in case of no-normalized and normalize PSD. In addition, I’m not sure that each value refer to the total power for specific epoch, channel, band. My doubt is because, in your example, it is defined as “relative power”. Moreover, some author, after apply the natural log to each of this values ( to reduce skewness), are you agree?

  1. I’m focusing on Welch because it is highly used, but maybe the use of the default multitaper would be much better in this case because do not require to specify “n_fft”. Are you agree that it is more accurate than welch in term of frequency resolution and variance, where the only disadvange is the computational speed?

  2. Finally, according to my objective

For example, I may compare the theta power of the epoch (from 0 to 1 sec) with the theta power of the epoch (from 1 to 2 sec)

I did not specify any overlap while creating the epochs. Nevertheless many authors tend to overlap epochs. What would you suggest in my case?

  1. Whatever further suggestion is extremely appreciated

Many thanks in advance, I apologize for all this question in a unique post but I believe that are all helpful to clarify some basic aspects underlying the Welch PSD.

2

hi @BaggioMarco

yes Sleep stage classification from polysomnography (PSG) data — MNE 1.3.0 documentation is the best entry point here.

does this example work for you?

We call this relative power as we do

psds /= np.sum(psds, axis=-1, keepdims=True)

but you don’t need to do this to compare the PSD in one band within the same subject.

hope this helps
Alex

2 Likes
3

Thank you @agramfort , first of all I would like to thank you so much for your effort in this project and for sharing so valuable knowledge in this discourse group.
Thanks also for the clarification about the relative power.
When possible could you please let me your thoughts also about the other questions?