How to remove artifacts using change cut offs

Hi, I have 60000 measurements of pupil diameter recorded at 120Hz, but it includes artifacts where the diameter rapidly drops in size due to blinks or etc. How could I design a filter, or write code to delete any data point that changed quickly between consecutive measurements? E.g. reduced or increased in value by 5 over two measurements. Thank you. Ben.

Answer 1

[Walter Roberson]

You could examine abs(diff(YourMeasurements)). Or perhaps you want to do median filtering.

[Ben Jaquiery]

Thanks Walter, I'm new to Matlab, but by using

>> pupildif = abs(diff(pupil);

I get a difference value between each pupil diameter measurement, 'pupildif'.

Is there then a way to use 'pupildiff' to remove the corresponding datapoints in 'pupil' (change the value to 0) if the value in 'pupildif' is 4 or larger, and then remove 3 datapoints on either side of the removed value?

Thanks

Ben

[Image Analyst]

When your pupil size changes, does it change in a way that there is a 4 or larger change from one element to the next? If so, diff is good. If it changed by 4 over 3 elements, then diff() would not be good. Why don't you give us some data during the bad/blink time plus some "good" measurements on either side of the "bad" data? Plus give us a screenshot of a plot of the data. Also removing 3 is not very robust. You should remove exactly how many you need, whether it's 1 or 3 or 5 or 7 or whatever.

[Walter Roberson]

diff(pupil,2) could be used for distances further apart. Or code

abs(pupil(3:end) - pupil(1:end-2))

[Ben Jaquiery]

I will post a screenshot tomorrow. The pupil diameter changes over time slowly (less than 2 over 1 element), but every so often an artifact will occur like this: 85 84 85 84 84 84 80 73 55 55 56 57 54 50 60 61 55 54 67 74 80 84 84 84 85 86 85. I'd like to be able to filter these artifacts out, i.e. change the numbers in the 50 and 60s to zero, and then also remove the numbers contaminated by the artifact, i.e. 80, 73 and 74 and 80 on each side.

Because the filter will eventually be used to on thousands of blinks, an estimate of removing 3 data points on each side of the artifact will be OK.

[Image Analyst]

One way to handle that is a Savitzky-Golay filter. It's like a moving window where the center element is replaced by the value of a polynomial fitted to values inside the window. It's done by sgolay() in the Signal Processing Toolbox. Let me know if you want a demo.

[Ben Jaquiery]

Great thanks Image Analyst, the Savitzky-Golay filter looks like it could be very useful. I will try it tomorrow and let you know if I need a demo. Thanks

Answer 2

[Image Analyst]

Do you have some known size, below which is junk data? If so, just threshold it out.

badDataIndexes = pupilDiameters < minAllowablePupilDiameter;
pupilDiameters(badDataIndexes) = []; % Remove bad data
% Need to remove those times also because now all 
% elements are not the same time difference apart
% like they (presumably) used to be.
timePoints(badDataIndexes) = [];