Hi All
I am following along on this web page which works through band power analysis and recommends welch with a window size of 4.0s. Reading the epochs compute psd documention I see i can specify additional keywords via **method_kw.
I would guess this is a dictionary e.g. {window:4.0} but this crashes because the underlying code is expecting an array.
Please could you let me know how this array should be configured.
Thank you and regards.
Robert
can you share a full code snippet to replicate your crash? so we easily see the full traceback
thanks
Alex
Hi Alex
Here is the code snippet. I am reading from the EDF dataset that you have loaded as a sample and is also available here.
EPOCH = 30.
FMIN = 1.
FMAX = 45.
print('PSG file: %s' % psgFileListing.iloc[0]['psgInputFile'])
print('Annotation file: %s' % hypnogramFileListing.iloc[0]['hypnogramOutputFile'] )
hngFile['onset'] = pd.to_datetime(hngFile['onset'])
firstOnset = hngFile.iloc[0]['onset']
hngFile['onset'] = hngFile['onset'] - firstOnset
hngFile['onset'] = hngFile['onset'].dt.total_seconds()
onsets = (hngFile['onset']).to_numpy()
durations = hngFile['duration'].to_numpy()
descriptions = hngFile['description'].to_numpy()
hngFile = mne.Annotations(onset=onsets, # in seconds
duration=durations, # in seconds, too
description=descriptions)
psgFile.set_annotations(hngFile, emit_warning=True)
tmax = EPOCH - 1. / psgFile.info['sfreq']
channels = psgFile.info['ch_names']
picks = list(filter(lambda c: match(EEG_CHANNEL_REGEX, c), channels))
events, _ = mne.events_from_annotations(psgFile,event_id=ANNOTATIONS,chunk_duration=EPOCH)
epochs = mne.Epochs(raw=psgFile,events=events,event_id=ANNOTATIONS,picks=picks,tmin=0., tmax=tmax, baseline=None, on_missing='ignore', preload=True)
###############################################
# remove frequencies above FMAX e.g. 50/60 Hz power outlet signal
# need to do for all picks
###############################################
epochs = epochs.filter(FMIN,FMAX,picks=picks)
###############################################
# compute power spectrum returns ndarray nEpochs x nChannels x nFrequencies
# following this article we need to set window=4.0 as a configuration parameter for welch
# https://raphaelvallat.com/bandpower.html
###############################################
epochPs = epochs.compute_psd(method='welch', picks=picks, fmin=FMIN, fmax=FMAX, **{'window': 4.0})
picks=['EEG Fpz-Cz', 'EEG Pz-Oz']
and here is the trace:
Connected to pydev debugger (build 212.5457.46)
## processing study edf
subject SC400
PSG file: data/sleep/sleep-edf/input/edf/SC4001E0-PSG.edf
Annotation file: data/sleep/sleep-edf/output/SC400_1_hypnogram.csv
Used Annotations descriptions: ['N1', 'N2', 'N3', 'R', 'W']
Not setting metadata
2650 matching events found
No baseline correction applied
0 projection items activated
Loading data for 2650 events and 3000 original time points ...
0 bad epochs dropped
Setting up band-pass filter from 1 - 45 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Upper passband edge: 45.00 Hz
- Upper transition bandwidth: 5.00 Hz (-6 dB cutoff frequency: 47.50 Hz)
- Filter length: 331 samples (3.310 sec)
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 0.0s remaining: 0.0s
[Parallel(n_jobs=1)]: Done 2 out of 2 | elapsed: 0.0s remaining: 0.0s
[Parallel(n_jobs=1)]: Done 3 out of 3 | elapsed: 0.0s remaining: 0.0s
[Parallel(n_jobs=1)]: Done 4 out of 4 | elapsed: 0.0s remaining: 0.0s
[Parallel(n_jobs=1)]: Done 5300 out of 5300 | elapsed: 1.2s finished
Effective window size : 2.560 (s)
Traceback (most recent call last):
File "/Users/rellis/Library/Application Support/JetBrains/IntelliJIdea2021.2/plugins/python/helpers/pydev/pydevd.py", line 1483, in _exec
pydev_imports.execfile(file, globals, locals) # execute the script
File "/Users/rellis/Library/Application Support/JetBrains/IntelliJIdea2021.2/plugins/python/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
File "/Users/rellis/dev/git/[project]/analysis/power-spectrum-measures.py", line 152, in <module>
epochPs = epochs.compute_psd(method='welch', picks=picks, fmin=FMIN, fmax=FMAX, **{'window': 4.0})
File "<decorator-gen-277>", line 12, in compute_psd
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/epochs.py", line 2031, in compute_psd
return EpochsSpectrum(
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/spectrum.py", line 1042, in __init__
self._compute_spectra(data, fmin, fmax, n_jobs, method_kw, verbose)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/spectrum.py", line 380, in _compute_spectra
result = self._psd_func(
File "<decorator-gen-199>", line 12, in psd_array_welch
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 184, in psd_array_welch
f_spect = parallel(my_spect_func(d, func=func, freq_sl=freq_sl,
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 184, in <genexpr>
f_spect = parallel(my_spect_func(d, func=func, freq_sl=freq_sl,
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 73, in _spect_func
spect = np.apply_along_axis(
File "<__array_function__ internals>", line 5, in apply_along_axis
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/numpy/lib/shape_base.py", line 379, in apply_along_axis
res = asanyarray(func1d(inarr_view[ind0], *args, **kwargs))
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 52, in _decomp_aggregate_mask
_, _, spect = func(epoch)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/spectral.py", line 737, in spectrogram
window, nperseg = _triage_segments(window, nperseg,
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/spectral.py", line 1969, in _triage_segments
raise ValueError('window must be 1-D')
ValueError: window must be 1-D
python-BaseException
Process finished with exit code 1
Hello,
epochs.compute_psd will call either Welch’s method or multitapers, depending on the argument method. Those 2 methods do not have the same ‘settings’, which have to be provided in the **method_kw as a dictionary.
Welch’s method is mne.time_frequency.psd_array_welch — MNE 1.3.dev0 documentation
The window size, the FFT window size and the overlap between window is defined by n_fft, n_overlap and n_per_segment. The size is defined in samples and not in seconds.
Thus, running Welch’s method on epochs with 4-second long windows requires:
epochs = ...
epochs.compute_psd(method="welch", n_fft=4*epochs.info["sfreq"])
Note that it also means that you need long epochs to fit many welch’s windows on which the spectrum can be estimated (or to increase the number of windows via overlap). The window length should be defined based on your frequencies of interest, 4 seconds feels a bit long.
Thank you very much Mathieu. I’m learning and appreciate the guidance. I was going by the earlier ref web page for 4 seconds which has the following text:
How do we define the optimal window duration then? A commonly used approach is to take a window sufficiently long to encompasses at least two full cycles of the lowest frequency of interest. In our case, our lowest frequency of interest is 0.5 Hz so we will choose a window of 2/0.5=4 seconds.
with 0.5Hz being the configured minimum for delta band.
I’m open to other views and ideas - thank you again very much.
i changed the code to run the provided snippet and now get a different crash:
Effective window size : 4.000 (s)
Traceback (most recent call last):
File "/Users/rellis/Library/Application Support/JetBrains/IntelliJIdea2021.2/plugins/python/helpers/pydev/pydevd.py", line 1483, in _exec
pydev_imports.execfile(file, globals, locals) # execute the script
File "/Users/rellis/Library/Application Support/JetBrains/IntelliJIdea2021.2/plugins/python/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
File "/Users/rellis/dev/git/[project]/analysis/power-spectrum-measures.py", line 152, in <module>
epochPs = epochs.compute_psd(method='welch', picks=picks, fmin=FMIN, fmax=FMAX, n_fft=4*epochs.info["sfreq"])
File "<decorator-gen-277>", line 12, in compute_psd
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/epochs.py", line 2031, in compute_psd
return EpochsSpectrum(
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/spectrum.py", line 1042, in __init__
self._compute_spectra(data, fmin, fmax, n_jobs, method_kw, verbose)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/spectrum.py", line 380, in _compute_spectra
result = self._psd_func(
File "<decorator-gen-199>", line 12, in psd_array_welch
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 184, in psd_array_welch
f_spect = parallel(my_spect_func(d, func=func, freq_sl=freq_sl,
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 184, in <genexpr>
f_spect = parallel(my_spect_func(d, func=func, freq_sl=freq_sl,
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 73, in _spect_func
spect = np.apply_along_axis(
File "<__array_function__ internals>", line 5, in apply_along_axis
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/numpy/lib/shape_base.py", line 379, in apply_along_axis
res = asanyarray(func1d(inarr_view[ind0], *args, **kwargs))
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 52, in _decomp_aggregate_mask
_, _, spect = func(epoch)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/spectral.py", line 737, in spectrogram
window, nperseg = _triage_segments(window, nperseg,
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/spectral.py", line 1965, in _triage_segments
win = get_window(window, nperseg)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 2153, in get_window
return winfunc(*params)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 1095, in hamming
return general_hamming(M, 0.54, sym)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 1017, in general_hamming
return general_cosine(M, [alpha, 1. - alpha], sym)
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 113, in general_cosine
fac = np.linspace(-np.pi, np.pi, M)
File "<__array_function__ internals>", line 5, in linspace
File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/numpy/core/function_base.py", line 120, in linspace
num = operator.index(num)
TypeError: 'float' object cannot be interpreted as an integer
python-BaseException
Process finished with exit code 1
thank you for any help you can provide.
If 0.5 Hz is your minimum frequency of interest, then I agree.
epochs.compute_psd(method="welch", n_fft=int(4*epochs.info["sfreq"]))
Those arguments are a number of samples and expect an integer. The multiplication with the sampling frequency returns a float.