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
417 stars 113 forks source link

Conversion to microvolts #59

Closed matiasandina closed 2 years ago

matiasandina commented 2 years ago

I am working with a np.array that I convert to mne.raw to pass to SleepStaging(). I was getting huge power values and I traced it back to this line:

https://github.com/raphaelvallat/yasa/blob/71c0a8245c61b328a82ab1197c34b673333120a3/yasa/staging.py#L200

In the comments for SleepStaging(), you have

    .. important:: The PSG data should be in micro-Volts. Do NOT transform (e.g. z-score) or filter
        the signal before running the sleep staging algorithm.

The data that I am entering is in microvolts, and the transformation to raw is not affecting that

eeg
array([[-76., -29., -30., ...,  -5., -12., -44.]])
raw_array.get_data()[0]
array([-76., -29., -30., ...,  -5., -12., -44.])

The sls object has the scaled data stored

sls.data[0]
array([-75994247.75916173, -37654850.96502215, -54616406.07807892, ...,
       -65886683.36596085, -81338678.52364798, -19189291.22413166])

If I calculate the spectrum by hand, I obtain the same result either using the eeg object or the raw_array.get_data(), but if I use the * 1e6 factor I get an obviously large number.

>>> simps(scipy.signal.welch(eeg, 512, nperseg=512*2)[1])
21057.190104166664
>>> simps(scipy.signal.welch(raw_array.get_data()[0], 512, nperseg=512*2)[1])
21057.208596347984
>>> simps(scipy.signal.welch(raw_array.get_data()[0] * 1e6, 512, nperseg=512*2)[1])
2.105720859634799e+16

Should there be a way to indicate to SleepStaging() that it shouldn't rescale the data ? I am thinking

sls = yasa.SleepStaging(..., units=None) # rescales as before
sls = yasa.SleepStaging(..., units="uV") # does not rescale

This issue might or might not be there for the other bands because the default behavior is:

yasa.bandpower_from_psd_ndarray(..., relative=True)

But if I understand the code correctly, it might be also calculating a scaled spectrum and then doing the relative scaling back.

raphaelvallat commented 2 years ago

Hi @matiasandina,

Thanks for opening the issue. Great point. I think that a quick and dirty fix is just to divide your Raw data by 1e6, e.g. raw._data /= 1e6 (I haven't tested this).

In the long term, you're right that we need a better solution. As a matter of fact, mne.Raw.get_data() now has an units parameter that will automatically do the conversion, if needed:

raw.get_data(units="uV")

raphaelvallat commented 2 years ago

TODO for next release of YASA