Authored by Anthony Boyer

Fix conversion error: "in_fread_edf.m" was the same function as in "in_fopen_edf…

…"; add "struct_fix_events.m" for FRE patients.
1 -function [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptions)  
2 -% IN_FOPEN_EDF: Open a BDF/EDF file (continuous recordings) 1 +function F = in_fread_edf(sFile, sfid, SamplesBounds, ChannelsRange)
  2 +% IN_FREAD_EDF: Read a block of recordings from a EDF/BDF file
3 % 3 %
4 -% USAGE: [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptions) 4 +% USAGE: F = in_fread_edf(sFile, sfid, SamplesBounds, ChannelsRange)
  5 +% F = in_fread_edf(sFile, sfid, SamplesBounds) : Read all channels
  6 +% F = in_fread_edf(sFile, sfid) : Read all channels, all the times
5 7
6 % @============================================================================= 8 % @=============================================================================
7 % This function is part of the Brainstorm software: 9 % This function is part of the Brainstorm software:
8 % https://neuroimage.usc.edu/brainstorm 10 % https://neuroimage.usc.edu/brainstorm
9 % 11 %
10 -% Copyright (c)2000-2020 University of Southern California & McGill University 12 +% Copyright (c)2000-2019 University of Southern California & McGill University
11 % This software is distributed under the terms of the GNU General Public License 13 % This software is distributed under the terms of the GNU General Public License
12 % as published by the Free Software Foundation. Further details on the GPLv3 14 % as published by the Free Software Foundation. Further details on the GPLv3
13 % license can be found at http://www.gnu.org/copyleft/gpl.html. 15 % license can be found at http://www.gnu.org/copyleft/gpl.html.
@@ -21,477 +23,203 @@ function [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptio @@ -21,477 +23,203 @@ function [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptio
21 % For more information type "brainstorm license" at command prompt. 23 % For more information type "brainstorm license" at command prompt.
22 % =============================================================================@ 24 % =============================================================================@
23 % 25 %
24 -% Authors: Francois Tadel, 2012-2019  
25 -  
26 -  
27 -% Parse inputs  
28 -if (nargin < 2) || isempty(ImportOptions)  
29 - ImportOptions = db_template('ImportOptions');  
30 -end  
31 -  
32 -  
33 -%% ===== READ HEADER =====  
34 -% Open file  
35 -fid = fopen(DataFile, 'r', 'ieee-le');  
36 -if (fid == -1)  
37 - error('Could not open file');  
38 -end  
39 -% Read all fields  
40 -hdr.version = fread(fid, [1 8], 'uint8=>char'); % Version of this data format ('0 ' for EDF, [255 'BIOSEMI'] for BDF)  
41 -hdr.patient_id = fread(fid, [1 80], '*char'); % Local patient identification  
42 -hdr.rec_id = fread(fid, [1 80], '*char'); % Local recording identification  
43 -hdr.startdate = fread(fid, [1 8], '*char'); % Startdate of recording (dd.mm.yy)  
44 -hdr.starttime = fread(fid, [1 8], '*char'); % Starttime of recording (hh.mm.ss)  
45 -hdr.hdrlen = str2double(fread(fid, [1 8], '*char')); % Number of bytes in header record  
46 -hdr.unknown1 = fread(fid, [1 44], '*char'); % Reserved ('24BIT' for BDF)  
47 -hdr.nrec = str2double(fread(fid, [1 8], '*char')); % Number of data records (-1 if unknown)  
48 -hdr.reclen = str2double(fread(fid, [1 8], '*char')); % Duration of a data record, in seconds  
49 -hdr.nsignal = str2double(fread(fid, [1 4], '*char')); % Number of signals in data record  
50 -% Check file integrity  
51 -if isnan(hdr.nsignal) || isempty(hdr.nsignal) || (hdr.nsignal ~= round(hdr.nsignal)) || (hdr.nsignal < 0)  
52 - error('File header is corrupted.');  
53 -end  
54 -% Read values for each nsignal  
55 -for i = 1:hdr.nsignal  
56 - hdr.signal(i).label = strtrim(fread(fid, [1 16], '*char'));  
57 -end  
58 -for i = 1:hdr.nsignal  
59 - hdr.signal(i).type = strtrim(fread(fid, [1 80], '*char'));  
60 -end  
61 -for i = 1:hdr.nsignal  
62 - hdr.signal(i).unit = strtrim(fread(fid, [1 8], '*char'));  
63 -end  
64 -for i = 1:hdr.nsignal  
65 - hdr.signal(i).physical_min = str2double(fread(fid, [1 8], '*char'));  
66 -end  
67 -for i = 1:hdr.nsignal  
68 - hdr.signal(i).physical_max = str2double(fread(fid, [1 8], '*char'));  
69 -end  
70 -for i = 1:hdr.nsignal  
71 - hdr.signal(i).digital_min = str2double(fread(fid, [1 8], '*char'));  
72 -end  
73 -for i = 1:hdr.nsignal  
74 - hdr.signal(i).digital_max = str2double(fread(fid, [1 8], '*char'));  
75 -end  
76 -for i = 1:hdr.nsignal  
77 - hdr.signal(i).filters = strtrim(fread(fid, [1 80], '*char'));  
78 -end  
79 -for i = 1:hdr.nsignal  
80 - hdr.signal(i).nsamples = str2num(fread(fid, [1 8], '*char'));  
81 -end  
82 -for i = 1:hdr.nsignal  
83 - hdr.signal(i).unknown2 = fread(fid, [1 32], '*char');  
84 -end  
85 -% Unknown record size, determine correct nrec  
86 -if (hdr.nrec == -1)  
87 - datapos = ftell(fid);  
88 - fseek(fid, 0, 'eof');  
89 - endpos = ftell(fid);  
90 - hdr.nrec = floor((endpos - datapos) / (sum([hdr.signal.nsamples]) * 2));  
91 -end  
92 -% Close file  
93 -fclose(fid);  
94 -  
95 -  
96 -%% ===== RECONSTRUCT INFO =====  
97 -% Individual signal gain  
98 -for i = 1:hdr.nsignal  
99 - % Interpet units  
100 - switch (hdr.signal(i).unit)  
101 - case 'mV', unit_gain = 1e3;  
102 - case {'uV', char([166 204 86])}, unit_gain = 1e6;  
103 - otherwise, unit_gain = 1;  
104 - end  
105 - % Check min/max values  
106 - if isempty(hdr.signal(i).digital_min) || isnan(hdr.signal(i).digital_min)  
107 - disp(['EDF> Warning: The digitial minimum is not set for channel "' hdr.signal(i).label '".']);  
108 - hdr.signal(i).digital_min = -2^15;  
109 - end  
110 - if isempty(hdr.signal(i).digital_max) || isnan(hdr.signal(i).digital_max)  
111 - disp(['EDF> Warning: The digitial maximum is not set for channel "' hdr.signal(i).label '".']);  
112 - hdr.signal(i).digital_max = -2^15;  
113 - end  
114 - if isempty(hdr.signal(i).physical_min) || isnan(hdr.signal(i).physical_min)  
115 - disp(['EDF> Warning: The physical minimum is not set for channel "' hdr.signal(i).label '".']);  
116 - hdr.signal(i).physical_min = hdr.signal(i).digital_min;  
117 - end  
118 - if isempty(hdr.signal(i).physical_max) || isnan(hdr.signal(i).physical_max)  
119 - disp(['EDF> Warning: The physical maximum is not set for channel "' hdr.signal(i).label '".']);  
120 - hdr.signal(i).physical_max = hdr.signal(i).digital_max;  
121 - end  
122 - if (hdr.signal(i).physical_min >= hdr.signal(i).physical_max)  
123 - disp(['EDF> Warning: Physical maximum larger than minimum for channel "' hdr.signal(i).label '".']);  
124 - hdr.signal(i).physical_min = hdr.signal(i).digital_min;  
125 - hdr.signal(i).physical_max = hdr.signal(i).digital_max;  
126 - end  
127 - % Calculate and save channel gain  
128 - hdr.signal(i).gain = unit_gain ./ (hdr.signal(i).physical_max - hdr.signal(i).physical_min) .* (hdr.signal(i).digital_max - hdr.signal(i).digital_min);  
129 - hdr.signal(i).offset = hdr.signal(i).physical_min ./ unit_gain - hdr.signal(i).digital_min ./ hdr.signal(i).gain;  
130 - % Error: The number of samples is not specified  
131 - if isempty(hdr.signal(i).nsamples)  
132 - % If it is not the first electrode: try to use the previous one  
133 - if (i > 1)  
134 - disp(['EDF> Warning: The number of samples is not specified for channel "' hdr.signal(i).label '".']);  
135 - hdr.signal(i).nsamples = hdr.signal(i-1).nsamples;  
136 - else  
137 - error(['The number of samples is not specified for channel "' hdr.signal(i).label '".']);  
138 - end  
139 - end  
140 - hdr.signal(i).sfreq = hdr.signal(i).nsamples ./ hdr.reclen;  
141 -end  
142 -% Find annotations channel  
143 -iAnnotChans = find(strcmpi({hdr.signal.label}, 'EDF Annotations')); % Mutliple "EDF Annotation" channels allowed in EDF+  
144 -iStatusChan = find(strcmpi({hdr.signal.label}, 'Status'), 1); % Only one "Status" channel allowed in BDF  
145 -iOtherChan = setdiff(1:hdr.nsignal, [iAnnotChans iStatusChan]);  
146 -% % Remove channels with lower sampling rates  
147 -% iIgnoreChan = find([hdr.signal(iOtherChan).sfreq] < max([hdr.signal(iOtherChan).sfreq])); % Ignore all the channels with lower sampling rate  
148 -% if ~isempty(iIgnoreChan)  
149 -% iOtherChan = setdiff(iOtherChan, iIgnoreChan);  
150 -% end  
151 -% Get all the other channels  
152 -if isempty(iOtherChan)  
153 - error('This file does not contain any data channel.');  
154 -end  
155 -% Read events preferencially from the EDF Annotations track  
156 -if ~isempty(iAnnotChans)  
157 - iEvtChans = iAnnotChans;  
158 -elseif ~isempty(iStatusChan)  
159 - iEvtChans = iStatusChan; 26 +% Authors: Francois Tadel, 2012-2017
  27 +
  28 +%% ===== PARSE INPUTS =====
  29 +nChannels = sFile.header.nsignal;
  30 +iChanAnnot = find(strcmpi({sFile.header.signal.label}, 'EDF Annotations'));
  31 +iBadChan = find(sFile.channelflag == -1);
  32 +iChanSignal = setdiff(1:nChannels, iChanAnnot);
  33 +iChanWrongRate = find([sFile.header.signal(iChanSignal).sfreq] ~= max([sFile.header.signal(setdiff(iChanSignal,iBadChan)).sfreq]));
  34 +iChanSkip = union(iChanAnnot, iChanWrongRate);
  35 +if (nargin < 4) || isempty(ChannelsRange)
  36 + ChannelsRange = [1, nChannels];
  37 +end
  38 +if (nargin < 3) || isempty(SamplesBounds)
  39 + SamplesBounds = [0, sFile.header.nrec * sFile.header.signal(ChannelsRange(1)).nsamples - 1];
  40 +end
  41 +nTimes = sFile.header.reclen * sFile.header.signal(ChannelsRange(1)).sfreq;
  42 +iTimes = SamplesBounds(1):SamplesBounds(2);
  43 +% Block of times/channels to extract
  44 +nReadChannels = double(ChannelsRange(2) - ChannelsRange(1) + 1);
  45 +% Read annotations instead of real data ?
  46 +isAnnotOnly = ~isempty(iChanAnnot) && (ChannelsRange(1) == ChannelsRange(2)) && ismember(ChannelsRange(1), iChanAnnot);
  47 +isBdfStatus = strcmpi(sFile.format, 'EEG-BDF') && (ChannelsRange(1) == ChannelsRange(2)) && strcmpi(sFile.header.signal(ChannelsRange(1)).label, 'Status');
  48 +% Data channels to read
  49 +if isAnnotOnly
  50 + % Read only on annotation channel
  51 + iChanF = 1;
160 else 52 else
161 - iEvtChans = [];  
162 -end  
163 -% % Detect channels with inconsistent sampling frenquency  
164 -% iErrChan = find([hdr.signal(iOtherChan).sfreq] ~= hdr.signal(iOtherChan(1)).sfreq);  
165 -% iErrChan = setdiff(iErrChan, iAnnotChans);  
166 -% if ~isempty(iErrChan)  
167 -% error('Files with mixed sampling rates are not supported yet.');  
168 -% end  
169 -% Detect interrupted signals (time non-linear)  
170 -hdr.interrupted = ischar(hdr.unknown1) && (length(hdr.unknown1) >= 5) && isequal(hdr.unknown1(1:5), 'EDF+D');  
171 -if hdr.interrupted  
172 - if ImportOptions.DisplayMessages  
173 - [res, isCancel] = java_dialog('question', ...  
174 - ['Interrupted EDF file ("EDF+D") detected. It is recommended to convert it' 10 ...  
175 - 'to a continuous ("EDF+C") file first. Do you want to continue reading this' 10 ...  
176 - 'file as continuous and attempt to fix the timing of event markers?' 10 ...  
177 - 'NOTE: This may not work as intended, use at your own risk!']);  
178 - hdr.fixinterrupted = ~isCancel && strcmpi(res, 'yes'); 53 + % Remove all the annotation channels from the list of channels to read
  54 + iChanF = setdiff(ChannelsRange(1):ChannelsRange(2), iChanSkip) - ChannelsRange(1) + 1;
  55 +% if any(diff(iChanF) ~= 1)
  56 +% error('All the data channels to read from the file must be contiguous (EDF Annotation channels must be at the end of the list).');
  57 +% end
  58 + if any(diff(iChanF) ~= 1)
  59 + iChanLast = find(diff(iChanF) ~= 1, 1);
  60 + iChanF = iChanF(1:iChanLast);
179 else 61 else
180 - hdr.fixinterrupted = 1; 62 + iChanLast = length(iChanF);
181 end 63 end
182 - if ~hdr.fixinterrupted  
183 - warning(['Interrupted EDF file ("EDF+D"): requires conversion to "EDF+C". ' 10 ...  
184 - 'Brainstorm will read this file as a continuous file ("EDF+C"), the timing of the samples after the first discontinuity will be wrong.' 10 ...  
185 - 'This may not cause any major problem unless there are time markers in the file, they might be inaccurate in all the segments >= 2.']); 64 + ChannelsRange = [iChanF(1), iChanF(iChanLast)] + ChannelsRange(1) - 1;
  65 +end
  66 +% Cannot read channels with different sampling rates at the same time
  67 +if (ChannelsRange(1) ~= ChannelsRange(2))
  68 + allFreq = [sFile.header.signal(ChannelsRange(1):ChannelsRange(2)).nsamples];
  69 + if any(allFreq ~= allFreq(1))
  70 + error(['Cannot read channels with mixed sampling rates at the same time.' 10 ...
  71 + 'Mark as bad channels with different sampling rates than EEG.' 10 ...
  72 + '(right-click on data file > Good/bad channels > Edit good/bad channels)']);
186 end 73 end
187 -else  
188 - hdr.fixinterrupted = 0;  
189 end 74 end
190 75
191 76
192 -%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====  
193 -% Initialize returned file structure  
194 -sFile = db_template('sfile');  
195 -% Add information read from header  
196 -sFile.byteorder = 'l';  
197 -sFile.filename = DataFile;  
198 -if (uint8(hdr.version(1)) == uint8(255))  
199 - sFile.format = 'EEG-BDF';  
200 - sFile.device = 'BDF'; 77 +%% ===== READ ALL NEEDED EPOCHS =====
  78 +% Detect which epochs are necessary for the range of data selected
  79 +epochRange = floor(SamplesBounds ./ nTimes);
  80 +epochsToRead = epochRange(1) : epochRange(2);
  81 +% Initialize block of data to read
  82 +if isAnnotOnly
  83 + F = zeros(nReadChannels, 2 * length(iTimes));
201 else 84 else
202 - sFile.format = 'EEG-EDF';  
203 - sFile.device = 'EDF'; 85 + F = zeros(nReadChannels, length(iTimes));
  86 +end
  87 +% Marker that we increment when we add data to F
  88 +iF = 1;
  89 +% Read all the needed epochs
  90 +for i = 1:length(epochsToRead)
  91 + % Find the samples to read from this epoch
  92 + BoundsEpoch = nTimes * epochsToRead(i) + [0, nTimes-1];
  93 + BoundsRead = [max(BoundsEpoch(1), iTimes(1)), ...
  94 + min(BoundsEpoch(2), iTimes(end))];
  95 + iTimeRead = BoundsRead(1):BoundsRead(2);
  96 + % Convert this samples into indices in this very epoch
  97 + iTimeRead = iTimeRead - nTimes * epochsToRead(i);
  98 + % New indices to read
  99 + if isAnnotOnly
  100 + iNewF = iF:(iF + 2*length(iTimeRead) - 1);
  101 + else
  102 + iNewF = iF:(iF + length(iTimeRead) - 1);
  103 + end
  104 + % Read epoch (full or partial)
  105 + F(iChanF,iNewF) = edf_read_epoch(sFile, sfid, epochsToRead(i), iTimeRead, ChannelsRange, isAnnotOnly, isBdfStatus);
  106 + % Increment marker
  107 + iF = iF + length(iTimeRead);
204 end 108 end
205 -sFile.header = hdr;  
206 -% Comment: short filename  
207 -[tmp__, sFile.comment, tmp__] = bst_fileparts(DataFile);  
208 -% No info on bad channels  
209 -sFile.channelflag = ones(hdr.nsignal,1);  
210 -% Acquisition date  
211 -sFile.acq_date = str_date(hdr.startdate);  
212 109
213 110
214 111
215 -%% ===== PROCESS CHANNEL NAMES/TYPES =====  
216 -% Try to split the channel names in "TYPE NAME"  
217 -SplitType = repmat({''}, 1, hdr.nsignal);  
218 -SplitName = repmat({''}, 1, hdr.nsignal);  
219 -for i = 1:hdr.nsignal  
220 - % Removing trailing dots (eg. "Fc5." instead of "FC5", as in: https://www.physionet.org/pn4/eegmmidb/)  
221 - if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)  
222 - hdr.signal(i).label(end) = [];  
223 - if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)  
224 - hdr.signal(i).label(end) = [];  
225 - if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)  
226 - hdr.signal(i).label(end) = [];  
227 - end  
228 - end  
229 - end  
230 - % Remove extra spaces  
231 - signalLabel = strrep(hdr.signal(i).label, ' - ', '-');  
232 - % Find space chars (label format "Type Name")  
233 - iSpace = find(signalLabel == ' ');  
234 - % Only if there is one space only  
235 - if (length(iSpace) == 1) && (iSpace >= 3)  
236 - SplitName{i} = signalLabel(iSpace+1:end);  
237 - SplitType{i} = signalLabel(1:iSpace-1);  
238 - % Accept also 2 spaces  
239 - elseif (length(iSpace) == 2) && (iSpace(1) >= 3)  
240 - SplitName{i} = strrep(signalLabel(iSpace(1)+1:end), ' ', '_');  
241 - SplitType{i} = signalLabel(1:iSpace(1)-1);  
242 - end  
243 -end  
244 -% Remove the classification if it makes some names non unique  
245 -uniqueNames = unique(SplitName);  
246 -for i = 1:length(uniqueNames)  
247 - if ~isempty(uniqueNames{i})  
248 - iName = find(strcmpi(SplitName, uniqueNames{i}));  
249 - if (length(iName) > 1)  
250 - [SplitName{iName}] = deal('');  
251 - [SplitType{iName}] = deal('');  
252 - end  
253 - end  
254 end 112 end
255 113
256 114
257 -%% ===== CREATE EMPTY CHANNEL FILE =====  
258 -ChannelMat = db_template('channelmat');  
259 -ChannelMat.Comment = [sFile.device ' channels'];  
260 -ChannelMat.Channel = repmat(db_template('channeldesc'), [1, hdr.nsignal]);  
261 -chRef = {};  
262 -% For each channel  
263 -for i = 1:hdr.nsignal  
264 - % If is the annotation channel  
265 - if ~isempty(iAnnotChans) && ismember(i, iAnnotChans)  
266 - ChannelMat.Channel(i).Type = 'EDF';  
267 - ChannelMat.Channel(i).Name = 'Annotations';  
268 - elseif ~isempty(iStatusChan) && (i == iStatusChan)  
269 - ChannelMat.Channel(i).Type = 'BDF';  
270 - ChannelMat.Channel(i).Name = 'Status';  
271 - % Regular channels  
272 - else  
273 - % If there is a pair name/type already detected  
274 - if ~isempty(SplitName{i}) && ~isempty(SplitType{i})  
275 - ChannelMat.Channel(i).Name = SplitName{i};  
276 - ChannelMat.Channel(i).Type = SplitType{i}; 115 +
  116 +%% ===== READ ONE EPOCH =====
  117 +function F = edf_read_epoch(sFile, sfid, iEpoch, iTimes, ChannelsRange, isAnnotOnly, isBdfStatus)
  118 + % ===== COMPUTE OFFSETS =====
  119 + nTimes = sFile.header.reclen * [sFile.header.signal.sfreq];
  120 + nReadTimes = length(iTimes);
  121 + nReadChannels = double(ChannelsRange(2) - ChannelsRange(1) + 1);
  122 + iChannels = ChannelsRange(1):ChannelsRange(2);
  123 + % Check that all the channels selected have the same freq rate
  124 + if any(nTimes(iChannels) ~= nTimes(iChannels(1)))
  125 + error('Cannot read at the same signals with different sampling frequency.');
  126 + end
  127 + % Size of one value
  128 + if strcmpi(sFile.format, 'EEG-BDF')
  129 + % BDF: int24 => 3 bytes
  130 + bytesPerVal = 3;
  131 + % Reading status or regular channel
  132 + if isBdfStatus
  133 + dataClass = 'ubit24';
