read_raw_xdf(): loading assumes evenly sampling

[DominiqueMakowski]

I was adding a wrapper around pyxdf to load and tidy up xdf data in Neurokit when I stumbled upon these issues:

xdf-modules/pyxdf#79
xdf-modules/pyxdf#1

So I decided to give mnelab a try to see if what I was doing was correct, but I found an issue that happens when there was interruptions in the streaming:

import neurokit2 as nk

df, info = nk.read_xdf(
    "raw/physio/sub-01/ses-S001/eeg/sub-01_ses-S001_task-HCT_run-001_eeg.xdf"
)

from mnelab.io.xdf import read_raw_xdf

df2 = read_raw_xdf(
    "raw/physio/sub-01/ses-S001/eeg/sub-01_ses-S001_task-HCT_run-001_eeg.xdf",
    stream_ids=[1, 2, 3, 4, 5],
    fs_new=2000,
).to_data_frame()
df2.index = df.index

plt.close()
fig, ax = plt.subplots(nrows=3)
for c in ["AF7"]:
    ax[0].plot(info["data"][3][c], label="original", color="blue")
    ax[1].plot(df[c], label="NK", color="green")
    ax[2].plot(df2[c], label="mnelab", color="orange")
plt.legend()
plt.show()

Essentially I plotted:

• in blue, the original signal from Muse EEG sampled at ~250He
• in orange, the resampling at 2000 Hz by NK based on pandas that is timestamp-based
• in blue the one by mnelab

There is an interruption at the beginning at the stream, but mnelab probably interpolates linearly which distorts the whole thing.

I can send you that xdf file by email if you need :)

(also tagging my friend @pmavros as this might be relevant to our eeg processing)

[cbrnr]

Thanks @DominiqueMakowski, you are correct that we currently assume regularly spaced samples per stream. Using pandas to handle interpolation is actually very clever, I wonder why I haven't thought of it before 😆 (the tradeoff is of course that it is a rather large dependency, but maybe worth it even for MNELAB).

Could you share the file with me so that I can play around with it to get a better grasp of the problem?

Regarding your function, how can I avoid linear interpolation for such a long interruption? By setting fillmissing=1/fs (where fs is the expected sampling frequency of the stream)?

Finally, I think it would be beneficial to get resampling directly into pyXDF. Did you check the implementation in xdf-modules/pyxdf#1 by any chance?

[DominiqueMakowski]

Could you share the file with me

Dropped you an email

Regarding your function, how can I avoid linear interpolation for such a long interruption?

I added a fillmissing argument that sets a limit in seconds on wether to prevent interpolations of long periods. It seem to work nicely and leaves NaNs, which is not the best for MNE but that's another issue

Finally, I think it would be beneficial to get resampling directly into pyXDF.

I agree...

[cbrnr]

I added a fillmissing argument that sets a limit in seconds on wether to prevent interpolations of long periods

So I set it to 1/fs? The default None seems to turn off NaNs completely if I'm not mistaken.

[DominiqueMakowski]

It does the transformation above so you just need to specify it in seconds, like 0.5 or 0.1. And yes None by default leaves all interruptions (normally)

[DominiqueMakowski]

If we are to have similar code in NK and mnelab, we might want to outsource some of it to pyxdf.

What I could see is a pyxdf.xdf_to_dataframe() function that tries to import pandas (and if not, errors saying that pandas is required for this function) - i.e., pandas is an optional dependency, and then converts and resamples the xdf to a dataframe (like what we do in NK now). Then, mnelab can extract the array from it and put that into mne.Raw and NK can use it as-is.

[cbrnr]

If we are to have similar code in NK and mnelab, we might want to outsource some of it to pyxdf.

Agreed. Resampling should really be handled by pyXDF, and a proposed solution already exists (although I'm not sure how easy it will be to rebase and if it is still working). But this should be discussed directly with the pyXDF people (maybe in xdf-modules/pyxdf#1).

[DominiqueMakowski]

For future reference: that method might suffer from some loss of precision, from my small experimentations using the union of existing and new indices was giving the best results

[cbrnr]

@DominiqueMakowski is the description of the signals in your top post correct? I think you might have mixed up the colors.

Just to be sure, the correct (expected) signal should contain a segment with missing data in the first second?

When loading just stream 4 (with or without resampling), I get this time series:

So I'm wondering if the import worked, and the problem is maybe in the to_data_frame() method? I can't investigate further for a couple of days (weeks?), so if you find anything in the meantime please LMK.

[cbrnr]

One more observation, the MNELAB GUI doesn't let you choose sampling frequencies greater than the highest sampling frequency in the file, i.e. 1000Hz in this example. Even then, the signals look exactly like in the screenshot, so maybe it's because you resample to 2000Hz (I doubt it, but still worth checking)?

[cbrnr]

Although I did spot another (probably unrelated) issue: the suggested resampling when selecting all streams is 52Hz, but it should be 1000Hz. Not sure what's going on here, but this seems like a separate issue.

[DominiqueMakowski]

Just to be sure, the correct (expected) signal should contain a segment with missing data in the first second?

No, I think the whole recording is like several minutes so it should be within the first minute or so (the time axis is messed up in my fig)

so maybe it's because you resample to 2000Hz (I doubt it, but still worth checking)?

The upsampling is done to avoid aliasing when merging signals with uneven sampling rates, but it should have fairly minimal impact

[cbrnr]

No, I think the whole recording is like several minutes so it should be within the first minute or so (the time axis is messed up in my fig)

So your three example plots do not actually show the problem? Sorry, I'm confused now, but now I don't understand what the problem with MNELAB is...

[DominiqueMakowski]

Can you zoom out in your fig to see all the signal horizontally?

[cbrnr]

Yes:

Looks the same whether I resample to 256Hz or not (in which case the sampling rate is 256.021Hz).

[DominiqueMakowski]

can you share the code to reproduce this fig?