Critical Bug in _jitter_removal

[ollie-d]

Hello,

Please forgive me if I am posting this in an undesired format. I simply want to explain the bug, how I found it, and what the solution is.

The bug When jitter removal is enabled, the current code creates an invalid mapping, thus distorting time. I am not sure if the effect is as dramatic as what I have experienced, but in my case it renders my markers completely useless.

Problematic code (line 563 of pyxdf.py): indices = np.hstack((0, indices, indices, nsamples - 1))

How I Found it I collected some pilot data and realized that my classifier was at chance. After going through all of the relevant code, I determined that my filtering and classifying were being conducted correctly. I then ran my analysis on simulated data, and the result was perfect. I then started analyzing my code for epoching the data. Luckily, I have started collecting photosensor channels so as not to rely too heavily on software markers. When I examined the markers position relative to the photosensor, I realized they were not aligned: I then decided to look at my data in MATLAB via the load_xdf() function shipped with EEGLAB, and my data looked as expected: After that, I decided to individually change all of the options in python and discovered that when dejitter_timestamps was set to False, my data looked as expected:

After determining that the removing of jitter was the problem, I went through each line of the relevant code in MATLAB and python and discovered that the ranges being used were inconsistent. For example, with my data, MATLAB first finds a problem in the range [2 11]. I expected to see that python would show the same problem in [1 10], but instead it finds a problem at [0 10]. I noticed that every index that python found had a lower bound by 1 index compared to MATLAB, even when taking 1-indexing into account: This subsequently caused a discrepancy in the mapping values between python and MATLAB.

Solution Increment the first index by 1 to match MATLAB. indices = np.hstack((0, indices+1, indices, nsamples - 1))

Conclusion An easy fix for a problem that had me question everything. Please let me know if you would like my data to test this solution on. Again, I'm not certain this is always a problem, but it certainly was in my case. The setup I'm using is an actiCHamp streaming data via the RDA with 64 + 8AUX at 500 Hz. Markers are sent using pylsl. Photosensor is activated on the same frame as my stimulus display, and I have done extensive testing to validate its precision and accuracy.

Thank you, Ollie

[cbrnr]