277 else 134 else
278 - % Channel name  
279 - ChannelMat.Channel(i).Name = hdr.signal(i).label(hdr.signal(i).label ~= ' ');  
280 - % Channel type  
281 - if ~isempty(hdr.signal(i).type)  
282 - if (length(hdr.signal(i).type) == 3)  
283 - ChannelMat.Channel(i).Type = hdr.signal(i).type(hdr.signal(i).type ~= ' ');  
284 - elseif isequal(hdr.signal(i).type, 'Active Electrode') || isequal(hdr.signal(i).type, 'AgAgCl electrode')  
285 - ChannelMat.Channel(i).Type = 'EEG';  
286 - else  
287 - ChannelMat.Channel(i).Type = 'Misc';  
288 - end  
289 - else  
290 - ChannelMat.Channel(i).Type = 'EEG';  
291 - end  
292 - end  
293 - % Extract reference name (at the end of the channel name, separated with a "-", eg. "-REF")  
294 - iDash = find(ChannelMat.Channel(i).Name == '-');  
295 - if ~isempty(iDash) && (iDash(end) < length(ChannelMat.Channel(i).Name))  
296 - chRef{end+1} = ChannelMat.Channel(i).Name(iDash(end):end); 135 + dataClass = 'bit24';
297 end 136 end
  137 + else
  138 + % EDF: int16 => 2 bytes
  139 + bytesPerVal = 2;
  140 + dataClass = 'int16';
  141 + isBdfStatus = 0;
