1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
function gdf_remove_eog( inputfile, outputfile, reffile, eegchan, eogchan )
% gdf_remove_eog( inputfile, outputfile, reffile )
%
% Remove eog artifacts by regression.
%
% Inputs:
% inputfile : filename of the original raw data
% outputfile : output filename
% reffile : file that contains eog reference data
% eegchan : EEG channel numbers
% eogchan : EOG channel numbers
% apply spatial filter matrix to a gdf
% load gdf files
[s, h, e] = gdf_reader( inputfile, 'dataformat', 'single' );
input.s = s;
input.h = h;
input.e = e;
[s, h, e] = gdf_reader( reffile, 'dataformat', 'single' );
ref.s = s;
ref.h = h;
ref.e = e;
% calculate eog regression
refeeg = [ref.s{eegchan}];
refeog = [ref.s{eogchan}];
refeeg( any( isnan( refeeg ), 2 ), : ) = [];
refeog( any( isnan( refeog ), 2 ), : ) = [];
Y = refeeg;
U = refeog;
CNN = U'*U/size(U,1);
CNY = U'*Y/size(U,1);
eog_weights = CNN\CNY;
% apply eog correction
eeg = [input.s{eegchan}];
eog = [input.s{eogchan}];
eeg = eeg - eog * eog_weights;
% 4. save new gdf
handle = gdf_writer( 'init' );
for c = 1 : input.h.file.num_signals
gdf_writer( 'createsignal', handle, c );
end
gdf_writer( 'setheader', handle, input.h );
gdf_writer( 'recordduration', handle, 0 ); % automatic record duration
gdf_writer( 'eventconfig', handle, input.e.mode, input.e.sample_rate );
gdf_writer( 'open', handle, outputfile );
chunksize = 1; % seconds
num_chunks = ceil( length(input.s{1}) / (chunksize*input.h.signals(1).sampling_rate) );
datapos = ones(1,input.h.file.num_signals);
for d = 1 : num_chunks
for c = 1 : input.h.file.num_signals
datalen = chunksize * input.h.signals(c).sampling_rate;
if find( eegchan==c )
data = input.s{c}( datapos(c) + (1:datalen) - 1 );
else
data = input.s{c}( datapos(c) + (1:datalen) - 1 );
end
gdf_writer( 'blitsamples', handle, c, data );
datapos(c) = datapos(c) + datalen;
end
end
if input.e.mode == 1
for e = 1 : length( input.e.position )
gdf_writer( 'mode1ev', handle, input.e.position(e), input.e.event_code(e) );
end
elseif input.e.mode == 3
for e = 1 : length( input.e.position )
gdf_writer( 'mode3ev', handle, input.e.position(e), input.e.event_code(e), input.e.channel(e), input.e.duration(e) );
end
end
gdf_writer( 'close', handle );
gdf_writer( 'clear', handle );
end
|