Thanks for this important find! Could you share a minimal test example (i.e. data and code) to reproduce the problem? I'd be very interested why there are such large discrepancies in your case, because I've never seen this with my data (although I have to admit that I don't have a lot of data in the first place). It seems like this could be a copy and paste error (I guess the MATLAB code was ported to Python), so we should make sure that there are no additional indexing errors - but since your data looks OK after the change this shouldn't be the case.

[ollie-d]

Hey Clemens,

Attached below you will find the dataset and Python code I used to discover the problem. The reason why I was so confused with this error and why it took me so long to suspect something was amiss was this issue has not been present in other datasets I've collected and analyzed with pyxdf. The precursor to this study used the BrainAmp connector along with SNAP, integrated with LabRecorder, and the distortion is not present; although it cannot be tested the same way since there is no photosensor data. With this iteration, I'm using actiChamp via PyCorder + BrainVision RDA integrated with LabRecorder, and I'm using PsychoPy for stimulus display along with a pylsl marker stream. It is worth noting that I set the marker stream to be sampled because when I compared jitter between unsampled/sampled in EEGLAB, sampled was more stable.

Thank you for taking the time to assist with this bug, and please let me know if you need any further information.

Best, Ollie

Link to data/code: https://drive.google.com/open?id=1Q4avYcydH-x_aNRLFn1D2LXJ764S12JV (Too big for GitHub)

[cbrnr]

Thanks a lot! I can reproduce the problem, and I think I understand what's going on.

Your marker stream is regularly sampled at 10kHz according to the info dict. However, within the recorded time range of about 686s there are only 396 data points, which translates to a true sampling frequency of roughly 0.58Hz (i.e. a marker occurs once every 1.73s on average). The dejitter algorithm works on the basis of the reported nominal sampling frequency of 10kHz, which means that when it detects breaks greater than the maximum of 1s or 500 samples (which is just 0.05s), it will correct/linearize the time stamps in between. In your case, there are 71 segments between markers greater than 1s, therefore the correction kicks in quite often. If you had 10000 samples per second as per the reported sampling rate, the algorithm would most likely never kick in.

Still, the issue you found is most likely legit. In cases where the reported sampling frequency really matches what is actually in the data, the error will be very small. This is probably why no one has noticed this yet.

Side question: why didn't you use an irregular marker stream (i.e. sampling frequency 0)? Most marker streams are irregularly sampled (and thus do not get dejittered), which might be an another reason why no one noticed this problem so far. I tested this with your data set, that is, when I change the <nominal_srate> tag to 0 (make sure the number of characters stays the same), the markers are at the expected locations even if dejitter_timestamps=True.

It would be great if you submitted a PR with the one line fix so we can properly credit you with this contribution. Furthermore, it would be very helpful if @cboulay or @tstenner double-checked that this fix really does the right thing.

[ollie-d]

Thank you for your explanation behind the onset of the bug in this dataset versus others. I suspected it could be something related to the sampled marker stream, but was not certain.

In response to your side question, I have been using unsampled marker streams for the last few years, and never had a problem with them. However, when deciding which software I was going to use to create my LSL streams for my current study, I discovered a significant delay between software marker onset and photosensor offset. I stumbled upon the ActiChamp specific LSL app and a document suggesting that sampled markers performed better than unsampled markers. I understood that this mode was viable with hardware triggered markers, but decided to use sampled markers to see if I could have any improvements. I ran some simple tests and the sampled markers seemed to perform better, however, I just re-conducted the test and there does not seem to be any affect: (Time series photosensor data collected/sent with actiCHamp via PyCorder + BrainVision RDA, marker data sent using pylsl, integrated using LabRecorder, epochs extracted via EEGLAB)

I now realize that my previous assumption about sampled markers was incorrect, i.e. I thought that the stream would send an empty marker every 1/srate and my marker stream would subsequently be functionally sampled at 10kHz, but this is of course not the case.

35

I will submit the PR after making this post. Thank you again for all of your assistance, and for all of your combined efforts for making such high quality software openly available to the masses :)

[cbrnr]

I never saw such large variability in irregular marker streams as described in the document you linked. Specifically, in your marker stream it makes no difference whether or not dejittering is performed (after your bugfix of course), so it is already very accurate. However, this might be due to the fact that you do have a sampled marker stream. On the other hand, the contents of the marker stream is clearly irregular, so IMO the specified nominal sampling rate of 10 kHz is incorrect.

This might be an issue specific to BrainProducts amplifiers - maybe high accuracy is only possible with sampled marker streams, but this technical detail should be hidden from the user and the XDF file should denote such streams as irregular marker streams. Otherwise, there should be data points even when no markers are present.

I'm thinking of whether it might be a good idea to issue a warning if we detect streams with vastly different nominal vs. effective (i.e. computed via number of samples divided by length in seconds) sampling frequencies.

@dmedine what do you think?

[tstenner]

On the other hand, the contents of the marker stream is clearly irregular, so IMO the specified nominal sampling rate of 10 kHz is incorrect.

There's three marker types: read from the device alongside the EEG data in the same stream (1), in a different stream (2) and sampled independently and sent in a different stream (3). (1) needs a sampling rate of 10KHz, but the others don't.

[cbrnr]

I don't understand. Why does (1) need fs=10kHz when the EEG is sampled at 500Hz? A stream can only have one sampling rate. In any case, if the stream says that it has fs=10kHz I expect 10000 samples per second and not samples at irregular intervals as is the case in this example data.

[tstenner]

Ah, via the RDA. I meant the EEG stream's sampling rate.

[cbrnr]

OK, this makes sense. But do you agree that for the marker stream in this example, the nominal sample frequency should be 0 and not 10000?

[cboulay]

I'm thinking of whether it might be a good idea to issue a warning if we detect streams with vastly different nominal vs. effective (i.e. computed via number of samples divided by length in seconds) sampling frequencies.

Yes that should happen.

[cbrnr]

Fixed in #35.