298 end 142 end
299 - ChannelMat.Channel(i).Loc = [0; 0; 0];  
300 - ChannelMat.Channel(i).Orient = [];  
301 - ChannelMat.Channel(i).Weight = 1;  
302 - % ChannelMat.Channel(i).Comment = hdr.signal(i).type;  
303 -end  
304 -  
305 -% If the same reference is indicated for all the channels: remove it  
306 -if (length(chRef) >= 2)  
307 - % Get the shortest reference tag  
308 - lenRef = cellfun(@length, chRef);  
309 - minLen = min(lenRef);  
310 - % Check if all the ref names are equal (up to the max length - some might be cut because the channel name is too long)  
311 - if all(cellfun(@(c)strcmpi(c(1:minLen), chRef{1}(1:minLen)), chRef))  
312 - % Remove the reference tag from all the channel names  
313 - for i = 1:length(ChannelMat.Channel)  
314 - ChannelMat.Channel(i).Name = strrep(ChannelMat.Channel(i).Name, chRef{1}, '');  
315 - ChannelMat.Channel(i).Name = strrep(ChannelMat.Channel(i).Name, chRef{1}(1:minLen), ''); 143 + % Offset of the beginning of the recordings in the file
  144 + offsetHeader = round(sFile.header.hdrlen);
  145 + % Offset of epoch
  146 + offsetEpoch = round(iEpoch * sum(nTimes) * bytesPerVal);
  147 + % Channel offset
  148 + offsetChannel = round(sum(nTimes(1:ChannelsRange(1)-1)) * bytesPerVal);
  149 + % Time offset at the beginning and end of each channel block
  150 + offsetTimeStart = round(iTimes(1) * bytesPerVal);
  151 + offsetTimeEnd = round((nTimes(ChannelsRange(1)) - iTimes(end) - 1) * bytesPerVal);
  152 + % ALL THE "ROUND" CALLS WERE ADDED AFTER DISCOVERING THAT THERE WERE SOMETIMES ROUNDING ERRORS IN THE MULTIPLICATIONS
  153 +
  154 + % Where to start reading in the file ?
  155 + % => After the header, the number of skipped epochs, channels and time samples
  156 + offsetStart = offsetHeader + offsetEpoch + offsetChannel + offsetTimeStart;
  157 + % Number of time samples to skip after each channel
  158 + offsetSkip = offsetTimeStart + offsetTimeEnd;
  159 +
  160 +
  161 + % ===== READ DATA BLOCK =====
  162 + % Position file at the beginning of the trial
  163 + fseek(sfid, offsetStart, 'bof');
  164 + % Read annotation data (char)
  165 + if isAnnotOnly
  166 + dataClass = 'char';
  167 + nReadTimes = bytesPerVal * nReadTimes; % 1 byte instead of 2
  168 + end
  169 + % Read trial data
  170 + % => WARNING: CALL TO FREAD WITH SKIP=0 DOES NOT WORK PROPERLY
  171 + if (offsetSkip == 0)
  172 + F = fread(sfid, [nReadTimes, nReadChannels], dataClass)';
  173 + elseif (bytesPerVal == 2)
  174 + precision = sprintf('%d*%s', nReadTimes, dataClass);
  175 + F = fread(sfid, [nReadTimes, nReadChannels], precision, offsetSkip)';
  176 + % => WARNING: READING USING ubit24 SOMETIMES DOESNT WORK => DOING IT MANUALLY
  177 + elseif (bytesPerVal == 3)
  178 + % Reading each bit independently
  179 + precision = sprintf('%d*%s', 3*nReadTimes, 'uint8');
  180 + F = fread(sfid, [3*nReadTimes, nReadChannels], precision, offsetSkip)';
  181 + % Grouping the 3 bits together
  182 + F = F(:,1:3:end) + F(:,2:3:end)*256 + F(:,3:3:end)*256*256;
  183 + % 2-Complement (negative value indicated by most significant bit)
  184 + if strcmpi(dataClass, 'bit24')
  185 + iNeg = (F >= 256*256*128);
  186 + F(iNeg) = F(iNeg) - 256*256*256;
