Regarding mne.time_frequency.psd_welch function, I cannot seem to understand how the average argument is affecting the Welchâ€™s power spectra estimate - the results are the same for me with â€˜meanâ€™ and â€˜medianâ€™ passed as averaging methods. I was hoping that using median would be more effective in case some epoch still has artefacts in them and maybe these artefacts wouldnâ€™t affect the final PSD for that epoch that much.

And in overall, is this the correct way to estimate PSD for epochs in EEG or do you recommend any other ways?

sfreq = epochs.info['sfreq'] # sampling frequency in my case is 1024
epoch_duration = 2 # each epoch is 2s long
psds, freqs = mne.time_frequency.psd_welch(epochs,n_fft=int(sfreq*epoch_duration),fmin=1,fmax=100,average='mean')

Thanks for the quick answer and thank you for creating this immersive library!
So you mean for example to have the epoch duration 4 seconds (4096 samples) and in psd_welch function argument n_fft 2 seconds (2048 samples)? Is there a specific rationale for this which you could explain me?

if the epochs are short enough that only 1 welch window is used, then there is no averaging to be done, so both mean and median will return the same thing, e.g., the estimate from the one single welch window.