Human neuronal firing varies with the frequency of local field potential oscillations.

1. function [bad_data,bad_percent ] = Find_Epileptic_signal(signal,params )
2. %identifies portions of signal that show epilepsy-related artifuct
3. %inputs: signal is raw signal (usually LFP)
4. %inputs: srate is sampling rate in Hz
7. %outputs bad_data, a boolean indicating badness(==1) at each sample
8. %outputs bad_percent, percent of signal that was detected as bad
10. %based on the methods described by Aghajan Ambulatory Theta preprint
11. %also based on Gelinas....Buzsaki 2016 Nature Comms paper
13. srate = params.srate;
14. unfilt_sd_thresh = params.bad_unfiltered_stdev_thresh; %5 by default
15. filt_sd_thresh = params.bad_filtered_stdev_thresh; %6 by default from aghajan
19. %mean center signal
20. signal = signal-nanmean(signal);
22. %Nov 6
23. %low pass filter signal first to get rid of any possible spike artificats
24. %these very brief spike artifacts threw off the algorithm and detected all
25. %spikes as bad events in some cases!
26. %[signal] = buttfilt(signal,[100],srate,'low',4);
27. [signal] = buttfilt(signal,[80],srate,'low',4);
31. %filter signal 25-80Hz
32. [filter_signal] = buttfilt(signal,[25 80],srate,'stop',4);
34. %add 10/31
35. filter_signal = filter_signal-nanmean(filter_signal);
38. %get envelope of signal and filter signal
39. hilb_sig = hilbert(signal);
40. signal_envelope = abs(hilb_sig);
42. hilb_filt = hilbert(filter_signal);
43. filt_envelope = abs(hilb_filt);
45. %
46. % %trust but verify
47. % figure;plot(filter_signal(1:2000));
48. % hold on;
49. % plot(filt_envelope(1:2000),'r')
51. crit_unfilt = std(signal_envelope)*unfilt_sd_thresh;
52. crit_filt = std(filt_envelope)*filt_sd_thresh;
54. bad_unfilt = find(signal_envelope>crit_unfilt);
55. bad_filt = find(filt_envelope>crit_filt);
57. %output a binary vector of badness
58. bad_data = zeros(1,length(signal));
59. bad_data(bad_unfilt)=1;
60. bad_data(bad_filt)=1;
64. %also throw out very small gaps in bad data by finding these little "good"
65. %windows between bad windows and making sure they are at least one second
66. %long.
67. %if they aren't at least a second, mark em bad
68. good_windows = bwlabel(~bad_data);
69. un_win = unique(good_windows);
70. for oG = 2:length(un_win); %start at 2 b/c 0 is first value
71. good_win_idx = find(good_windows==un_win(oG));
72. good_win_length(oG) =sum(good_windows==un_win(oG));
73. if good_win_length(oG)<srate; %if its not at least one second of good data, make it a bad window
74. bad_data(good_win_idx) = 1;
75. end
76. end
78. bad_percent = (sum(bad_data)./length(signal)).*100;
80. %
81. % %wanna visualize what got caught?
82. % bad_win = bwlabel(bad_data);
83. % unbad = unique(bad_win);
84. % figure;
85. % for iBad = 2:17%length(unbad)
86. % subplot(4,4,iBad-1)
87. % bad_idx = find(bad_win==unbad(iBad));
88. % badd = bad_idx(1)-srate/2:bad_idx(end)+srate/2;
89. % [ax]=plotyy(1:length(badd),signal(badd),1:length(badd),bad_data(badd))
90. % xlim(ax(1),[1 length(badd)])
91. % xlim(ax(2),[1 length(badd)])
92. %
93. % end

Find_Epileptic_signal.m, under CC-BY-4.0 · at the source