316 end 187 end
317 end 188 end
318 -end  
319 -  
320 -% If there are only "Misc" and no "EEG" channels: rename to "EEG"  
321 -iMisc = find(strcmpi({ChannelMat.Channel.Type}, 'Misc'));  
322 -iEeg = find(strcmpi({ChannelMat.Channel.Type}, 'EEG'));  
323 -if ~isempty(iMisc) && isempty(iEeg)  
324 - [ChannelMat.Channel(iMisc).Type] = deal('EEG');  
325 - iEeg = iMisc;  
326 -end  
327 -  
328 -  
329 -%% ===== DETECT MULTIPLE SAMPLING RATES =====  
330 -% Use the first "EEG" channel as the reference sampling rate (or the first channel if no "EEG" channels available)  
331 -if ~isempty(iEeg) && ismember(iEeg(1), iOtherChan)  
332 - iChanFreqRef = iEeg(1);  
333 -else  
334 - iChanFreqRef = iOtherChan(1);  
335 -end  
336 -% Mark as bad channels with sampling rates different from EEG  
337 -iChanWrongRate = find([sFile.header.signal.sfreq] ~= sFile.header.signal(iChanFreqRef).sfreq);  
338 -iChanWrongRate = intersect(iChanWrongRate, iOtherChan);  
339 -if ~isempty(iChanWrongRate)  
340 - sFile.channelflag(iChanWrongRate) = -1;  
341 -end  
342 -  
343 -% Consider that the sampling rate of the file is the sampling rate of the first signal  
344 -sFile.prop.sfreq = hdr.signal(iChanFreqRef).sfreq;  
345 -sFile.prop.times = [0, hdr.signal(iChanFreqRef).nsamples * hdr.nrec - 1] ./ sFile.prop.sfreq;  
346 -sFile.prop.nAvg = 1;  
347 -  
348 -  
349 - 189 + % Check that data block was fully read
  190 + if (numel(F) < nReadTimes * nReadChannels)
  191 + % Number of full time samples that were read
  192 + nTimeTrunc = max(0, floor(numel(F) / nReadChannels) - 1);
  193 + % Error message
  194 + disp(sprintf('EDF> ERROR: File is truncated (%d time samples were read instead of %d). Padding with zeros...', nTimeTrunc, nReadTimes));
  195 + % Pad data with zeros
  196 + Ftmp = zeros(nReadTimes, nReadChannels);
  197 + F = F';
  198 + Ftmp(1:numel(F)) = F(:);
  199 + F = Ftmp';
  200 + end
