lrkrol / SEREEGA

SEREEGA: Simulating Event-Related EEG Activity
64 stars 17 forks source link

Timing is off by one sample #10

Closed behinger closed 4 years ago

behinger commented 5 years ago

Still, great toolbox. But the timing is off by one sample I think.

Reproducible code:

lf = lf_generate_fromnyhead('montage', 'S64');
p1_comp = utl_get_component_fromtemplate('visual_p100_erp',lf);

epochs = struct();
epochs.n  = 10;             % the number of epochs to simulate
epochs.srate = 20;
epochs.length = 1000*0.5;
data = generate_scalpdata(p1_comp, lf, epochs,'showprogress',0);
EEG = utl_create_eeglabdataset(data, epochs, lf, 'marker', 'event1');

p1_comp(1).signal{1}.peakLatency
mean(EEG.data(63,:,:),3)
EEG.times

Output:

p1_comp(1).signal{1}.peakLatency

   100

That is, I would expect to see a peak at 100ms.

mean(EEG.data(63,:,:),3)

  1×10 single row vector

    0.0005192576       0.7032224    2.322604e-05    7.963054e-19    7.701536e-42               0               0               0               0               0

I find a peak at the second data entry

EEG.times

     0    50   100   150   200   250   300   350   400   450

But the second data entry is 50 not 100ms. It seems there is a mistake here

Best, Bene

EDIT: I noticed there might be jitter in the p1, but I don't think the jitter has influence on the problem

lrkrol commented 4 years ago

Good find! Seems I forgot compensate for MATLAB's one-based indexing.

Testing this with a peak of a single sample, we see that an ERP at 1000 Hz with a 1-ms-wide peak at 100 ms has this peak at index 100:

[~, i] = max(erp_generate_signal(100, 1, 1, 1000, 1000))

i =
   100

But indeed we tend to start counting at 0 when dealing with epochs, so this should be 101.

Therefore, I moved the ERP latency over one sample in erp_generate_signal on line 59:

peakLatency = (peakLatency/1000)*srate + 1;

That should fix the issue.

The burst latency of the ERSP had the same issue, so I applied the same fix in ersp_generate_signal on line 156.

[~, i] = max(ersp_generate_signal(1, 1, 0, 1000, 1000, 'modulation', 'burst', 'modLatency', 250, 'modWidth', 1))

I believe those are the only scripts that were affected by this.

Thanks for the report -- the issue is significant enough to make a new release out of the current code. If you found anything else, now would be a good time to let me know. :D

behinger commented 4 years ago

great :-) No, this is the only thing I found