Is it possible to calculate the data in epochs.plot_psd separate from plotting?

💬 Support & Discussions
, ,
1

All,

I am trying to extract the power spectrum from hundreds of PSG files for an analysis project. I have the pipeline built to extract the sleep staging into a csv per subject. I am now trying to extract the power spectrum into a csv file for later analysis.

given an edf raw file and annotations file and defined frequency bands I know I can generate charts using this code:

FREQ_BANDS = {"delta": [0.5, 4.5],
              "theta": [4.5, 8.5],
              "alpha": [8.5, 11.5],
              "sigma": [11.5, 15.5],
              "beta": [15.5, 30.0]}

psgEpochs.plot_psd(fmin=0., fmax=FREQ_BANDS['beta'][1], average=True, spatial_colors=False)

and get a chart like this:

image
image1264×430 58.6 KB

however the data is encapsulated in the chart object. Realizing that the plot_psd function uses a multi-taper I then run this code:

 psds, freqs = mne.time_frequency.psd_multitaper(
        psgFile,
        fmin=0,
        fmax=FREQ_BANDS['beta'][1],
        picks=psgInfo['ch_names'][0:2]
    )

    fig, subplot = plt.subplots()
    fig.set_size_inches(11., 8., forward=True)
    for psd in psds:
        subplot.plot(freqs,psd,color='k')
    plt.show()

and get a chart like this:

image
image1170×886 31.4 KB

I realize I am new to such work and very much learning as I go. Accepting this limitation, please could someone let me know how I can get the data shown in psgEpochs.plot_psd independent of plotting a chart.

Thank you for any help you can provide.

regards

Robert

2

Hello, I think you’re almost there, you just need to convert the data to a dB scale via

psds_dB = 10 * np.log10(psds)

and then plot these values. I believe this should give you the correct results.

Best wishes,
Richard

3

Thank you very much @richard - most helpful. I had not appreciated that the last step was as simple as it is.

regards

Robert

1 Like