350 201
351 -%% ===== READ EDF ANNOTATION CHANNEL =====  
352 -if ~isempty(iEvtChans) % && ~isequal(ImportOptions.EventsMode, 'ignore')  
353 - % Set reading options  
354 - ImportOptions.ImportMode = 'Time';  
355 - ImportOptions.UseSsp = 0;  
356 - ImportOptions.UseCtfComp = 0;  
357 - % Read EDF annotations  
358 - if strcmpi(sFile.format, 'EEG-EDF')  
359 - evtList = {};  
360 - % In EDF+, the first annotation channel has epoch time stamps (EDF  
361 - % calls epochs records). So read all annotation channels per epoch.  
362 - for irec = 1:hdr.nrec  
363 - for ichan = 1:length(iEvtChans)  
364 - bst_progress('text', sprintf('Reading annotations... [%d%%]', round((ichan + (irec-1)*length(iEvtChans))/length(iEvtChans)/hdr.nrec*100)));  
365 - % Sample indices for the current epoch (=record)  
366 - SampleBounds = [irec-1,irec] * sFile.header.signal(iEvtChans(ichan)).nsamples - [0,1];  
367 - % Read record  
368 - F = in_fread(sFile, ChannelMat, 1, SampleBounds, iEvtChans(ichan), ImportOptions);  
369 - % Find all the TALs (Time-stamped Annotations Lists): separated with [20][0]  
370 - % Possible configurations:  
371 - % Onset[21]Duration[20]Annot1[20]Annot2[20]..AnnotN[20][0]  
372 - % Onset[20]Annot1[20]Annot2[20]..AnnotN[20][0]  
373 - % Onset[20]Annot1[20][0]  
374 - iSeparator = [-1, find((F(1:end-1) == 20) & (F(2:end) == 0))];  
375 - % Loop on all the TALs  
376 - for iAnnot = 1:length(iSeparator)-1  
377 - % Get annotation  
378 - strAnnot = char(F(iSeparator(iAnnot)+2:iSeparator(iAnnot+1)-1));  
379 - % Split in blocks with [20]  
380 - splitAnnot = str_split(strAnnot, 20);  
381 - % The first TAL in a record should be indicating the timing of the first sample of the record  
382 - if (iAnnot == 1) && (length(splitAnnot) == 1) && ~any(splitAnnot{1} == 21)  
383 - % Ignore if this is not the first channel  
384 - if (ichan > 1)  
385 - continue;  
386 - end  
387 - % Get record time stamp  
388 - t0_rec = str2double(splitAnnot{1});  
389 - if isempty(t0_rec) || isnan(t0_rec)  
390 - continue;  
391 - end  
392 - if (irec == 1)  
393 - t0_file = t0_rec;  
394 - % Find discontinuities larger than 1 sample  
395 - elseif abs(t0_rec - prev_rec - (irec - prev_irec) * hdr.reclen) > (1 / sFile.prop.sfreq)  
396 - % Brainstorm fills partial/interrupted records with zeros  
397 - bstTime = prev_rec + hdr.reclen;  
398 - timeDiff = bstTime - t0_rec;  
399 - % If we want to fix timing, apply skip to initial timestamp  
400 - if hdr.fixinterrupted  
401 - t0_file = t0_file - timeDiff;  
402 - end  
403 - % Warn user of discontinuity  
404 - if timeDiff > 0  
405 - expectMsg = 'blank data';  
406 - else  
407 - expectMsg = 'skipped data';  
408 - end  
409 - startTime = min(t0_rec - t0_file - [0, timeDiff]); % before and after t0_file adjustment  
410 - endTime = max(t0_rec - t0_file - [0, timeDiff]);  
411 - fprintf('WARNING: Found discontinuity between %.3fs and %.3fs, expect %s in between.\n', startTime, endTime, expectMsg);  
412 - % Create event for users information  
413 - if timeDiff < 0  
414 - endTime = startTime; % no extent in this case, there is skipped time.  
415 - end  
416 - evtList(end+1,:) = {'EDF+D Discontinuity', [startTime; endTime]};  
417 - end  
418 - prev_rec = t0_rec;  
419 - prev_irec = irec;  
420 - % Regular TAL: indicating a marker  
421 - else  
422 - % Split time in onset/duration  
423 - splitTime = str_split(splitAnnot{1}, 21);  
424 - % Get time  
425 - t = str2double(splitTime{1});  
426 - if isempty(t) || isnan(t)  
427 - continue;  
428 - end  
429 - % Get duration  
430 - if (length(splitTime) > 1)  
431 - duration = str2double(splitTime{2});  
432 - % Exclude 1-sample long events  
433 - if isempty(duration) || isnan(duration) || (round(duration .* sFile.prop.sfreq) <= 1)  
434 - duration = 0;  
435 - end  
436 - else  
437 - duration = 0;  
438 - end  
439 - % Unnamed event  
440 - if (length(splitAnnot) == 1) || isempty(splitAnnot{2})  
441 - splitAnnot{2} = 'Unnamed';  
442 - end  
443 - % Create one event for each label in the TAL  
444 - for iLabel = 2:length(splitAnnot)  
445 - evtList(end+1,:) = {splitAnnot{iLabel}, (t-t0_file) + [0;duration]};  
446 - end  
447 - end  
448 - end  
449 - end  
450 - end  
451 -  
452 - % If there are events: create a create an events structure  
453 - if ~isempty(evtList)  
454 - % Initialize events list  
455 - sFile.events = repmat(db_template('event'), 0);  
456 - % Events list  
457 - [uniqueEvt, iUnique] = unique(evtList(:,1));  
458 - uniqueEvt = evtList(sort(iUnique),1);  
459 - % Build events list  
460 - for iEvt = 1:length(uniqueEvt)  
461 - % Find all the occurrences of this event  
462 - iOcc = find(strcmpi(uniqueEvt{iEvt}, evtList(:,1)));  
463 - % Concatenate all times  
464 - t = [evtList{iOcc,2}];  
465 - % If second row is equal to the first one (no extended events): delete it  
466 - if all(t(1,:) == t(2,:))  
467 - t = t(1,:);  
468 - end  
469 - % Set event  
470 - sFile.events(iEvt).label = strtrim(uniqueEvt{iEvt});  
471 - sFile.events(iEvt).times = round(t .* sFile.prop.sfreq) ./ sFile.prop.sfreq;  
472 - sFile.events(iEvt).epochs = 1 + 0*t(1,:);  
473 - sFile.events(iEvt).select = 1;  
474 - sFile.events(iEvt).channels = cell(1, size(sFile.events(iEvt).times, 2));  
475 - sFile.events(iEvt).notes = cell(1, size(sFile.events(iEvt).times, 2));  
476 - end  
477 - end  
478 -  
479 - % BDF Status line  
480 - elseif strcmpi(sFile.format, 'EEG-BDF') && ~strcmpi(ImportOptions.EventsTrackMode, 'ignore')  
481 - % Ask how to read the events  
482 - [events, ImportOptions.EventsTrackMode] = process_evt_read('Compute', sFile, ChannelMat, ChannelMat.Channel(iEvtChans).Name, ImportOptions.EventsTrackMode);  
483 - if isequal(events, -1)  
484 - sFile = [];  
485 - ChannelMat = [];  
486 - return;  
487 - end  
488 - % Report the events in the file structure  
489 - sFile.events = events;  
490 - % Remove the 'Status: ' string in front of the events  
491 - for i = 1:length(sFile.events)  
492 - sFile.events(i).label = strrep(sFile.events(i).label, 'Status: ', '');  
493 - end  
494 - % Group events by time  
495 - % sFile.events = process_evt_grouptime('Compute', sFile.events); 202 + % Processing for BDF status file
  203 + if isBdfStatus
  204 + % Mask to keep only the first 15 bits (Triggers bits)
  205 + % Bit 16 : High when new Epoch is started
  206 + % Bit 17-19 : Speed bits 0 1 2
  207 + % Bit 20 : High when CMS is within range
  208 + % Bit 21 : Speed bit 3
  209 + % Bit 22 : High when battery is low
  210 + % Bit 23 : High if ActiveTwo MK2
  211 + F = bitand(F, bin2dec('000000000111111111111111'));
  212 + % Processing for real data
  213 + elseif ~isAnnotOnly
  214 + % Convert to double
  215 + F = double(F);
  216 + % Apply gains
  217 + F = bst_bsxfun(@rdivide, F, [sFile.header.signal(iChannels).gain]');
  218 +% IN THEORY: THIS OFFSET SECTION IS USEFUL, BUT IN PRACTICE IT LOOKS LIKE THE VALUES IN ALL THE FILES ARE CENTERED ON ZERO
  219 +% % Add offset
  220 +% if isfield(sFile.header.signal, 'offset') && ~isempty(sFile.header.signal(1).offset)
  221 +% % ...
  222 +% end
496 end 223 end
497 end 224 end
  225 +