def edf2psd(edf, config, *, stats=None, verbose=None, logger=pu.stdout):
"""Extract the bands from the Raw EDF data.
Remove flatlines using annotate_amplitude() and expand them.
Process each channel one-by-one so that a flatline on one channel
does not affect the others.
Limit processing to channels.
Use welch method, limiting frequencies to [fmin,fmax].
Ensure frequency step is at most .25Hz."""
start_time = pu.start_time()
edf_df = edf.to_data_frame()
sfreq = edf.info["sfreq"]
psdd = {}
freqs = None
bad = {}
data_end = edf.tmax
dropped = []
for ch in config.channels or edf.ch_names:
ra = mne.io.RawArray(edf_df[ch].values.reshape((-1, edf_df.shape[0])),
mne.create_info(ch_names=[ch], sfreq=sfreq, ch_types="eeg"),
verbose=verbose)
# https://mne.tools/stable/generated/mne.preprocessing.annotate_amplitude.html
anns, bads = mne.preprocessing.annotate_amplitude(
ra, flat=0, bad_percent=100, verbose=verbose)
if bads:
bad.setdefault("annotate_amplitude",[]).append(ch)
dropped.append(1)
continue
if anns: # non-trivial
anns = expand_annotations(anns, data_end=data_end, min_good=config.min_good)
if anns is None:
bad.setdefault("expand_annotations",[]).append(ch)
dropped.append(1)
continue
ra.set_annotations(anns)
dropped.append(anns.duration.sum()/(edf.tmax - edf.tmin))
else:
dropped.append(0) # nothing was dropped!
# https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw.filter
ra = ra.filter(l_freq=config.l_freq, h_freq=config.h_freq,
skip_by_annotation="BAD_", verbose=verbose)
# https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw.compute_psd
# n_jobs > #channels -- does not yield any performance benefit
# https://github.com/mne-tools/mne-python/pull/11298
# https://github.com/mne-tools/mne-python/issues/11297
spectrum = ra.compute_psd(method="welch", fmin=config.fmin, fmax=config.fmax,
n_fft=2048, n_jobs=1, verbose=verbose)
psds, freqs0 = spectrum.get_data(return_freqs=True)
psds0 = np.isnan(psds)
if psds0.any():
bad.setdefault("spectrum NaN",[]).append(ch)
if not psds0.all():
raise ValueError("edf2psd: some but not all PSD are NaN",edf,ch,psds)
else:
if freqs is None:
freqs = freqs0
if np.diff(freqs).mean() > 0.25:
raise ValueError("edf2psd: freq step too high",
edf.info, np.diff(freqs).mean())
elif not np.array_equal(freqs,freqs0):
raise ValueError("edf2psd: inconsistent freq",edf,ch,freqs,freqs0)
psdd[ch] = psds.flatten()
if not psdd:
raise ValueError("edf2psd: no useful data", edf, bad)
df = pd.DataFrame(psdd, index=freqs)
df.index.name = "Freq"
elapsed = pu.elapsed(start_time)
logger.info("edf2psd: %s [%s]", pu.shape(df), elapsed)
if stats is not None:
for d in dropped:
stats.dropped.add(d)
stats.run_time.add(elapsed)
stats.edf_len.add(len(edf))
if bad:
logger.info(" %d failed:%s", sum(map(len,bad.values())), "".join(
f"\n {kind} {len(chans)}: {chans}" for kind,chans in bad.items()))
return df
