raphaelvallat / yasa

YASA (Yet Another Spindle Algorithm): a Python package to analyze polysomnographic sleep recordings.
https://raphaelvallat.com/yasa/
BSD 3-Clause "New" or "Revised" License
428 stars 115 forks source link

IRASA- periodic power values #136

Closed apoorva6262 closed 1 year ago

apoorva6262 commented 1 year ago

Hi @raphaelvallat ,

I used the following function to determine the periodic power from the signal : freqs, psd_aperiodic, psd_osc, fit_params = yasa.irasa(sig, Fs, band=[3,60],kwargs_welch={'average': 'median'})

This gave me the output for average periodic power (np.nanmean(psd_osc)) as 4.913666773031852e-14.

I tried neurodsp package (https://neurodsp-tools.github.io/neurodsp/auto_tutorials/aperiodic/plot_IRASA.html?highlight=irasa) with same settings and I got the average periodic power to be 1.3315638777465248e-13. Wondering what could cause the difference? I used the same welch and median method to determine the psd.

Additionally how do you determine the periodic power for a subset of band frequencies like alpha from the initial fit ? Right now I have set the band parameter to 3-60hz, but how do I extract alpha frequencies from the fit itself ? I did not want to set the initial band paramaters to 8,12hz as that is a small range to find the aperiodic fit. I tried indexing but that gave me some shape issue.

Thanks, Apoorva.

raphaelvallat commented 1 year ago

Hi @apoorva6262,

Regarding your first question, these two values are actually both zero:

np.isclose(4.913666773031852e-14, 1.3315638777465248e-13)  # True
np.isclose(4.913666773031852e-14, 0)  # Also true

For your second question, I would recommend using the yasa.bandpower_from_psd function (https://raphaelvallat.com/yasa/build/html/generated/yasa.bandpower_from_psd.html):

yasa.bandpower_from_psd(psd_osc, freqs, band=[(8, 12, "Alpha")])

Hope this helps, Raphael

apoorva6262 commented 1 year ago

@raphaelvallat , For the second question, I replaced psd for 'psd_osc' in 'yasa.bandpower_from_psd' as I am interested in periodic peaks only :

freqs, psd_aperiodic, psd_osc, fit_params = yasa.irasa(sig, Fs, band=[3,60],kwargs_welch={'average': 'median'})
yasa.bandpower_from_psd(psd_osc, freqs, ch_names=['F7'], bands=[ 4, 8, 'Theta', 8, 12, 'Alpha',  13, 30, 
                                                         'Beta'], relative=True) 

The shape of psd_osc is (1,229) , freqs is (229), both are numpy arrays. But I get this error: TypeError: 'int' object is not subscriptable

Screen Shot 2023-02-13 at 11 18 35 AM
raphaelvallat commented 1 year ago

Please try the following:

yasa.bandpower_from_psd(psd_osc, freqs, ch_names=['F7'], bands=[(4, 8, 'Theta'), (8, 12, 'Alpha'), (13, 30,  'Beta')], relative=True) 
apoorva6262 commented 1 year ago

Perfect thanks! One last question, a lot of the PSD values are negative and I did see your post on how to handle it. But I was wondering if you had come up with any other solution besides the ones you have mentioned here https://github.com/raphaelvallat/yasa/issues/29 ?

raphaelvallat commented 1 year ago

Unfortunately I do not have new solutions to this problem :/ Have you tried using FOOOF? It is also a great alternative to IRASA.

apoorva6262 commented 1 year ago

No worries, I did try FOOOF and that is great. I am actually comparing IRASA to FOOOF

apoorva6262 commented 1 year ago

@raphaelvallat , Another question, I am running this IRASA function to determine periodic power. For some files, I see that the fit(r2 value) is poor. Do you have any suggestions as to how to exclude poor fits or any outlier methods you employ for further statistical analysis?

raphaelvallat commented 1 year ago

I have not used IRASA in a long time for my own analyses so it's hard to say. I would probably not exclude a night just based on a poor R^2 value. But I do recommend plotting the full, aperiodic, and periodic spectra to get a better idea of what may be happening.

apoorva6262 commented 1 year ago

Thanks! With my dataset, I am determining power across all 32 eeg channels for three different groups. For consistency sake, I have made sure to use the same settings of IRASA for all the channels and for all subjects in different groups. But there are 1-2 subjects with poor fits due to noise in higher frequencies (r2=0.002) , would you recommend just reporting the noise? It wouldn't make sense to use different settings for each file/subject to get a better fit.

raphaelvallat commented 1 year ago

Can you try to filter out the high-freq noise for these 1-2 participants? If the power spectrum is really too noisy then I don't think it's worth calculating IRASA honestly, the results will likely be biased.

apoorva6262 commented 1 year ago
Screen Shot 2023-03-01 at 10 37 53 AM

Also , I have this issue of negative psd issue with vast majority of subjects, i have used the default resampling method.