|
|
|
function [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptions)
|
|
|
|
% IN_FOPEN_EDF: Open a BDF/EDF file (continuous recordings)
|
|
|
|
function F = in_fread_edf(sFile, sfid, SamplesBounds, ChannelsRange)
|
|
|
|
% IN_FREAD_EDF: Read a block of recordings from a EDF/BDF file
|
|
|
|
%
|
|
|
|
% USAGE: [sFile, ChannelMat, ImportOptions] = in_fopen_edf(DataFile, ImportOptions)
|
|
|
|
% USAGE: F = in_fread_edf(sFile, sfid, SamplesBounds, ChannelsRange)
|
|
|
|
% F = in_fread_edf(sFile, sfid, SamplesBounds) : Read all channels
|
|
|
|
% F = in_fread_edf(sFile, sfid) : Read all channels, all the times
|
|
|
|
|
|
|
|
% @=============================================================================
|
|
|
|
% This function is part of the Brainstorm software:
|
|
|
|
% https://neuroimage.usc.edu/brainstorm
|
|
|
|
%
|
|
|
|
% Copyright (c)2000-2020 University of Southern California & McGill University
|
|
|
|
% Copyright (c)2000-2019 University of Southern California & McGill University
|
|
|
|
% This software is distributed under the terms of the GNU General Public License
|
|
|
|
% as published by the Free Software Foundation. Further details on the GPLv3
|
|
|
|
% 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 |
|
|
|
% For more information type "brainstorm license" at command prompt.
|
|
|
|
% =============================================================================@
|
|
|
|
%
|
|
|
|
% Authors: Francois Tadel, 2012-2019
|
|
|
|
|
|
|
|
|
|
|
|
% Parse inputs
|
|
|
|
if (nargin < 2) || isempty(ImportOptions)
|
|
|
|
ImportOptions = db_template('ImportOptions');
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
%% ===== READ HEADER =====
|
|
|
|
% Open file
|
|
|
|
fid = fopen(DataFile, 'r', 'ieee-le');
|
|
|
|
if (fid == -1)
|
|
|
|
error('Could not open file');
|
|
|
|
end
|
|
|
|
% Read all fields
|
|
|
|
hdr.version = fread(fid, [1 8], 'uint8=>char'); % Version of this data format ('0 ' for EDF, [255 'BIOSEMI'] for BDF)
|
|
|
|
hdr.patient_id = fread(fid, [1 80], '*char'); % Local patient identification
|
|
|
|
hdr.rec_id = fread(fid, [1 80], '*char'); % Local recording identification
|
|
|
|
hdr.startdate = fread(fid, [1 8], '*char'); % Startdate of recording (dd.mm.yy)
|
|
|
|
hdr.starttime = fread(fid, [1 8], '*char'); % Starttime of recording (hh.mm.ss)
|
|
|
|
hdr.hdrlen = str2double(fread(fid, [1 8], '*char')); % Number of bytes in header record
|
|
|
|
hdr.unknown1 = fread(fid, [1 44], '*char'); % Reserved ('24BIT' for BDF)
|
|
|
|
hdr.nrec = str2double(fread(fid, [1 8], '*char')); % Number of data records (-1 if unknown)
|
|
|
|
hdr.reclen = str2double(fread(fid, [1 8], '*char')); % Duration of a data record, in seconds
|
|
|
|
hdr.nsignal = str2double(fread(fid, [1 4], '*char')); % Number of signals in data record
|
|
|
|
% Check file integrity
|
|
|
|
if isnan(hdr.nsignal) || isempty(hdr.nsignal) || (hdr.nsignal ~= round(hdr.nsignal)) || (hdr.nsignal < 0)
|
|
|
|
error('File header is corrupted.');
|
|
|
|
end
|
|
|
|
% Read values for each nsignal
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).label = strtrim(fread(fid, [1 16], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).type = strtrim(fread(fid, [1 80], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).unit = strtrim(fread(fid, [1 8], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).physical_min = str2double(fread(fid, [1 8], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).physical_max = str2double(fread(fid, [1 8], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).digital_min = str2double(fread(fid, [1 8], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).digital_max = str2double(fread(fid, [1 8], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).filters = strtrim(fread(fid, [1 80], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).nsamples = str2num(fread(fid, [1 8], '*char'));
|
|
|
|
end
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
hdr.signal(i).unknown2 = fread(fid, [1 32], '*char');
|
|
|
|
end
|
|
|
|
% Unknown record size, determine correct nrec
|
|
|
|
if (hdr.nrec == -1)
|
|
|
|
datapos = ftell(fid);
|
|
|
|
fseek(fid, 0, 'eof');
|
|
|
|
endpos = ftell(fid);
|
|
|
|
hdr.nrec = floor((endpos - datapos) / (sum([hdr.signal.nsamples]) * 2));
|
|
|
|
end
|
|
|
|
% Close file
|
|
|
|
fclose(fid);
|
|
|
|
|
|
|
|
|
|
|
|
%% ===== RECONSTRUCT INFO =====
|
|
|
|
% Individual signal gain
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
% Interpet units
|
|
|
|
switch (hdr.signal(i).unit)
|
|
|
|
case 'mV', unit_gain = 1e3;
|
|
|
|
case {'uV', char([166 204 86])}, unit_gain = 1e6;
|
|
|
|
otherwise, unit_gain = 1;
|
|
|
|
end
|
|
|
|
% Check min/max values
|
|
|
|
if isempty(hdr.signal(i).digital_min) || isnan(hdr.signal(i).digital_min)
|
|
|
|
disp(['EDF> Warning: The digitial minimum is not set for channel "' hdr.signal(i).label '".']);
|
|
|
|
hdr.signal(i).digital_min = -2^15;
|
|
|
|
end
|
|
|
|
if isempty(hdr.signal(i).digital_max) || isnan(hdr.signal(i).digital_max)
|
|
|
|
disp(['EDF> Warning: The digitial maximum is not set for channel "' hdr.signal(i).label '".']);
|
|
|
|
hdr.signal(i).digital_max = -2^15;
|
|
|
|
end
|
|
|
|
if isempty(hdr.signal(i).physical_min) || isnan(hdr.signal(i).physical_min)
|
|
|
|
disp(['EDF> Warning: The physical minimum is not set for channel "' hdr.signal(i).label '".']);
|
|
|
|
hdr.signal(i).physical_min = hdr.signal(i).digital_min;
|
|
|
|
end
|
|
|
|
if isempty(hdr.signal(i).physical_max) || isnan(hdr.signal(i).physical_max)
|
|
|
|
disp(['EDF> Warning: The physical maximum is not set for channel "' hdr.signal(i).label '".']);
|
|
|
|
hdr.signal(i).physical_max = hdr.signal(i).digital_max;
|
|
|
|
end
|
|
|
|
if (hdr.signal(i).physical_min >= hdr.signal(i).physical_max)
|
|
|
|
disp(['EDF> Warning: Physical maximum larger than minimum for channel "' hdr.signal(i).label '".']);
|
|
|
|
hdr.signal(i).physical_min = hdr.signal(i).digital_min;
|
|
|
|
hdr.signal(i).physical_max = hdr.signal(i).digital_max;
|
|
|
|
end
|
|
|
|
% Calculate and save channel gain
|
|
|
|
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);
|
|
|
|
hdr.signal(i).offset = hdr.signal(i).physical_min ./ unit_gain - hdr.signal(i).digital_min ./ hdr.signal(i).gain;
|
|
|
|
% Error: The number of samples is not specified
|
|
|
|
if isempty(hdr.signal(i).nsamples)
|
|
|
|
% If it is not the first electrode: try to use the previous one
|
|
|
|
if (i > 1)
|
|
|
|
disp(['EDF> Warning: The number of samples is not specified for channel "' hdr.signal(i).label '".']);
|
|
|
|
hdr.signal(i).nsamples = hdr.signal(i-1).nsamples;
|
|
|
|
else
|
|
|
|
error(['The number of samples is not specified for channel "' hdr.signal(i).label '".']);
|
|
|
|
end
|
|
|
|
end
|
|
|
|
hdr.signal(i).sfreq = hdr.signal(i).nsamples ./ hdr.reclen;
|
|
|
|
end
|
|
|
|
% Find annotations channel
|
|
|
|
iAnnotChans = find(strcmpi({hdr.signal.label}, 'EDF Annotations')); % Mutliple "EDF Annotation" channels allowed in EDF+
|
|
|
|
iStatusChan = find(strcmpi({hdr.signal.label}, 'Status'), 1); % Only one "Status" channel allowed in BDF
|
|
|
|
iOtherChan = setdiff(1:hdr.nsignal, [iAnnotChans iStatusChan]);
|
|
|
|
% % Remove channels with lower sampling rates
|
|
|
|
% iIgnoreChan = find([hdr.signal(iOtherChan).sfreq] < max([hdr.signal(iOtherChan).sfreq])); % Ignore all the channels with lower sampling rate
|
|
|
|
% if ~isempty(iIgnoreChan)
|
|
|
|
% iOtherChan = setdiff(iOtherChan, iIgnoreChan);
|
|
|
|
% end
|
|
|
|
% Get all the other channels
|
|
|
|
if isempty(iOtherChan)
|
|
|
|
error('This file does not contain any data channel.');
|
|
|
|
end
|
|
|
|
% Read events preferencially from the EDF Annotations track
|
|
|
|
if ~isempty(iAnnotChans)
|
|
|
|
iEvtChans = iAnnotChans;
|
|
|
|
elseif ~isempty(iStatusChan)
|
|
|
|
iEvtChans = iStatusChan;
|
|
|
|
% Authors: Francois Tadel, 2012-2017
|
|
|
|
|
|
|
|
%% ===== PARSE INPUTS =====
|
|
|
|
nChannels = sFile.header.nsignal;
|
|
|
|
iChanAnnot = find(strcmpi({sFile.header.signal.label}, 'EDF Annotations'));
|
|
|
|
iBadChan = find(sFile.channelflag == -1);
|
|
|
|
iChanSignal = setdiff(1:nChannels, iChanAnnot);
|
|
|
|
iChanWrongRate = find([sFile.header.signal(iChanSignal).sfreq] ~= max([sFile.header.signal(setdiff(iChanSignal,iBadChan)).sfreq]));
|
|
|
|
iChanSkip = union(iChanAnnot, iChanWrongRate);
|
|
|
|
if (nargin < 4) || isempty(ChannelsRange)
|
|
|
|
ChannelsRange = [1, nChannels];
|
|
|
|
end
|
|
|
|
if (nargin < 3) || isempty(SamplesBounds)
|
|
|
|
SamplesBounds = [0, sFile.header.nrec * sFile.header.signal(ChannelsRange(1)).nsamples - 1];
|
|
|
|
end
|
|
|
|
nTimes = sFile.header.reclen * sFile.header.signal(ChannelsRange(1)).sfreq;
|
|
|
|
iTimes = SamplesBounds(1):SamplesBounds(2);
|
|
|
|
% Block of times/channels to extract
|
|
|
|
nReadChannels = double(ChannelsRange(2) - ChannelsRange(1) + 1);
|
|
|
|
% Read annotations instead of real data ?
|
|
|
|
isAnnotOnly = ~isempty(iChanAnnot) && (ChannelsRange(1) == ChannelsRange(2)) && ismember(ChannelsRange(1), iChanAnnot);
|
|
|
|
isBdfStatus = strcmpi(sFile.format, 'EEG-BDF') && (ChannelsRange(1) == ChannelsRange(2)) && strcmpi(sFile.header.signal(ChannelsRange(1)).label, 'Status');
|
|
|
|
% Data channels to read
|
|
|
|
if isAnnotOnly
|
|
|
|
% Read only on annotation channel
|
|
|
|
iChanF = 1;
|
|
|
|
else
|
|
|
|
iEvtChans = [];
|
|
|
|
end
|
|
|
|
% % Detect channels with inconsistent sampling frenquency
|
|
|
|
% iErrChan = find([hdr.signal(iOtherChan).sfreq] ~= hdr.signal(iOtherChan(1)).sfreq);
|
|
|
|
% iErrChan = setdiff(iErrChan, iAnnotChans);
|
|
|
|
% if ~isempty(iErrChan)
|
|
|
|
% error('Files with mixed sampling rates are not supported yet.');
|
|
|
|
% end
|
|
|
|
% Detect interrupted signals (time non-linear)
|
|
|
|
hdr.interrupted = ischar(hdr.unknown1) && (length(hdr.unknown1) >= 5) && isequal(hdr.unknown1(1:5), 'EDF+D');
|
|
|
|
if hdr.interrupted
|
|
|
|
if ImportOptions.DisplayMessages
|
|
|
|
[res, isCancel] = java_dialog('question', ...
|
|
|
|
['Interrupted EDF file ("EDF+D") detected. It is recommended to convert it' 10 ...
|
|
|
|
'to a continuous ("EDF+C") file first. Do you want to continue reading this' 10 ...
|
|
|
|
'file as continuous and attempt to fix the timing of event markers?' 10 ...
|
|
|
|
'NOTE: This may not work as intended, use at your own risk!']);
|
|
|
|
hdr.fixinterrupted = ~isCancel && strcmpi(res, 'yes');
|
|
|
|
% Remove all the annotation channels from the list of channels to read
|
|
|
|
iChanF = setdiff(ChannelsRange(1):ChannelsRange(2), iChanSkip) - ChannelsRange(1) + 1;
|
|
|
|
% if any(diff(iChanF) ~= 1)
|
|
|
|
% error('All the data channels to read from the file must be contiguous (EDF Annotation channels must be at the end of the list).');
|
|
|
|
% end
|
|
|
|
if any(diff(iChanF) ~= 1)
|
|
|
|
iChanLast = find(diff(iChanF) ~= 1, 1);
|
|
|
|
iChanF = iChanF(1:iChanLast);
|
|
|
|
else
|
|
|
|
hdr.fixinterrupted = 1;
|
|
|
|
iChanLast = length(iChanF);
|
|
|
|
end
|
|
|
|
if ~hdr.fixinterrupted
|
|
|
|
warning(['Interrupted EDF file ("EDF+D"): requires conversion to "EDF+C". ' 10 ...
|
|
|
|
'Brainstorm will read this file as a continuous file ("EDF+C"), the timing of the samples after the first discontinuity will be wrong.' 10 ...
|
|
|
|
'This may not cause any major problem unless there are time markers in the file, they might be inaccurate in all the segments >= 2.']);
|
|
|
|
ChannelsRange = [iChanF(1), iChanF(iChanLast)] + ChannelsRange(1) - 1;
|
|
|
|
end
|
|
|
|
% Cannot read channels with different sampling rates at the same time
|
|
|
|
if (ChannelsRange(1) ~= ChannelsRange(2))
|
|
|
|
allFreq = [sFile.header.signal(ChannelsRange(1):ChannelsRange(2)).nsamples];
|
|
|
|
if any(allFreq ~= allFreq(1))
|
|
|
|
error(['Cannot read channels with mixed sampling rates at the same time.' 10 ...
|
|
|
|
'Mark as bad channels with different sampling rates than EEG.' 10 ...
|
|
|
|
'(right-click on data file > Good/bad channels > Edit good/bad channels)']);
|
|
|
|
end
|
|
|
|
else
|
|
|
|
hdr.fixinterrupted = 0;
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
|
|
|
|
% Initialize returned file structure
|
|
|
|
sFile = db_template('sfile');
|
|
|
|
% Add information read from header
|
|
|
|
sFile.byteorder = 'l';
|
|
|
|
sFile.filename = DataFile;
|
|
|
|
if (uint8(hdr.version(1)) == uint8(255))
|
|
|
|
sFile.format = 'EEG-BDF';
|
|
|
|
sFile.device = 'BDF';
|
|
|
|
%% ===== READ ALL NEEDED EPOCHS =====
|
|
|
|
% Detect which epochs are necessary for the range of data selected
|
|
|
|
epochRange = floor(SamplesBounds ./ nTimes);
|
|
|
|
epochsToRead = epochRange(1) : epochRange(2);
|
|
|
|
% Initialize block of data to read
|
|
|
|
if isAnnotOnly
|
|
|
|
F = zeros(nReadChannels, 2 * length(iTimes));
|
|
|
|
else
|
|
|
|
sFile.format = 'EEG-EDF';
|
|
|
|
sFile.device = 'EDF';
|
|
|
|
F = zeros(nReadChannels, length(iTimes));
|
|
|
|
end
|
|
|
|
% Marker that we increment when we add data to F
|
|
|
|
iF = 1;
|
|
|
|
% Read all the needed epochs
|
|
|
|
for i = 1:length(epochsToRead)
|
|
|
|
% Find the samples to read from this epoch
|
|
|
|
BoundsEpoch = nTimes * epochsToRead(i) + [0, nTimes-1];
|
|
|
|
BoundsRead = [max(BoundsEpoch(1), iTimes(1)), ...
|
|
|
|
min(BoundsEpoch(2), iTimes(end))];
|
|
|
|
iTimeRead = BoundsRead(1):BoundsRead(2);
|
|
|
|
% Convert this samples into indices in this very epoch
|
|
|
|
iTimeRead = iTimeRead - nTimes * epochsToRead(i);
|
|
|
|
% New indices to read
|
|
|
|
if isAnnotOnly
|
|
|
|
iNewF = iF:(iF + 2*length(iTimeRead) - 1);
|
|
|
|
else
|
|
|
|
iNewF = iF:(iF + length(iTimeRead) - 1);
|
|
|
|
end
|
|
|
|
% Read epoch (full or partial)
|
|
|
|
F(iChanF,iNewF) = edf_read_epoch(sFile, sfid, epochsToRead(i), iTimeRead, ChannelsRange, isAnnotOnly, isBdfStatus);
|
|
|
|
% Increment marker
|
|
|
|
iF = iF + length(iTimeRead);
|
|
|
|
end
|
|
|
|
sFile.header = hdr;
|
|
|
|
% Comment: short filename
|
|
|
|
[tmp__, sFile.comment, tmp__] = bst_fileparts(DataFile);
|
|
|
|
% No info on bad channels
|
|
|
|
sFile.channelflag = ones(hdr.nsignal,1);
|
|
|
|
% Acquisition date
|
|
|
|
sFile.acq_date = str_date(hdr.startdate);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
%% ===== PROCESS CHANNEL NAMES/TYPES =====
|
|
|
|
% Try to split the channel names in "TYPE NAME"
|
|
|
|
SplitType = repmat({''}, 1, hdr.nsignal);
|
|
|
|
SplitName = repmat({''}, 1, hdr.nsignal);
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
% Removing trailing dots (eg. "Fc5." instead of "FC5", as in: https://www.physionet.org/pn4/eegmmidb/)
|
|
|
|
if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)
|
|
|
|
hdr.signal(i).label(end) = [];
|
|
|
|
if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)
|
|
|
|
hdr.signal(i).label(end) = [];
|
|
|
|
if (hdr.signal(i).label(end) == '.') && (length(hdr.signal(i).label) > 1)
|
|
|
|
hdr.signal(i).label(end) = [];
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
% Remove extra spaces
|
|
|
|
signalLabel = strrep(hdr.signal(i).label, ' - ', '-');
|
|
|
|
% Find space chars (label format "Type Name")
|
|
|
|
iSpace = find(signalLabel == ' ');
|
|
|
|
% Only if there is one space only
|
|
|
|
if (length(iSpace) == 1) && (iSpace >= 3)
|
|
|
|
SplitName{i} = signalLabel(iSpace+1:end);
|
|
|
|
SplitType{i} = signalLabel(1:iSpace-1);
|
|
|
|
% Accept also 2 spaces
|
|
|
|
elseif (length(iSpace) == 2) && (iSpace(1) >= 3)
|
|
|
|
SplitName{i} = strrep(signalLabel(iSpace(1)+1:end), ' ', '_');
|
|
|
|
SplitType{i} = signalLabel(1:iSpace(1)-1);
|
|
|
|
end
|
|
|
|
end
|
|
|
|
% Remove the classification if it makes some names non unique
|
|
|
|
uniqueNames = unique(SplitName);
|
|
|
|
for i = 1:length(uniqueNames)
|
|
|
|
if ~isempty(uniqueNames{i})
|
|
|
|
iName = find(strcmpi(SplitName, uniqueNames{i}));
|
|
|
|
if (length(iName) > 1)
|
|
|
|
[SplitName{iName}] = deal('');
|
|
|
|
[SplitType{iName}] = deal('');
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
%% ===== CREATE EMPTY CHANNEL FILE =====
|
|
|
|
ChannelMat = db_template('channelmat');
|
|
|
|
ChannelMat.Comment = [sFile.device ' channels'];
|
|
|
|
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, hdr.nsignal]);
|
|
|
|
chRef = {};
|
|
|
|
% For each channel
|
|
|
|
for i = 1:hdr.nsignal
|
|
|
|
% If is the annotation channel
|
|
|
|
if ~isempty(iAnnotChans) && ismember(i, iAnnotChans)
|
|
|
|
ChannelMat.Channel(i).Type = 'EDF';
|
|
|
|
ChannelMat.Channel(i).Name = 'Annotations';
|
|
|
|
elseif ~isempty(iStatusChan) && (i == iStatusChan)
|
|
|
|
ChannelMat.Channel(i).Type = 'BDF';
|
|
|
|
ChannelMat.Channel(i).Name = 'Status';
|
|
|
|
% Regular channels
|
|
|
|
else
|
|
|
|
% If there is a pair name/type already detected
|
|
|
|
if ~isempty(SplitName{i}) && ~isempty(SplitType{i})
|
|
|
|
ChannelMat.Channel(i).Name = SplitName{i};
|
|
|
|
ChannelMat.Channel(i).Type = SplitType{i};
|
|
|
|
|
|
|
|
%% ===== READ ONE EPOCH =====
|
|
|
|
function F = edf_read_epoch(sFile, sfid, iEpoch, iTimes, ChannelsRange, isAnnotOnly, isBdfStatus)
|
|
|
|
% ===== COMPUTE OFFSETS =====
|
|
|
|
nTimes = sFile.header.reclen * [sFile.header.signal.sfreq];
|
|
|
|
nReadTimes = length(iTimes);
|
|
|
|
nReadChannels = double(ChannelsRange(2) - ChannelsRange(1) + 1);
|
|
|
|
iChannels = ChannelsRange(1):ChannelsRange(2);
|
|
|
|
% Check that all the channels selected have the same freq rate
|
|
|
|
if any(nTimes(iChannels) ~= nTimes(iChannels(1)))
|
|
|
|
error('Cannot read at the same signals with different sampling frequency.');
|
|
|
|
end
|
|
|
|
% Size of one value
|
|
|
|
if strcmpi(sFile.format, 'EEG-BDF')
|
|
|
|
% BDF: int24 => 3 bytes
|
|
|
|
bytesPerVal = 3;
|
|
|
|
% Reading status or regular channel
|
|
|
|
if isBdfStatus
|
|
|
|
dataClass = 'ubit24';
|
|
|
|
else
|
|
|
|
% Channel name
|
|
|
|
ChannelMat.Channel(i).Name = hdr.signal(i).label(hdr.signal(i).label ~= ' ');
|
|
|
|
% Channel type
|
|
|
|
if ~isempty(hdr.signal(i).type)
|
|
|
|
if (length(hdr.signal(i).type) == 3)
|
|
|
|
ChannelMat.Channel(i).Type = hdr.signal(i).type(hdr.signal(i).type ~= ' ');
|
|
|
|
elseif isequal(hdr.signal(i).type, 'Active Electrode') || isequal(hdr.signal(i).type, 'AgAgCl electrode')
|
|
|
|
ChannelMat.Channel(i).Type = 'EEG';
|
|
|
|
else
|
|
|
|
ChannelMat.Channel(i).Type = 'Misc';
|
|
|
|
end
|
|
|
|
else
|
|
|
|
ChannelMat.Channel(i).Type = 'EEG';
|
|
|
|
end
|
|
|
|
end
|
|
|
|
% Extract reference name (at the end of the channel name, separated with a "-", eg. "-REF")
|
|
|
|
iDash = find(ChannelMat.Channel(i).Name == '-');
|
|
|
|
if ~isempty(iDash) && (iDash(end) < length(ChannelMat.Channel(i).Name))
|
|
|
|
chRef{end+1} = ChannelMat.Channel(i).Name(iDash(end):end);
|
|
|
|
dataClass = 'bit24';
|
|
|
|
end
|
|
|
|
else
|
|
|
|
% EDF: int16 => 2 bytes
|
|
|
|
bytesPerVal = 2;
|
|
|
|
dataClass = 'int16';
|
|
|
|
isBdfStatus = 0;
|
|
|
|
end
|
|
|
|
ChannelMat.Channel(i).Loc = [0; 0; 0];
|
|
|
|
ChannelMat.Channel(i).Orient = [];
|
|
|
|
ChannelMat.Channel(i).Weight = 1;
|
|
|
|
% ChannelMat.Channel(i).Comment = hdr.signal(i).type;
|
|
|
|
end
|
|
|
|
|
|
|
|
% If the same reference is indicated for all the channels: remove it
|
|
|
|
if (length(chRef) >= 2)
|
|
|
|
% Get the shortest reference tag
|
|
|
|
lenRef = cellfun(@length, chRef);
|
|
|
|
minLen = min(lenRef);
|
|
|
|
% Check if all the ref names are equal (up to the max length - some might be cut because the channel name is too long)
|
|
|
|
if all(cellfun(@(c)strcmpi(c(1:minLen), chRef{1}(1:minLen)), chRef))
|
|
|
|
% Remove the reference tag from all the channel names
|
|
|
|
for i = 1:length(ChannelMat.Channel)
|
|
|
|
ChannelMat.Channel(i).Name = strrep(ChannelMat.Channel(i).Name, chRef{1}, '');
|
|
|
|
ChannelMat.Channel(i).Name = strrep(ChannelMat.Channel(i).Name, chRef{1}(1:minLen), '');
|
|
|
|
% Offset of the beginning of the recordings in the file
|
|
|
|
offsetHeader = round(sFile.header.hdrlen);
|
|
|
|
% Offset of epoch
|
|
|
|
offsetEpoch = round(iEpoch * sum(nTimes) * bytesPerVal);
|
|
|
|
% Channel offset
|
|
|
|
offsetChannel = round(sum(nTimes(1:ChannelsRange(1)-1)) * bytesPerVal);
|
|
|
|
% Time offset at the beginning and end of each channel block
|
|
|
|
offsetTimeStart = round(iTimes(1) * bytesPerVal);
|
|
|
|
offsetTimeEnd = round((nTimes(ChannelsRange(1)) - iTimes(end) - 1) * bytesPerVal);
|
|
|
|
% ALL THE "ROUND" CALLS WERE ADDED AFTER DISCOVERING THAT THERE WERE SOMETIMES ROUNDING ERRORS IN THE MULTIPLICATIONS
|
|
|
|
|
|
|
|
% Where to start reading in the file ?
|
|
|
|
% => After the header, the number of skipped epochs, channels and time samples
|
|
|
|
offsetStart = offsetHeader + offsetEpoch + offsetChannel + offsetTimeStart;
|
|
|
|
% Number of time samples to skip after each channel
|
|
|
|
offsetSkip = offsetTimeStart + offsetTimeEnd;
|
|
|
|
|
|
|
|
|
|
|
|
% ===== READ DATA BLOCK =====
|
|
|
|
% Position file at the beginning of the trial
|
|
|
|
fseek(sfid, offsetStart, 'bof');
|
|
|
|
% Read annotation data (char)
|
|
|
|
if isAnnotOnly
|
|
|
|
dataClass = 'char';
|
|
|
|
nReadTimes = bytesPerVal * nReadTimes; % 1 byte instead of 2
|
|
|
|
end
|
|
|
|
% Read trial data
|
|
|
|
% => WARNING: CALL TO FREAD WITH SKIP=0 DOES NOT WORK PROPERLY
|
|
|
|
if (offsetSkip == 0)
|
|
|
|
F = fread(sfid, [nReadTimes, nReadChannels], dataClass)';
|
|
|
|
elseif (bytesPerVal == 2)
|
|
|
|
precision = sprintf('%d*%s', nReadTimes, dataClass);
|
|
|
|
F = fread(sfid, [nReadTimes, nReadChannels], precision, offsetSkip)';
|
|
|
|
% => WARNING: READING USING ubit24 SOMETIMES DOESNT WORK => DOING IT MANUALLY
|
|
|
|
elseif (bytesPerVal == 3)
|
|
|
|
% Reading each bit independently
|
|
|
|
precision = sprintf('%d*%s', 3*nReadTimes, 'uint8');
|
|
|
|
F = fread(sfid, [3*nReadTimes, nReadChannels], precision, offsetSkip)';
|
|
|
|
% Grouping the 3 bits together
|
|
|
|
F = F(:,1:3:end) + F(:,2:3:end)*256 + F(:,3:3:end)*256*256;
|
|
|
|
% 2-Complement (negative value indicated by most significant bit)
|
|
|
|
if strcmpi(dataClass, 'bit24')
|
|
|
|
iNeg = (F >= 256*256*128);
|
|
|
|
F(iNeg) = F(iNeg) - 256*256*256;
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
% If there are only "Misc" and no "EEG" channels: rename to "EEG"
|
|
|
|
iMisc = find(strcmpi({ChannelMat.Channel.Type}, 'Misc'));
|
|
|
|
iEeg = find(strcmpi({ChannelMat.Channel.Type}, 'EEG'));
|
|
|
|
if ~isempty(iMisc) && isempty(iEeg)
|
|
|
|
[ChannelMat.Channel(iMisc).Type] = deal('EEG');
|
|
|
|
iEeg = iMisc;
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
%% ===== DETECT MULTIPLE SAMPLING RATES =====
|
|
|
|
% Use the first "EEG" channel as the reference sampling rate (or the first channel if no "EEG" channels available)
|
|
|
|
if ~isempty(iEeg) && ismember(iEeg(1), iOtherChan)
|
|
|
|
iChanFreqRef = iEeg(1);
|
|
|
|
else
|
|
|
|
iChanFreqRef = iOtherChan(1);
|
|
|
|
end
|
|
|
|
% Mark as bad channels with sampling rates different from EEG
|
|
|
|
iChanWrongRate = find([sFile.header.signal.sfreq] ~= sFile.header.signal(iChanFreqRef).sfreq);
|
|
|
|
iChanWrongRate = intersect(iChanWrongRate, iOtherChan);
|
|
|
|
if ~isempty(iChanWrongRate)
|
|
|
|
sFile.channelflag(iChanWrongRate) = -1;
|
|
|
|
end
|
|
|
|
|
|
|
|
% Consider that the sampling rate of the file is the sampling rate of the first signal
|
|
|
|
sFile.prop.sfreq = hdr.signal(iChanFreqRef).sfreq;
|
|
|
|
sFile.prop.times = [0, hdr.signal(iChanFreqRef).nsamples * hdr.nrec - 1] ./ sFile.prop.sfreq;
|
|
|
|
sFile.prop.nAvg = 1;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
% Check that data block was fully read
|
|
|
|
if (numel(F) < nReadTimes * nReadChannels)
|
|
|
|
% Number of full time samples that were read
|
|
|
|
nTimeTrunc = max(0, floor(numel(F) / nReadChannels) - 1);
|
|
|
|
% Error message
|
|
|
|
disp(sprintf('EDF> ERROR: File is truncated (%d time samples were read instead of %d). Padding with zeros...', nTimeTrunc, nReadTimes));
|
|
|
|
% Pad data with zeros
|
|
|
|
Ftmp = zeros(nReadTimes, nReadChannels);
|
|
|
|
F = F';
|
|
|
|
Ftmp(1:numel(F)) = F(:);
|
|
|
|
F = Ftmp';
|
|
|
|
end
|
|
|
|
|
|
|
|
%% ===== READ EDF ANNOTATION CHANNEL =====
|
|
|
|
if ~isempty(iEvtChans) % && ~isequal(ImportOptions.EventsMode, 'ignore')
|
|
|
|
% Set reading options
|
|
|
|
ImportOptions.ImportMode = 'Time';
|
|
|
|
ImportOptions.UseSsp = 0;
|
|
|
|
ImportOptions.UseCtfComp = 0;
|
|
|
|
% Read EDF annotations
|
|
|
|
if strcmpi(sFile.format, 'EEG-EDF')
|
|
|
|
evtList = {};
|
|
|
|
% In EDF+, the first annotation channel has epoch time stamps (EDF
|
|
|
|
% calls epochs records). So read all annotation channels per epoch.
|
|
|
|
for irec = 1:hdr.nrec
|
|
|
|
for ichan = 1:length(iEvtChans)
|
|
|
|
bst_progress('text', sprintf('Reading annotations... [%d%%]', round((ichan + (irec-1)*length(iEvtChans))/length(iEvtChans)/hdr.nrec*100)));
|
|
|
|
% Sample indices for the current epoch (=record)
|
|
|
|
SampleBounds = [irec-1,irec] * sFile.header.signal(iEvtChans(ichan)).nsamples - [0,1];
|
|
|
|
% Read record
|
|
|
|
F = in_fread(sFile, ChannelMat, 1, SampleBounds, iEvtChans(ichan), ImportOptions);
|
|
|
|
% Find all the TALs (Time-stamped Annotations Lists): separated with [20][0]
|
|
|
|
% Possible configurations:
|
|
|
|
% Onset[21]Duration[20]Annot1[20]Annot2[20]..AnnotN[20][0]
|
|
|
|
% Onset[20]Annot1[20]Annot2[20]..AnnotN[20][0]
|
|
|
|
% Onset[20]Annot1[20][0]
|
|
|
|
iSeparator = [-1, find((F(1:end-1) == 20) & (F(2:end) == 0))];
|
|
|
|
% Loop on all the TALs
|
|
|
|
for iAnnot = 1:length(iSeparator)-1
|
|
|
|
% Get annotation
|
|
|
|
strAnnot = char(F(iSeparator(iAnnot)+2:iSeparator(iAnnot+1)-1));
|
|
|
|
% Split in blocks with [20]
|
|
|
|
splitAnnot = str_split(strAnnot, 20);
|
|
|
|
% The first TAL in a record should be indicating the timing of the first sample of the record
|
|
|
|
if (iAnnot == 1) && (length(splitAnnot) == 1) && ~any(splitAnnot{1} == 21)
|
|
|
|
% Ignore if this is not the first channel
|
|
|
|
if (ichan > 1)
|
|
|
|
continue;
|
|
|
|
end
|
|
|
|
% Get record time stamp
|
|
|
|
t0_rec = str2double(splitAnnot{1});
|
|
|
|
if isempty(t0_rec) || isnan(t0_rec)
|
|
|
|
continue;
|
|
|
|
end
|
|
|
|
if (irec == 1)
|
|
|
|
t0_file = t0_rec;
|
|
|
|
% Find discontinuities larger than 1 sample
|
|
|
|
elseif abs(t0_rec - prev_rec - (irec - prev_irec) * hdr.reclen) > (1 / sFile.prop.sfreq)
|
|
|
|
% Brainstorm fills partial/interrupted records with zeros
|
|
|
|
bstTime = prev_rec + hdr.reclen;
|
|
|
|
timeDiff = bstTime - t0_rec;
|
|
|
|
% If we want to fix timing, apply skip to initial timestamp
|
|
|
|
if hdr.fixinterrupted
|
|
|
|
t0_file = t0_file - timeDiff;
|
|
|
|
end
|
|
|
|
% Warn user of discontinuity
|
|
|
|
if timeDiff > 0
|
|
|
|
expectMsg = 'blank data';
|
|
|
|
else
|
|
|
|
expectMsg = 'skipped data';
|
|
|
|
end
|
|
|
|
startTime = min(t0_rec - t0_file - [0, timeDiff]); % before and after t0_file adjustment
|
|
|
|
endTime = max(t0_rec - t0_file - [0, timeDiff]);
|
|
|
|
fprintf('WARNING: Found discontinuity between %.3fs and %.3fs, expect %s in between.\n', startTime, endTime, expectMsg);
|
|
|
|
% Create event for users information
|
|
|
|
if timeDiff < 0
|
|
|
|
endTime = startTime; % no extent in this case, there is skipped time.
|
|
|
|
end
|
|
|
|
evtList(end+1,:) = {'EDF+D Discontinuity', [startTime; endTime]};
|
|
|
|
end
|
|
|
|
prev_rec = t0_rec;
|
|
|
|
prev_irec = irec;
|
|
|
|
% Regular TAL: indicating a marker
|
|
|
|
else
|
|
|
|
% Split time in onset/duration
|
|
|
|
splitTime = str_split(splitAnnot{1}, 21);
|
|
|
|
% Get time
|
|
|
|
t = str2double(splitTime{1});
|
|
|
|
if isempty(t) || isnan(t)
|
|
|
|
continue;
|
|
|
|
end
|
|
|
|
% Get duration
|
|
|
|
if (length(splitTime) > 1)
|
|
|
|
duration = str2double(splitTime{2});
|
|
|
|
% Exclude 1-sample long events
|
|
|
|
if isempty(duration) || isnan(duration) || (round(duration .* sFile.prop.sfreq) <= 1)
|
|
|
|
duration = 0;
|
|
|
|
end
|
|
|
|
else
|
|
|
|
duration = 0;
|
|
|
|
end
|
|
|
|
% Unnamed event
|
|
|
|
if (length(splitAnnot) == 1) || isempty(splitAnnot{2})
|
|
|
|
splitAnnot{2} = 'Unnamed';
|
|
|
|
end
|
|
|
|
% Create one event for each label in the TAL
|
|
|
|
for iLabel = 2:length(splitAnnot)
|
|
|
|
evtList(end+1,:) = {splitAnnot{iLabel}, (t-t0_file) + [0;duration]};
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
% If there are events: create a create an events structure
|
|
|
|
if ~isempty(evtList)
|
|
|
|
% Initialize events list
|
|
|
|
sFile.events = repmat(db_template('event'), 0);
|
|
|
|
% Events list
|
|
|
|
[uniqueEvt, iUnique] = unique(evtList(:,1));
|
|
|
|
uniqueEvt = evtList(sort(iUnique),1);
|
|
|
|
% Build events list
|
|
|
|
for iEvt = 1:length(uniqueEvt)
|
|
|
|
% Find all the occurrences of this event
|
|
|
|
iOcc = find(strcmpi(uniqueEvt{iEvt}, evtList(:,1)));
|
|
|
|
% Concatenate all times
|
|
|
|
t = [evtList{iOcc,2}];
|
|
|
|
% If second row is equal to the first one (no extended events): delete it
|
|
|
|
if all(t(1,:) == t(2,:))
|
|
|
|
t = t(1,:);
|
|
|
|
end
|
|
|
|
% Set event
|
|
|
|
sFile.events(iEvt).label = strtrim(uniqueEvt{iEvt});
|
|
|
|
sFile.events(iEvt).times = round(t .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
|
|
|
|
sFile.events(iEvt).epochs = 1 + 0*t(1,:);
|
|
|
|
sFile.events(iEvt).select = 1;
|
|
|
|
sFile.events(iEvt).channels = cell(1, size(sFile.events(iEvt).times, 2));
|
|
|
|
sFile.events(iEvt).notes = cell(1, size(sFile.events(iEvt).times, 2));
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
% BDF Status line
|
|
|
|
elseif strcmpi(sFile.format, 'EEG-BDF') && ~strcmpi(ImportOptions.EventsTrackMode, 'ignore')
|
|
|
|
% Ask how to read the events
|
|
|
|
[events, ImportOptions.EventsTrackMode] = process_evt_read('Compute', sFile, ChannelMat, ChannelMat.Channel(iEvtChans).Name, ImportOptions.EventsTrackMode);
|
|
|
|
if isequal(events, -1)
|
|
|
|
sFile = [];
|
|
|
|
ChannelMat = [];
|
|
|
|
return;
|
|
|
|
end
|
|
|
|
% Report the events in the file structure
|
|
|
|
sFile.events = events;
|
|
|
|
% Remove the 'Status: ' string in front of the events
|
|
|
|
for i = 1:length(sFile.events)
|
|
|
|
sFile.events(i).label = strrep(sFile.events(i).label, 'Status: ', '');
|
|
|
|
end
|
|
|
|
% Group events by time
|
|
|
|
% sFile.events = process_evt_grouptime('Compute', sFile.events);
|
|
|
|
% Processing for BDF status file
|
|
|
|
if isBdfStatus
|
|
|
|
% Mask to keep only the first 15 bits (Triggers bits)
|
|
|
|
% Bit 16 : High when new Epoch is started
|
|
|
|
% Bit 17-19 : Speed bits 0 1 2
|
|
|
|
% Bit 20 : High when CMS is within range
|
|
|
|
% Bit 21 : Speed bit 3
|
|
|
|
% Bit 22 : High when battery is low
|
|
|
|
% Bit 23 : High if ActiveTwo MK2
|
|
|
|
F = bitand(F, bin2dec('000000000111111111111111'));
|
|
|
|
% Processing for real data
|
|
|
|
elseif ~isAnnotOnly
|
|
|
|
% Convert to double
|
|
|
|
F = double(F);
|
|
|
|
% Apply gains
|
|
|
|
F = bst_bsxfun(@rdivide, F, [sFile.header.signal(iChannels).gain]');
|
|
|
|
% IN THEORY: THIS OFFSET SECTION IS USEFUL, BUT IN PRACTICE IT LOOKS LIKE THE VALUES IN ALL THE FILES ARE CENTERED ON ZERO
|
|
|
|
% % Add offset
|
|
|
|
% if isfield(sFile.header.signal, 'offset') && ~isempty(sFile.header.signal(1).offset)
|
|
|
|
% % ...
|
|
|
|
% end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
...
|
...
|
|
