/公众号：好玩的Matlab
小编是个从事脑电、肌电、心电方面的，在处理数据时候遇到edf文件格式的处理，经过查询资料终于找到了数据处理方法。
EEGLAB工具箱：EEGLAB
目录
Matlab读取edf文件方法(博瑞康数据放大器)
Matlab读取edf文件方法
function [hdr, record] = edfread(fname, varargin)
% Read European Data Format file into MATLAB
%
% [hdr, record] = edfread(fname)
% Reads data from ALL RECORDS of file fname ('*.edf'). Header
% information is returned in structure hdr, and the signals
% (waveforms) are returned in structure record, with waveforms
% associated with the records returned as fields titled 'data' of
% structure record.
%
% [...] = edfread(fname, 'assignToVariables', assignToVariables)
% If assignToVariables is true, triggers writing of individual
% output variables, as defined by field 'labels', into the caller
% workspace.
%
% [...] = edfread(...,'targetSignals',targetSignals)
% Allows user to specify the names (or position numbers) of the
% subset of signals to be read. |targetSignals| may be either a
% string, a cell array of comma-separated strings, or a vector of
% numbers. (Default behavior is to read all signals.)
% E.g.:
% data = edfread(mydata.edf,'targetSignals','Thoracic');
% data = edfread(mydata.edf,'targetSignals',{'Thoracic1','Abdominal'});
% or
% data = edfread(mydata.edf,'targetSignals',[2,4,6:13]);
%
% FORMAT SPEC: Source: http://www.edfplus.info/specs/edf.html SEE ALSO:
% http://www.dpmi.tu-graz.ac.at/~schloegl/matlab/eeg/edf_spec.htm
%
% The first 256 bytes of the header record specify the version number of
% this format, local patient and recording identification, time information
% about the recording, the number of data records and finally the number of
% signals (ns) in each data record. Then for each signal another 256 bytes
% follow in the header record, each specifying the type of signal (e.g.
% EEG, body temperature, etc.), amplitude calibration and the number of
% samples in each data record (from which the sampling frequency can be
% derived since the duration of a data record is also known). In this way,
% the format allows for different gains and sampling frequencies for each
% signal. The header record contains 256 + (ns * 256) bytes.
%
% Following the header record, each of the subsequent data records contains
% 'duration' seconds of 'ns' signals, with each signal being represented by
% the specified (in the header) number of samples. In order to reduce data
% size and adapt to commonly used software for acquisition, processing and
% graphical display of polygraphic signals, each sample value is
% represented as a 2-byte integer in 2's complement format. Figure 1 shows
% the detailed format of each data record.
%
% DATA SOURCE: Signals of various types (including the sample signal used
% below) are available from PHYSIONET: http://www.physionet.org/
%
%
% % EXAMPLE 1:
% % Read all waveforms/data associated with file 'ecgca998.edf':
%
% [header, recorddata] = edfread('ecgca998.edf');
%
% % EXAMPLE 2:
% % Read records 3 and 5, associated with file 'ecgca998.edf':
%
% % E
% header = edfread('ecgca998.edf','AssignToVariables',true);
% % Header file specifies data labels 'label_1'...'label_n'; these are
% % created as variables in the caller workspace.
%
% Coded 8/27/09 by Brett Shoelson, PhD
% brett.shoelson@mathworks.com
% Copyright 2009 - 2012 MathWorks, Inc.
%
% Modifications:
% 5/6/13 Fixed a problem with a poorly subscripted variable. (Under certain
% conditions, data were being improperly written to the 'records' variable.
% Thanks to Hisham El Moaqet for reporting the problem and for sharing a
% file that helped me track it down.)
%
% 5/22/13 Enabled import of a user-selected subset of signals. Thanks to
% Farid and Cindy for pointing out the deficiency. Also fixed the import of
% signals that had "bad" characters (spaces, etc) in their names.
%
% 10/30/14 Now outputs frequency field directly, and (thanks to Olga Imas
% for the suggestion.)
%
% 3/23/2015 Fixed a doc problem with an improperly named variable
% (desiredSignals->targetSignals).
% HEADER RECORD
% 8 ascii : version of this data format (0)
% 80 ascii : local patient identification
% 80 ascii : local recording identification
% 8 ascii : startdate of recording (dd.mm.yy)
% 8 ascii : starttime of recording (hh.mm.ss)
% 8 ascii : number of bytes in header record
% 44 ascii : reserved
% 8 ascii : number of data records (-1 if unknown)
% 8 ascii : duration of a data record, in seconds
% 4 ascii : number of signals (ns) in data record
% ns * 16 ascii : ns * label (e.g. EEG FpzCz or Body temp)
% ns * 80 ascii : ns * transducer type (e.g. AgAgCl electrode)
% ns * 8 ascii : ns * physical dimension (e.g. uV or degreeC)
% ns * 8 ascii : ns * physical minimum (e.g. -500 or 34)
% ns * 8 ascii : ns * physical maximum (e.g. 500 or 40)
% ns * 8 ascii : ns * digital minimum (e.g. -2048)
% ns * 8 ascii : ns * digital maximum (e.g. 2047)
% ns * 80 ascii : ns * prefiltering (e.g. HP:0.1Hz LP:75Hz)
% ns * 8 ascii : ns * nr of samples in each data record
% ns * 32 ascii : ns * reserved
% DATA RECORD
% nr of samples[1] * integer : first signal in the data record
% nr of samples[2] * integer : second signal
% ..
% ..
% nr of samples[ns] * integer : last signal
if nargin > 5
error('EDFREAD: Too many input arguments.');
end
if ~nargin
error('EDFREAD: Requires at least one input argument (filename to read).');
end
[fid,msg] = fopen(fname,'r');
if fid == -1
error(msg)
end
assignToVariables = false; %Default
targetSignals = []; %Default
for ii = 1:2:numel(varargin)
switch lower(varargin{ii})
case 'assigntovariables'
assignToVariables = varargin{ii+1};
case 'targetsignals'
targetSignals = varargin{ii+1};
otherwise
error('EDFREAD: Unrecognized parameter-value pair specified. Valid values are ''assignToVariables'' and ''targetSignals''.')
end
end
% HEADER
hdr.ver = str2double(char(fread(fid,8)'));
hdr.patientID = fread(fid,80,'*char')';
hdr.recordID = fread(fid,80,'*char')';
hdr.startdate = fread(fid,8,'*char')';% (dd.mm.yy)
% hdr.startdate = datestr(datenum(fread(fid,8,'*char')','dd.mm.yy'), 29); %'yyyy-mm-dd' (ISO 8601)
hdr.starttime = fread(fid,8,'*char')';% (hh.mm.ss)
% hdr.starttime = datestr(datenum(fread(fid,8,'*char')','hh.mm.ss'), 13); %'HH:MM:SS' (ISO 8601)
hdr.bytes = str2double(fread(fid,8,'*char')');
reserved = fread(fid,44);%#ok
hdr.records = str2double(fread(fid,8,'*char')');
if hdr.records == -1
beep
disp('There appears to be a problem with this file; it returns an out-of-spec value of -1 for ''numberOfRecords.''')
disp('Attempting to read the file with ''edfReadUntilDone'' instead....');
[hdr, record] = edfreadUntilDone(fname, varargin);
return
end
hdr.duration = str2double(fread(fid,8,'*char')');
% Number of signals
hdr.ns = str2double(fread(fid,4,'*char')');
for ii = 1:hdr.ns
hdr.label{ii} = regexprep(fread(fid,16,'*char')','\W','');
end
if isempty(targetSignals)
targetSignals = 1:numel(hdr.label);
elseif iscell(targetSignals)||ischar(targetSignals)
targetSignals = find(ismember(hdr.label,regexprep(targetSignals,'\W','')));
end
if isempty(targetSignals)
error('EDFREAD: The signal(s) you requested were not detected.')
end
for ii = 1:hdr.ns
hdr.transducer{ii} = fread(fid,80,'*char')';
end
% Physical dimension
for ii = 1:hdr.ns
hdr.units{ii} = fread(fid,8,'*char')';
end
% Physical minimum
for ii = 1:hdr.ns
hdr.physicalMin(ii) = str2double(fread(fid,8,'*char')');
end
% Physical maximum
for ii = 1:hdr.ns
hdr.physicalMax(ii) = str2double(fread(fid,8,'*char')');
end
% Digital minimum
for ii = 1:hdr.ns
hdr.digitalMin(ii) = str2double(fread(fid,8,'*char')');
end
% Digital maximum
for ii = 1:hdr.ns
hdr.digitalMax(ii) = str2double(fread(fid,8,'*char')');
end
for ii = 1:hdr.ns
hdr.prefilter{ii} = fread(fid,80,'*char')';
end
for ii = 1:hdr.ns
hdr.samples(ii) = str2double(fread(fid,8,'*char')');
end
for ii = 1:hdr.ns
reserved = fread(fid,32,'*char')';%#ok
end
hdr.label = hdr.label(targetSignals);
hdr.label = regexprep(hdr.label,'\W','');
hdr.units = regexprep(hdr.units,'\W','');
hdr.frequency = hdr.samples./hdr.duration;
disp('Step 1 of 2: Reading requested records. (This may take a few minutes.)...');
if nargout > 1 || assignToVariables
% Scale data (linear scaling)
scalefac = (hdr.physicalMax - hdr.physicalMin)./(hdr.digitalMax - hdr.digitalMin);
dc = hdr.physicalMax - scalefac .* hdr.digitalMax;
% RECORD DATA REQUESTED
tmpdata = struct;
for recnum = 1:hdr.records
for ii = 1:hdr.ns
% Read or skip the appropriate number of data points
if ismember(ii,targetSignals)
% Use a cell array for DATA because number of samples may vary
% from sample to sample
tmpdata(recnum).data{ii} = fread(fid,hdr.samples(ii),'int16') * scalefac(ii) + dc(ii);
else
fseek(fid,hdr.samples(ii)*2,0);
end
end
end
hdr.units = hdr.units(targetSignals);
hdr.physicalMin = hdr.physicalMin(targetSignals);
hdr.physicalMax = hdr.physicalMax(targetSignals);
hdr.digitalMin = hdr.digitalMin(targetSignals);
hdr.digitalMax = hdr.digitalMax(targetSignals);
hdr.prefilter = hdr.prefilter(targetSignals);
hdr.transducer = hdr.transducer(targetSignals);
record = zeros(numel(hdr.label), hdr.samples(1)*hdr.records);
% NOTE: 5/6/13 Modified for loop below to change instances of hdr.samples to
% hdr.samples(ii). I think this underscored a problem with the reader.
% Sort by requested targetSignals order:
disp('Step 2 of 2: Parsing data...');
recnum = 1;
for ii = 1:hdr.ns
if ismember(ii,targetSignals)
ctr = 1;
for jj = 1:hdr.records
try %#ok
record(recnum, ctr : ctr + hdr.samples(ii) - 1) = tmpdata(jj).data{ii};
end
ctr = ctr + hdr.samples(ii);
end
recnum = recnum + 1;
end
end
hdr.ns = numel(hdr.label);
hdr.samples = hdr.samples(targetSignals);
if assignToVariables
for ii = 1:numel(hdr.label)
try %#ok
eval(['assignin(''caller'',''',hdr.label{ii},''',record(ii,:))'])
end
end
% Uncomment this line to duplicate output in a single matrix
% ''record''. (Could be memory intensive if dataset is large.)
record = [];% No need to output record data as a matrix?
end
end
fclose(fid);
Matlab读取bdf数据格式的方法
function EEG = readbdfdata(filename, pathname)
%
% Syntax: EEG = readbdfdata(pathname)
%
%
% Inputs:
% filename:
% pathname: path to data files
% Outputs:
% EEG data structure
%
% Example:
%
%
%
% Subfunctions: read_bdf, readLowLevel, readEvents
% MAT-files required: none
%
% Author: Xiaoshan Huang, hxs@neuracle.cn
%
% Versions:
% v0.1: 2017-09-27, orignal
% v0.2: 2018-08-07, fix update events bug by FANG Junying
% Copyright (c) 2017 Neuracle, Inc. All Rights Reserved. http://neuracle.cn
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find the index of data file
index = [];
if ~iscell(filename)
filename = {filename};
end
for i = 1:length(filename)
if length(filename{i}) >= 4 %length('data')==4
if isequal(filename{i}(1:4), 'data')
index = [index, i];
end
end
end
datafilelength = length(index);
%read data files
if datafilelength >0
%add information for eeglab structure
EEG.setname = 'BDF file';
EEG.filename = char(filename);
EEG.filepath = pathname;
EEG.subject = '';
EEG.group = '';
EEG.condition = '';
EEG.session = [];
EEG.comments = ['Original files:' pathname];
EEG.icaact = [];
EEG.icawinv = [];
EEG.icasphere = [];
EEG.icaweights = [];
EEG.icachansind = [];
chaninfo.plotrad = [];
chaninfo.shrink = [];
chaninfo.nosedir = '+X';
chaninfo.icachansind = [];
EEG.chaninfo = chaninfo;
EEG.ref = 'common';
EEG.specdata = [];
EEG.specicaact = [];
EEG.splinefile = '';
EEG.icasplinefile = '';
EEG.dipfit = [];
EEG.xmin = 0;
EEG.xmax = 0;
EEG.pnts = 0;
T0 = [];
EEG.urevent = [];
EEG.eventdescription = {};
EEG.epoch = [];
EEG.epochdescription = {};
EEG.data = [];
Data = {};
datapnts = [];
for i = 1:datafilelength
datafilename = [pathname filename{index(i)}];
%read header
hdr = read_bdf(datafilename);
EEG.srate = hdr.Fs;
EEG.nbchan = hdr.nChans;
datapnts = [datapnts,hdr.nSamples];
EEG.pnts = EEG.pnts+hdr.nSamples;
EEG.chanlocs= hdr.chanlocs;
EEG.trials = hdr.nTrials;
EEG.xmax = EEG.xmax+EEG.pnts/EEG.srate; %in seconds
T0 = [T0; hdr.T0];
%compatible with older version data file which has event channel
if hdr.Annotation == 1
evt = hdr.event;
evt = cell2mat(evt);
EEG.event = struct([]);
if isstruct(evt)
for kk = 1:size(evt,2)
event = struct('type',evt(kk).eventvalue,'latency',round(evt(kk).offset_in_sec*EEG.srate));
EEG.event = [EEG.event;event];
end
end
end
if datafilelength == 1
EEG.data = read_bdf(datafilename,hdr,1,hdr.nSamples);
else
Data{i} = read_bdf(datafilename,hdr,1,hdr.nSamples);
end
end
EEG.times = [1:EEG.pnts]*1000/EEG.srate;
for j=1:size(T0,1)
startsecs(j) = T0(j,1)*946080000+T0(j,2)*2592000+T0(j,3)*86400+T0(j,4)*3600+T0(j,5)*60+T0(j,6);
end
[startsecs,indexsorted] = sort(startsecs);
T0 = T0(indexsorted,:);
if datafilelength > 1
for k = 1:datafilelength
EEG.data = [EEG.data,Data{indexsorted(k)}];
end
end
etc.T0 = T0(1,:);
EEG.etc = etc;
%read event
ind = find(ismember(filename,'evt.bdf'));
if numel(ind)
eventfilename = [pathname filename{ind}];
hdr = read_bdf(eventfilename);
evt = hdr.event;
evt = cell2mat(evt);
EEG.event = struct([]);
if isstruct(evt)
for i = 1:size(evt,2)
event = struct('type',evt(i).eventvalue,'latency',round(evt(i).offset_in_sec*EEG.srate));
EEG.event = [EEG.event;event];
end
end
%shift event latency for multiple data files only
if datafilelength > 1
startpnts = (startsecs-startsecs(1))*EEG.srate; %relative time
endpnts = zeros(1,datafilelength);
pausepnts = 0;
Event = struct([]);
datapnts = datapnts(indexsorted);
for j = 1:datafilelength
endpnts(j) = startpnts(j)+datapnts(j);
if j > 1
pausepnts = pausepnts+(startpnts(j)-endpnts(j-1));
end
for k = 1:length(EEG.event)
if startpnts(j) <= EEG.event(k).latency && endpnts(j) >= EEG.event(k).latency
EEG.event(k).latency = EEG.event(k).latency-pausepnts;
Event = [Event;EEG.event(k)];
end
end
if j < datafilelength
eventBoundary = struct('type','boundary','latency',endpnts(j)+1);
Event = [Event;eventBoundary];
end
end
EEG.event = Event;
end
elseif exist([pathname 'evt.bdf'],'file') == 0 && hdr.Annotation == 0
EEG.event = struct([]);
end
%read mems data and combine with EEG data
memstype = {'acc.edf','gyro.edf','mag.edf'};
for gg = 1:length(memstype)
%initialize
indexmems = [];
T0mems = [];
memsData = {};
mems.data = [];
%find the index of the mems file
for ii = 1:length(filename)
lennametype = length(memstype{gg})-4; %ignore '.edf'
if length(filename{ii}) >= lennametype
if isequal(filename{ii}(1:lennametype), memstype{gg}(1:lennametype))
indexmems = [indexmems, ii];
end
end
end
filelength = length(indexmems);
if filelength >0
for ii = 1:filelength
memsfilename = [pathname filename{indexmems(ii)}];
%read header
hdrmems = read_bdf(memsfilename);
chanlocs= hdrmems.chanlocs;
nbchan = hdrmems.nChans;
T0mems = [T0mems; hdrmems.T0];
%read data
data = read_bdf(memsfilename,hdrmems,1,hdrmems.nSamples);
%resample
datarsp = [];
for jj = 1:nbchan
datamems = resample(data(jj,:), EEG.srate, hdrmems.Fs);
datarsp = [datarsp;datamems];
end
if filelength == 1
mems.data = datarsp;
else
memsData{ii} = datarsp;
end
end
if filelength > 1
%sort the mems data files based on start time (T0)
for aa=1:size(T0mems,1)
startsecs(aa) = T0mems(aa,1)*946080000+T0mems(aa,2)*2592000+T0mems(aa,3)*86400+T0mems(aa,4)*3600+T0mems(aa,5)*60+T0mems(aa,6);
end
[startsecs,indexsorted] = sort(startsecs);
for bb = 1:filelength
mems.data = [mems.data,memsData{indexsorted(bb)}];
end
end
%make sure EEG and mems data have same length
nbpnts = size(mems.data,2) ;
if nbpnts == EEG.pnts
EEG.data = [EEG.data; mems.data];
elseif nbpnts > EEG.pnts
EEG.data = [EEG.data; mems.data(:,1:EEG.pnts)];
elseif nbpnts < EEG.pnts
EEG.data = [EEG.data(:,1:hdr.nSamples); mems.data];
end
EEG.chanlocs = [EEG.chanlocs ; chanlocs ];
EEG.nbchan = EEG.nbchan + nbchan;
end
end
else
error('Please select EEG datasets with the filenames starting by "data" ')
end
function dat = read_bdf(filename, hdr, begsample, endsample, chanindx)
if nargin==1
% read the header, this code is from EEGLAB's openbdf
FILENAME = filename;
% defines Seperator for Subdirectories
SLASH='/';
BSLASH=char(92);
cname=computer;
if cname(1:2)=='PC' SLASH=BSLASH; end;
fid=fopen(FILENAME,'r','ieee-le');
if fid<0
fprintf(2,['Error LOADEDF: File ' FILENAME ' not found\n']);
return;
end;
EDF.FILE.FID=fid;
EDF.FILE.OPEN = 1;
EDF.FileName = FILENAME;
PPos=min([max(find(FILENAME=='.')) length(FILENAME)+1]);
SPos=max([0 find((FILENAME=='/') | (FILENAME==BSLASH))]);
EDF.FILE.Ext = FILENAME(PPos+1:length(FILENAME));
EDF.FILE.Name = FILENAME(SPos+1:PPos-1);
if SPos==0
EDF.FILE.Path = pwd;
else
EDF.FILE.Path = FILENAME(1:SPos-1);
end;
EDF.FileName = [EDF.FILE.Path SLASH EDF.FILE.Name '.' EDF.FILE.Ext];
H1=char(fread(EDF.FILE.FID,256,'char')'); %
EDF.FILETYPE = H1(1); % 1 Byte bdf or edf
hdr.filetype = str2num(EDF.FILETYPE);
EDF.VERSION=H1(1:8); % 8 Byte Versionsnummer
%if 0 fprintf(2,'LOADEDF: WARNING Version EDF Format %i',ver); end;
EDF.PID = deblank(H1(9:88)); % 80 Byte local patient identification
EDF.RID = deblank(H1(89:168)); % 80 Byte local recording identification
%EDF.H.StartDate = H1(169:176); % 8 Byte
%EDF.H.StartTime = H1(177:184); % 8 Byte
EDF.T0=[str2num(H1(168+[7 8])) str2num(H1(168+[4 5])) str2num(H1(168+[1 2])) str2num(H1(168+[9 10])) str2num(H1(168+[12 13])) str2num(H1(168+[15 16])) ];
hdr.T0 = EDF.T0;
% Y2K compatibility until year 2090
if EDF.VERSION(1)=='0'
if EDF.T0(1) < 91
EDF.T0(1)=2000+EDF.T0(1);
else
EDF.T0(1)=1900+EDF.T0(1);
end;
else ;
% in a future version, this is hopefully not needed
end;
EDF.HeadLen = str2num(H1(185:192)); % 8 Byte Length of Header
% reserved = H1(193:236); % 44 Byte
EDF.NRec = str2num(H1(237:244)); % 8 Byte # of data records
EDF.Dur = str2num(H1(245:252)); % 8 Byte # duration of data record in sec
EDF.NS = str2num(H1(253:256)); % 8 Byte # of channels
EDF.Label = char(fread(EDF.FILE.FID,[16,EDF.NS],'char')'); %labels of the channels
EDF.Transducer = char(fread(EDF.FILE.FID,[80,EDF.NS],'char')');
EDF.PhysDim = char(fread(EDF.FILE.FID,[8,EDF.NS],'char')');
EDF.PhysMin= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
EDF.PhysMax= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
EDF.DigMin = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
EDF.DigMax = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
% check validity of DigMin and DigMax
if (length(EDF.DigMin) ~= EDF.NS)
fprintf(2,'Warning OPENEDF: Failing Digital Minimum\n');
EDF.DigMin = -(2^15)*ones(EDF.NS,1);
end
if (length(EDF.DigMax) ~= EDF.NS)
fprintf(2,'Warning OPENEDF: Failing Digital Maximum\n');
EDF.DigMax = (2^15-1)*ones(EDF.NS,1);
end
if (any(EDF.DigMin >= EDF.DigMax))
fprintf(2,'Warning OPENEDF: Digital Minimum larger than Maximum\n');
end
% check validity of PhysMin and PhysMax
if (length(EDF.PhysMin) ~= EDF.NS)
fprintf(2,'Warning OPENEDF: Failing Physical Minimum\n');
EDF.PhysMin = EDF.DigMin;
end
if (length(EDF.PhysMax) ~= EDF.NS)
fprintf(2,'Warning OPENEDF: Failing Physical Maximum\n');
EDF.PhysMax = EDF.DigMax;
end
if (any(EDF.PhysMin >= EDF.PhysMax))
fprintf(2,'Warning OPENEDF: Physical Minimum larger than Maximum\n');
EDF.PhysMin = EDF.DigMin;
EDF.PhysMax = EDF.DigMax;
end
EDF.PreFilt= char(fread(EDF.FILE.FID,[80,EDF.NS],'char')'); %
tmp = fread(EDF.FILE.FID,[8,EDF.NS],'char')'; % samples per data record
EDF.SPR = str2num(char(tmp)); % samples per data record
fseek(EDF.FILE.FID,32*EDF.NS,0);
EDF.Cal = (EDF.PhysMax-EDF.PhysMin)./(EDF.DigMax-EDF.DigMin);
EDF.Off = EDF.PhysMin - EDF.Cal .* EDF.DigMin;
tmp = find(EDF.Cal < 0);
EDF.Cal(tmp) = ones(size(tmp));
EDF.Off(tmp) = zeros(size(tmp));
EDF.Calib=[EDF.Off';(diag(EDF.Cal))];
%EDF.Calib=sparse(diag([1; EDF.Cal]));
%EDF.Calib(1,2:EDF.NS+1)=EDF.Off';
EDF.SampleRate = EDF.SPR / EDF.Dur;
EDF.FILE.POS = ftell(EDF.FILE.FID);
if EDF.NRec == -1 % unknown record size, determine correct NRec
fseek(EDF.FILE.FID, 0, 'eof');
endpos = ftell(EDF.FILE.FID);
EDF.NRec = floor((endpos - EDF.FILE.POS) / (sum(EDF.SPR) * 2));
fseek(EDF.FILE.FID, EDF.FILE.POS, 'bof');
H1(237:244)=sprintf('%-8i',EDF.NRec); % write number of records
end;
EDF.Chan_Select=(EDF.SPR==max(EDF.SPR));
for k=1:EDF.NS
if EDF.Chan_Select(k)
EDF.ChanTyp(k)='N';
else
EDF.ChanTyp(k)=' ';
end;
if findstr(upper(EDF.Label(k,:)),'ECG')
EDF.ChanTyp(k)='C';
elseif findstr(upper(EDF.Label(k,:)),'EKG')
EDF.ChanTyp(k)='C';
elseif findstr(upper(EDF.Label(k,:)),'EEG')
EDF.ChanTyp(k)='E';
elseif findstr(upper(EDF.Label(k,:)),'EOG')
EDF.ChanTyp(k)='O';
elseif findstr(upper(EDF.Label(k,:)),'EMG')
EDF.ChanTyp(k)='M';
end;
end;
EDF.AS.spb = sum(EDF.SPR); % Samples per Block
hdr.Annotation = 0;
hdr.AnnotationChn = -1;
if any(EDF.SampleRate~=EDF.SampleRate(1))
chn_indx = find(EDF.SampleRate~=EDF.SampleRate(1));
if(length(chn_indx) > 0 && (strcmp(EDF.Label(chn_indx,:),repmat('BDF Annotations',length(chn_indx),1))) || strcmp(EDF.Label(chn_indx,:),repmat('BDF Annotations ',length(chn_indx),1)))
disp('BDF+ file format detected, the BDF Annotation channel will be skipped when reading data');
hdr.Annotation = 1;
hdr.AnnotationChn = chn_indx;
EDF.SampleRateOrg = EDF.SampleRate;
EDF.SampleRate(chn_indx) = EDF.SampleRate(1);%force it to be the same as others, in order to pass the FieldTrip test
EDF.Label(chn_indx,:) = [];%delete the info about the Annotation channel
EDF.NS = EDF.NS - length(chn_indx);
else
error('Channels with different sampling rate but not BDF Annotation found');
end
end
%read all the events as well
if(hdr.Annotation)
epoch_total = EDF.Dur * sum(EDF.SampleRateOrg); %the number of total data per record, usually one second
% allocate memory to hold the data
n_annotation = sum(EDF.Dur*EDF.SampleRateOrg(hdr.AnnotationChn)); %the number of annotation per record
dat_event = zeros(EDF.NRec,n_annotation*3);
% read and concatenate all required data epochs
for i=1:EDF.NRec
%offset: the bit index right before the current first event
offset = EDF.HeadLen + (i-1)*epoch_total*3 + EDF.Dur*sum(EDF.SampleRateOrg(1:EDF.NS))*3;
%extract the bits for events
buf = readEvents(filename, offset, n_annotation); % see below in subfunction
% buf = readLowLevel(filename, offset, n_annotation);
dat_event(i,:) = buf;
end
%decoding the events
event_cnt = 1;
event = [];
for i = 1:EDF.NRec %last index
char20_index = find(dat_event(i,:)==20);
TAL_start = char20_index(2) + 1;%the first event right after the second char 20
char20_index = char20_index(3:end);%remove the first two 20 (belonging to the TAL time index)
num_event = length(char20_index)/2;%two char 20 per event
if(num_event == 0)
continue;
end
for j = 1:num_event
%structure of a single event: [offset in sec] [char21][Duration in sec] [char20][Trigger Code] [char20][char0]
%the field [char21][Duration in sec] can be skipped
%multiple events can be stored within one Annotation block, the last ends with [char0][char0]
singleTAL = dat_event(i,TAL_start:char20_index(2*j)-1);
char21_one = find(singleTAL == 21);
char20_one = find(singleTAL == 20);
event{event_cnt}.eventtype = 'trigger';%fixed info
event{event_cnt}.eventvalue = char(singleTAL(char20_one+1:end));%trigger code, char
if(~isempty(char21_one))%if Duration field exist
event{event_cnt}.offset_in_sec = str2num(char(singleTAL(1:char21_one-1)));%in sec
event{event_cnt}.duration = str2num(char(singleTAL(char21_one+1:char20_one-1)));%in sec;
else%if no Duration field
event{event_cnt}.offset_in_sec = str2num(char(singleTAL(1:char20_one-1)));%in sec
event{event_cnt}.duration = 0;
end
event{event_cnt}.offset = round(event{event_cnt}.offset_in_sec*EDF.SampleRateOrg(1));%in sampling points
event_cnt = event_cnt + 1;
TAL_start = char20_index(2*j) + 1;
if(dat_event(i,char20_index(2*j)+1) ~= 0)%check if the TAL is complete
error('BDF+ event error');
end
end
%check if all TALs are read
if ~(dat_event(i,char20_index(2*num_event)+1) == 0 && dat_event(i,char20_index(2*num_event)+2) == 0)
error('BDF+ final event error');
end
end
hdr.event = event;
end
% close the file
fclose(EDF.FILE.FID);
% convert the header to Fieldtrip-style
hdr.Fs = EDF.SampleRate(1);
hdr.nChans = EDF.NS;
hdr.label = cellstr(EDF.Label);
hdr.chanlocs = struct([]);
for i = 1:hdr.nChans
chanlocs = struct('labels',hdr.label(i),'ref',[],'theta',[],'radius',[],'X',[],'Y',[],'Z',[],'sph_theta',[],'sph_phi',[],'sph_radius',[],'type',[],'urchan',[] );
hdr.chanlocs = [hdr.chanlocs;chanlocs];
end
% it is continuous data, therefore append all records in one trial
hdr.nTrials = 1;
hdr.nSamples = EDF.NRec * EDF.Dur * EDF.SampleRate(1);
hdr.nSamplesPre = 0;
hdr.orig = EDF;
% return the header
dat = hdr;
else
% read the data
% retrieve the original header
EDF = hdr.orig;
% determine the trial containing the begin and end sample
epochlength = EDF.Dur * EDF.SampleRate(1);
begepoch = floor((begsample-1)/epochlength) + 1;
endepoch = floor((endsample-1)/epochlength) + 1;
nepochs = endepoch - begepoch + 1;
nchans = EDF.NS;
if(hdr.Annotation)
epoch_total = sum(EDF.Dur * EDF.SampleRateOrg);
else
epoch_total = EDF.Dur * EDF.SampleRate(1)*nchans;
end
epoch_data = EDF.Dur * EDF.SampleRate(1)*nchans;
%TODO HERE
if nargin<5
chanindx = 1:nchans;
end
if(hdr.Annotation)
for tmp_indx = 1:length(chanindx)
%make this revision as the Annotation channel is invisible to the
%users, hereby channel index (by user) larger than the Annotation
%channel should increase their index by 1
if(chanindx(tmp_indx) >= hdr.AnnotationChn)
chanindx(tmp_indx) = chanindx(tmp_indx) + 1;
end
end
end
% allocate memory to hold the data
dat = zeros(length(chanindx),nepochs*epochlength);
% read and concatenate all required data epochs
for i=begepoch:endepoch%number of records
if hdr.filetype == 0
offset = EDF.HeadLen + (i-1)*epoch_total*2;
else
offset = EDF.HeadLen + (i-1)*epoch_total*3;
end
% read the data from all channels and then select the desired channels
buf = readLowLevel(filename, offset, epoch_data,hdr.filetype); % see below in subfunction
buf = reshape(buf, epochlength, nchans);
dat(:,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf(:,chanindx)';
end
% select the desired samples
begsample = begsample - (begepoch-1)*epochlength; % correct for the number of bytes that were skipped
endsample = endsample - (begepoch-1)*epochlength; % correct for the number of bytes that were skipped
dat = dat(:, begsample:endsample);
% Calibrate the data
calib = diag(EDF.Cal(chanindx));
if length(chanindx)>1
% using a sparse matrix speeds up the multiplication
dat = sparse(calib) * dat;
else
% in case of one channel the sparse multiplication would result in a sparse array
dat = calib * dat;
end
end
% SUBFUNCTION for reading the 24 bit values
function buf = readLowLevel(filename, offset, numwords,filetype)
% if offset < 2*1024^3
% % use the external mex file, only works for <2GB
% buf = read_24bit(filename, offset, numwords);
% % this would be the only difference between the bdf and edf implementation
% % buf = read_16bit(filename, offset, numwords);
% else
% use plain matlab, thanks to Philip van der Broek
fp = fopen(filename,'r','ieee-le');
status = fseek(fp, offset, 'bof');
if status
error(['failed seeking ' filename]);
end
if filetype == 0
[buf,num] = fread(fp,numwords,'bit16=>double'); %edf
else
[buf,num] = fread(fp,numwords,'bit24=>double'); %bdf
end
fclose(fp);
if (num<numwords)
error(['failed opening ' filename]);
return
% end
end
function buf = readEvents(filename, offset, numwords)
% use plain matlab, thanks to Philip van der Broek
fp = fopen(filename,'r','ieee-le');
status = fseek(fp, offset, 'bof');
if status
error(['failed seeking ' filename]);
end
[buf,num] = fread(fp,numwords*3,'uint8=>char');
fclose(fp);
if (num<numwords)
error(['failed opening ' filename]);
return
end
Matlab读取cnf文件方法
% loadcnt() - Load a Neuroscan continuous signal file.
%
% Usage:
% >> cnt = loadcnt(file, varargin)
%
% Inputs:
% filename - name of the file with extension
%
% Optional inputs:
% 't1' - start at time t1, default 0. Warning, events latency
% might be innacurate (this is an open issue).
% 'sample1' - start at sample1, default 0, overrides t1. Warning,
% events latency might be innacurate.
% 'lddur' - duration of segment to load, default = whole file
% 'ldnsamples' - number of samples to load, default = whole file,
% overrides lddur
% 'scale' - ['on'|'off'] scale data to microvolt (default:'on')
% 'dataformat' - ['int16'|'int32'] default is 'int16' for 16-bit data.
% Use 'int32' for 32-bit data.
% 'blockread' - [integer] by default it is automatically determined
% from the file header, though sometimes it finds an
% incorect value, so you may want to enter a value manually
% here (1 is the most standard value).
% 'memmapfile' - ['memmapfile_name'] use this option if the .cnt file
% is too large to read in conventially. The suffix of
% the memmapfile_name must be .fdt. The memmapfile
% functions process files based on their suffix, and an
% error will occur if you use a different suffix.
%
% Outputs:
% cnt - structure with the continuous data and other informations
% cnt.header
% cnt.electloc
% cnt.data
% cnt.tag
%
% Authors: Sean Fitzgibbon, Arnaud Delorme, 2000-
%
% Note: function original name was load_scan41.m
%
% Known limitations:
% For more see http://www.cnl.salk.edu/~arno/cntload/index.html
% Copyright (C) 2000 Sean Fitzgibbon, <psspf@id.psy.flinders.edu.au>
% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
function [f,lab,ev2p] = loadcnt(filename,varargin)
if ~isempty(varargin)
r=struct(varargin{:});
else r = [];
end;
try, r.t1; catch, r.t1=0; end
try, r.sample1; catch, r.sample1=[]; end
try, r.lddur; catch, r.lddur=[]; end
try, r.ldnsamples; catch, r.ldnsamples=[]; end
try, r.scale; catch, r.scale='on'; end
try, r.blockread; catch, r.blockread = []; end
try, r.dataformat; catch, r.dataformat = 'int32'; end
try, r.memmapfile; catch, r.memmapfile = ''; end
sizeEvent1 = 8 ; %%% 8 bytes for Event1
sizeEvent2 = 19 ; %%% 19 bytes for Event2
sizeEvent3 = 19 ; %%% 19 bytes for Event3
type='cnt';
if nargin ==1
scan=0;
end
fid = fopen(filename,'r','l');
disp(['Loading file ' filename ' ...'])
h.rev = fread(fid,12,'char');
h.nextfile = fread(fid,1,'long');
h.prevfile = fread(fid,1,'ulong');
h.type = fread(fid,1,'char');
h.id = fread(fid,20,'char');
h.oper = fread(fid,20,'char');
h.doctor = fread(fid,20,'char');
h.referral = fread(fid,20,'char');
h.hospital = fread(fid,20,'char');
h.patient = fread(fid,20,'char');
h.age = fread(fid,1,'short');
h.sex = fread(fid,1,'char');
h.hand = fread(fid,1,'char');
h.med = fread(fid,20, 'char');
h.category = fread(fid,20, 'char');
h.state = fread(fid,20, 'char');
h.label = fread(fid,20, 'char');
h.date = fread(fid,10, 'char');
h.time = fread(fid,12, 'char');
h.mean_age = fread(fid,1,'float');
h.stdev = fread(fid,1,'float');
h.n = fread(fid,1,'short');
h.compfile = fread(fid,38,'char');
h.spectwincomp = fread(fid,1,'float');
h.meanaccuracy = fread(fid,1,'float');
h.meanlatency = fread(fid,1,'float');
h.sortfile = fread(fid,46,'char');
h.numevents = fread(fid,1,'int');
h.compoper = fread(fid,1,'char');
h.avgmode = fread(fid,1,'char');
h.review = fread(fid,1,'char');
h.nsweeps = fread(fid,1,'ushort');
h.compsweeps = fread(fid,1,'ushort');
h.acceptcnt = fread(fid,1,'ushort');
h.rejectcnt = fread(fid,1,'ushort');
h.pnts = fread(fid,1,'ushort');
h.nchannels = fread(fid,1,'ushort');
h.avgupdate = fread(fid,1,'ushort');
h.domain = fread(fid,1,'char');
h.variance = fread(fid,1,'char');
h.rate = fread(fid,1,'ushort'); % A USER CLAIMS THAT SAMPLING RATE CAN BE
h.scale = fread(fid,1,'double'); % FRACTIONAL IN NEUROSCAN WHICH IS
h.veogcorrect = fread(fid,1,'char'); % OBVIOUSLY NOT POSSIBLE HERE (BUG 606)
h.heogcorrect = fread(fid,1,'char');
h.aux1correct = fread(fid,1,'char');
h.aux2correct = fread(fid,1,'char');
h.veogtrig = fread(fid,1,'float');
h.heogtrig = fread(fid,1,'float');
h.aux1trig = fread(fid,1,'float');
h.aux2trig = fread(fid,1,'float');
h.heogchnl = fread(fid,1,'short');
h.veogchnl = fread(fid,1,'short');
h.aux1chnl = fread(fid,1,'short');
h.aux2chnl = fread(fid,1,'short');
h.veogdir = fread(fid,1,'char');
h.heogdir = fread(fid,1,'char');
h.aux1dir = fread(fid,1,'char');
h.aux2dir = fread(fid,1,'char');
h.veog_n = fread(fid,1,'short');
h.heog_n = fread(fid,1,'short');
h.aux1_n = fread(fid,1,'short');
h.aux2_n = fread(fid,1,'short');
h.veogmaxcnt = fread(fid,1,'short');
h.heogmaxcnt = fread(fid,1,'short');
h.aux1maxcnt = fread(fid,1,'short');
h.aux2maxcnt = fread(fid,1,'short');
h.veogmethod = fread(fid,1,'char');
h.heogmethod = fread(fid,1,'char');
h.aux1method = fread(fid,1,'char');
h.aux2method = fread(fid,1,'char');
h.ampsensitivity = fread(fid,1,'float');
h.lowpass = fread(fid,1,'char');
h.highpass = fread(fid,1,'char');
h.notch = fread(fid,1,'char');
h.autoclipadd = fread(fid,1,'char');
h.baseline = fread(fid,1,'char');
h.offstart = fread(fid,1,'float');
h.offstop = fread(fid,1,'float');
h.reject = fread(fid,1,'char');
h.rejstart = fread(fid,1,'float');
h.rejstop = fread(fid,1,'float');
h.rejmin = fread(fid,1,'float');
h.rejmax = fread(fid,1,'float');
h.trigtype = fread(fid,1,'char');
h.trigval = fread(fid,1,'float');
h.trigchnl = fread(fid,1,'char');
h.trigmask = fread(fid,1,'short');
h.trigisi = fread(fid,1,'float');
h.trigmin = fread(fid,1,'float');
h.trigmax = fread(fid,1,'float');
h.trigdir = fread(fid,1,'char');
h.autoscale = fread(fid,1,'char');
h.n2 = fread(fid,1,'short');
h.dir = fread(fid,1,'char');
h.dispmin = fread(fid,1,'float');
h.dispmax = fread(fid,1,'float');
h.xmin = fread(fid,1,'float');
h.xmax = fread(fid,1,'float');
h.automin = fread(fid,1,'float');
h.automax = fread(fid,1,'float');
h.zmin = fread(fid,1,'float');
h.zmax = fread(fid,1,'float');
h.lowcut = fread(fid,1,'float');
h.highcut = fread(fid,1,'float');
h.common = fread(fid,1,'char');
h.savemode = fread(fid,1,'char');
h.manmode = fread(fid,1,'char');
h.ref = fread(fid,10,'char');
h.rectify = fread(fid,1,'char');
h.displayxmin = fread(fid,1,'float');
h.displayxmax = fread(fid,1,'float');
h.phase = fread(fid,1,'char');
h.screen = fread(fid,16,'char');
h.calmode = fread(fid,1,'short');
h.calmethod = fread(fid,1,'short');
h.calupdate = fread(fid,1,'short');
h.calbaseline = fread(fid,1,'short');
h.calsweeps = fread(fid,1,'short');
h.calattenuator = fread(fid,1,'float');
h.calpulsevolt = fread(fid,1,'float');
h.calpulsestart = fread(fid,1,'float');
h.calpulsestop = fread(fid,1,'float');
h.calfreq = fread(fid,1,'float');
h.taskfile = fread(fid,34,'char');
h.seqfile = fread(fid,34,'char');
h.spectmethod = fread(fid,1,'char');
h.spectscaling = fread(fid,1,'char');
h.spectwindow = fread(fid,1,'char');
h.spectwinlength = fread(fid,1,'float');
h.spectorder = fread(fid,1,'char');
h.notchfilter = fread(fid,1,'char');
h.headgain = fread(fid,1,'short');
h.additionalfiles = fread(fid,1,'int');
h.unused = fread(fid,5,'char');
h.fspstopmethod = fread(fid,1,'short');
h.fspstopmode = fread(fid,1,'short');
h.fspfvalue = fread(fid,1,'float');
h.fsppoint = fread(fid,1,'short');
h.fspblocksize = fread(fid,1,'short');
h.fspp1 = fread(fid,1,'ushort');
h.fspp2 = fread(fid,1,'ushort');
h.fspalpha = fread(fid,1,'float');
h.fspnoise = fread(fid,1,'float');
h.fspv1 = fread(fid,1,'short');
h.montage = fread(fid,40,'char');
h.eventfile = fread(fid,40,'char');
h.fratio = fread(fid,1,'float');
h.minor_rev = fread(fid,1,'char');
h.eegupdate = fread(fid,1,'short');
h.compressed = fread(fid,1,'char');
h.xscale = fread(fid,1,'float');
h.yscale = fread(fid,1,'float');
h.xsize = fread(fid,1,'float');
h.ysize = fread(fid,1,'float');
h.acmode = fread(fid,1,'char');
h.commonchnl = fread(fid,1,'uchar');
h.xtics = fread(fid,1,'char');
h.xrange = fread(fid,1,'char');
h.ytics = fread(fid,1,'char');
h.yrange = fread(fid,1,'char');
h.xscalevalue = fread(fid,1,'float');
h.xscaleinterval = fread(fid,1,'float');
h.yscalevalue = fread(fid,1,'float');
h.yscaleinterval = fread(fid,1,'float');
h.scaletoolx1 = fread(fid,1,'float');
h.scaletooly1 = fread(fid,1,'float');
h.scaletoolx2 = fread(fid,1,'float');
h.scaletooly2 = fread(fid,1,'float');
h.port = fread(fid,1,'short');
h.numsamples = fread(fid,1,'ulong');
h.filterflag = fread(fid,1,'char');
h.lowcutoff = fread(fid,1,'float');
h.lowpoles = fread(fid,1,'short');
h.highcutoff = fread(fid,1,'float');
h.highpoles = fread(fid,1,'short');
h.filtertype = fread(fid,1,'char');
h.filterdomain = fread(fid,1,'char');
h.snrflag = fread(fid,1,'char');
h.coherenceflag = fread(fid,1,'char');
h.continuoustype = fread(fid,1,'char');
h.eventtablepos = fread(fid,1,'ulong');
h.continuousseconds = fread(fid,1,'float');
h.channeloffset = fread(fid,1,'long');
h.autocorrectflag = fread(fid,1,'char');
h.dcthreshold = fread(fid,1,'uchar');
for n = 1:h.nchannels
e(n).lab = deblank(char(fread(fid,10,'char')'));
e(n).reference = fread(fid,1,'char');
e(n).skip = fread(fid,1,'char');
e(n).reject = fread(fid,1,'char');
e(n).display = fread(fid,1,'char');
e(n).bad = fread(fid,1,'char');
e(n).n = fread(fid,1,'ushort');
e(n).avg_reference = fread(fid,1,'char');
e(n).clipadd = fread(fid,1,'char');
e(n).x_coord = fread(fid,1,'float');
e(n).y_coord = fread(fid,1,'float');
e(n).veog_wt = fread(fid,1,'float');
e(n).veog_std = fread(fid,1,'float');
e(n).snr = fread(fid,1,'float');
e(n).heog_wt = fread(fid,1,'float');
e(n).heog_std = fread(fid,1,'float');
e(n).baseline = fread(fid,1,'short');
e(n).filtered = fread(fid,1,'char');
e(n).fsp = fread(fid,1,'char');
e(n).aux1_wt = fread(fid,1,'float');
e(n).aux1_std = fread(fid,1,'float');
e(n).senstivity = fread(fid,1,'float');
e(n).gain = fread(fid,1,'char');
e(n).hipass = fread(fid,1,'char');
e(n).lopass = fread(fid,1,'char');
e(n).page = fread(fid,1,'uchar');
e(n).size = fread(fid,1,'uchar');
e(n).impedance = fread(fid,1,'uchar');
e(n).physicalchnl = fread(fid,1,'uchar');
e(n).rectify = fread(fid,1,'char');
e(n).calib = fread(fid,1,'float');
end
% finding if 32-bits of 16-bits file
% ----------------------------------
begdata = ftell(fid);
if strcmpi(r.dataformat, 'auto')
r.dataformat = 'int16';
if (h.nextfile > 0)
fseek(fid,h.nextfile+52,'bof');
is32bit = fread(fid,1,'char');
if (is32bit == 1)
r.dataformat = 'int32'
end;
fseek(fid,begdata,'bof');
end;
end;
enddata = h.eventtablepos; % after data
if strcmpi(r.dataformat, 'int16')
nums = (enddata-begdata)/h.nchannels/2;
else nums = (enddata-begdata)/h.nchannels/4;
end;
% number of sample to read
% ------------------------
if ~isempty(r.sample1)
r.t1 = r.sample1/h.rate;
else
r.sample1 = r.t1*h.rate;
end;
if strcmpi(r.dataformat, 'int16')
startpos = r.t1*h.rate*2*h.nchannels;
else startpos = r.t1*h.rate*4*h.nchannels;
end;
if isempty(r.ldnsamples)
if ~isempty(r.lddur)
r.ldnsamples = round(r.lddur*h.rate);
else r.ldnsamples = nums;
end;
end;
% channel offset
% --------------
if ~isempty(r.blockread)
h.channeloffset = r.blockread;
end;
if h.channeloffset > 1
fprintf('WARNING: reading data in blocks of %d, if this fails, try using option "''blockread'', 1"\n', ...
h.channeloffset);
end;
disp('Reading data .....')
if type == 'cnt'
% while (ftell(fid) +1 < h.eventtablepos)
%d(:,i)=fread(fid,h.nchannels,'int16');
%end
fseek(fid, startpos, 0);
% **** This marks the beginning of the code modified for reading
% large .cnt files
% Switched to r.memmapfile for continuity. Check to see if the
% variable exists. If it does, then the user has indicated the
% file is too large to be processed in memory. If the variable
% is blank, the file is processed in memory.
if (~isempty(r.memmapfile))
% open a file for writing
foutid = fopen(r.memmapfile, 'w') ;
% This portion of the routine reads in a section of the EEG file
% and then writes it out to the harddisk.
samples_left = h.nchannels * r.ldnsamples ;
% the size of the data block to be read is limited to 4M
% samples. This equates to 16MB and 32MB of memory for
% 16 and 32 bit files, respectively.
data_block = 4000000 ;
max_rows = data_block / h.nchannels ;
%warning off ;
max_written = h.nchannels * uint32(max_rows) ;
%warning on ;
% This while look tracks the remaining samples. The
% data is processed in chunks rather than put into
% memory whole.
while (samples_left > 0)
% Check to see if the remaining data is smaller than
% the general processing block by looking at the
% remaining number of rows.
to_read = max_rows ;
if (data_block > samples_left)
to_read = samples_left / h.nchannels ;
end ;
% Read data in a relatively small chunk
temp_dat = fread(fid, [h.nchannels to_read], r.dataformat) ;
% The data is then scaled using the original routine.
% In the original routine, the entire data set was scaled
% after being read in. For this version, scaling occurs
% after every chunk is read.
if strcmpi(r.scale, 'on')
disp('Scaling data .....')
%%% scaling to microvolts
for i=1:h.nchannels
bas=e(i).baseline;sen=e(i).senstivity;cal=e(i).calib;
mf=sen*(cal/204.8);
temp_dat(i,:)=(temp_dat(i,:)-bas).*mf;
end
end
% Write out data in float32 form to the file name
% supplied by the user.
written = fwrite (foutid, temp_dat, 'float32') ;
if (written ~= max_written)
samples_left = 0 ;
else
samples_left = samples_left - written ;
end ;
end ;
fclose (foutid) ;
% Set the dat variable. This gets used later by other
% EEGLAB functions.
dat = r.memmapfile ;
% This variable tracks how the data should be read.
bReadIntoMemory = false ;
else
% The memmapfile variable is empty, read into memory.
bReadIntoMemory = true ;
end
% This ends the modifications made to read large files.
% Everything contained within the following if statement is the
% original code.
if (bReadIntoMemory == true)
if h.channeloffset <= 1
dat=fread(fid, [h.nchannels Inf], r.dataformat);
if size(dat,2) < r.ldnsamples
dat=single(dat);
r.ldnsamples = size(dat,2);
else
dat=single(dat(:,1:r.ldnsamples));
end;
else
h.channeloffset = h.channeloffset/2;
% reading data in blocks
dat = zeros( h.nchannels, r.ldnsamples, 'single');
dat(:, 1:h.channeloffset) = fread(fid, [h.channeloffset h.nchannels], r.dataformat)';
counter = 1;
while counter*h.channeloffset < r.ldnsamples
dat(:, counter*h.channeloffset+1:counter*h.channeloffset+h.channeloffset) = ...
fread(fid, [h.channeloffset h.nchannels], r.dataformat)';
counter = counter + 1;
end;
end ;
% ftell(fid)
if strcmpi(r.scale, 'on')
disp('Scaling data .....')
%%% scaling to microvolts
for i=1:h.nchannels
bas=e(i).baseline;sen=e(i).senstivity;cal=e(i).calib;
mf=sen*(cal/204.8);
dat(i,:)=(dat(i,:)-bas).*mf;
end % end for i=1:h.nchannels
end; % end if (strcmpi(r.scale, 'on')
end ;
ET_offset = (double(h.prevfile) * (2^32)) + double(h.eventtablepos); % prevfile contains high order bits of event table offset, eventtablepos contains the low order bits
fseek(fid, ET_offset, 'bof');
disp('Reading Event Table...')
eT.teeg = fread(fid,1,'uchar');
eT.size = fread(fid,1,'ulong');
eT.offset = fread(fid,1,'ulong');
if eT.teeg==2
nevents=eT.size/sizeEvent2;
if nevents > 0
ev2(nevents).stimtype = [];
for i=1:nevents
ev2(i).stimtype = fread(fid,1,'ushort');
ev2(i).keyboard = fread(fid,1,'char');
temp = fread(fid,1,'uint8');
ev2(i).keypad_accept = bitand(15,temp);
ev2(i).accept_ev1 = bitshift(temp,-4);
ev2(i).offset = fread(fid,1,'long');
ev2(i).type = fread(fid,1,'short');
ev2(i).code = fread(fid,1,'short');
ev2(i).latency = fread(fid,1,'float');
ev2(i).epochevent = fread(fid,1,'char');
ev2(i).accept = fread(fid,1,'char');
ev2(i).accuracy = fread(fid,1,'char');
end
else
ev2 = [];
end;
elseif eT.teeg==3 % type 3 is similar to type 2 except the offset field encodes the global sample frame
nevents=eT.size/sizeEvent3;
if nevents > 0
ev2(nevents).stimtype = [];
if r.dataformat == 'int32'
bytes_per_samp = 4; % I only have 32 bit data, unable to check whether this is necessary,
else % perhaps there is no type 3 file with 16 bit data
bytes_per_samp = 2;
end
for i=1:nevents
ev2(i).stimtype = fread(fid,1,'ushort');
ev2(i).keyboard = fread(fid,1,'char');
temp = fread(fid,1,'uint8');
ev2(i).keypad_accept = bitand(15,temp);
ev2(i).accept_ev1 = bitshift(temp,-4);
os = fread(fid,1,'ulong');
ev2(i).offset = os * bytes_per_samp * h.nchannels;
ev2(i).type = fread(fid,1,'short');
ev2(i).code = fread(fid,1,'short');
ev2(i).latency = fread(fid,1,'float');
ev2(i).epochevent = fread(fid,1,'char');
ev2(i).accept = fread(fid,1,'char');
ev2(i).accuracy = fread(fid,1,'char');
end
else
ev2 = [];
end;
elseif eT.teeg==1
nevents=eT.size/sizeEvent1;
if nevents > 0
ev2(nevents).stimtype = [];
for i=1:nevents
ev2(i).stimtype = fread(fid,1,'ushort');
ev2(i).keyboard = fread(fid,1,'char');
% modified by Andreas Widmann 2005/05/12 14:15:00
%ev2(i).keypad_accept = fread(fid,1,'char');
temp = fread(fid,1,'uint8');
ev2(i).keypad_accept = bitand(15,temp);
ev2(i).accept_ev1 = bitshift(temp,-4);
% end modification
ev2(i).offset = fread(fid,1,'long');
end;
else
ev2 = [];
end;
else
disp('Skipping event table (tag != 1,2,3 ; theoritically impossible)');
ev2 = [];
end
fseek(fid, -1, 'eof');
t = fread(fid,'char');
f.header = h;
f.electloc = e;
f.data = dat;
f.Teeg = eT;
f.event = ev2;
f.tag=t;
% Surgical addition of number of samples
f.ldnsamples = r.ldnsamples ;
%%%% channels labels
for i=1:h.nchannels
plab=sprintf('%c',f.electloc(i).lab);
if i>1
lab=str2mat(lab,plab);
else
lab=plab;
end
end
%%%% to change offest in bytes to points
if ~isempty(ev2)
if r.sample1 ~= 0
fprintf(2,'Warning: events imported with a time shift might be innacurate\n');
end;
ev2p=ev2;
ioff=900+(h.nchannels*75); %% initial offset : header + electordes desc
if strcmpi(r.dataformat, 'int16')
for i=1:nevents
ev2p(i).offset=(ev2p(i).offset-ioff)/(2*h.nchannels) - r.sample1; %% 2 short int end
end
else % 32 bits
for i=1:nevents
ev2p(i).offset=(ev2p(i).offset-ioff)/(4*h.nchannels) - r.sample1; %% 4 short int end
end
end;
f.event = ev2p;
end;
frewind(fid);
fclose(fid);
endtdms文件
function [ConvertedData,ConvertVer,ChanNames,GroupNames,ci]=convertTDMS(varargin)
%Function to load LabView TDMS data file(s) into variables in the MATLAB workspace.
%An *.MAT file can also be created. If called with one input, the user selects
%a data file.
%
% TDMS format is based on information provided by National Instruments at:
% http://zone.ni.com/devzone/cda/tut/p/id/5696
%
% [ConvertedData,ConvertVer,ChanNames]=convertTDMS(SaveConvertedFile,filename)
%
% Inputs:
% SaveConvertedFile (required) - Logical flag (true/false) that
% determines whether a MAT file is created. The MAT file's name
% is the same as 'filename' except that the 'TDMS' file extension is
% replaced with 'MAT'. The MAT file is saved in the same folder
% and will overwrite an existing file without warning. The
% MAT file contains all the output variables.
%
% filename (optional) - Filename (fully defined) to be converted.
% If not supplied, the user is provided a 'File Open' dialog box
% to navigate to a file. Can be a cell array of files for bulk
% conversion.
%
% Outputs:
% ConvertedData (required) - Structure of all of the data objects.
% ConvertVer (optional) - Version number of this function.
% ChanNames (optional) - Cell array of channel names
% GroupNames (optional) - Cell array of group names
% ci (optional) - Structure of the channel index (an index to
% where all of the information for a channel resides in a
% file.
%
%
%'ConvertedData' is a structure with 'FileName', 'FileFolder', 'SegTDMSVerNum',
%'NumOfSegments' and 'Data' fields'. The 'Data' field is a structure.
%
%'ConvertedData.SegTDMSVerNum' is a vector of the TDMS version number for each
%segment.
%
%'ConvertedData.Data' is a structure with 'Root' and 'MeasuredData' fields.
%
%'ConvertedData.Data.Root' is a structure with 'Name' and 'Property' fields.
%The 'Property' field is also a structure; it contains all the specified properties
%(1 entry for each 'Property) for the 'Root' group. For each 'Property' there are
%'Name' and 'Value' fields. To display a list of all the property names, input
%'{ConvertedData.Data.Root.Property.Name}'' in the Command Window.
%
%'ConvertedData.Data.MeasuredData' is a structure containing all the channel/group
%information. For each index (for example, 'ConvertedData.Data.MeasuredData(1)'),
%there are 'Name', 'Data' and 'Property' fields. The list of channel names can
%be displayed by typing 'ChanNames' in the Command Window. Similarly, the list
%of group names can be displayed by typing 'GroupNames' in the Command Window.
%The 'Property' field is also a structure; it contains all the specified properties
%for that index (1 entry in the structure for each 'Property'). Any LabView waveform
%attributes ('wf_start_time', 'wf_start_offset', 'wf_increment' and 'wf_samples') that
%may exist are also included in the properties. For each 'Property' there are 'Name'
%and 'Value' fields. To display a list of all the property names, input
%'{ConvertedData.Data.MeasuredData(#).Property.Name}'' in the Command Window
%where '#' is the index of interest.
%
%If you recieve an error that DAQmxRaw data cannot be converted, this is due
%to the details on parsing this data type not being published by NI. the
%work around is to use a VI that reads in the data file and then writes it
%to a new file. See: https://decibel.ni.com/content/docs/DOC-32817
%
% See Also: simpleconvertTDMS
%-------------------------------------------------------------------------
%Brad Humphreys - v1.0 2008-04-23
%ZIN Technologies
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Brad Humphreys - v1.1 2008-07-03
%ZIN Technologies
%-Added abilty for timestamp to be a raw data type, not just meta data.
%-Addressed an issue with having a default nsmaples entry for new objects.
%-Added Error trap if file name not found.
%-Corrected significant problem where it was assumed that once an object
% existsed, it would in in every subsequent segement. This is not true.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Grant Lohsen - v1.2 2009-11-15
%Georgia Tech Research Institute
%-Converts TDMS v2 files
%Folks, it's not pretty but I don't have time to make it pretty. Enjoy.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Jeff Sitterle - v1.3 2010-01-10
%Georgia Tech Research Institute
%Modified to return all information stored in the TDMS file to inlcude
%name, start time, start time offset, samples per read, total samples, unit
%description, and unit string. Also provides event time and event
%description in text form
%Vast speed improvement as save was the previous longest task
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Grant Lohsen - v1.4 2009-04-15
%Georgia Tech Research Institute
%Reads file header info and stores in the Root Structure.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Robert Seltzer - v1.5 2010-07-14
%BorgWarner Morse TEC
%-Tested in MATLAB 2007b and 2010a.
%-APPEARS to now be compatible with TDMS version 1.1 (a.k.a 4712) files;
% although, this has not been extensively tested. For some unknown
% reason, the version 1.2 (4713) files process noticeably faster. I think
% that it may be related to the 'TDSm' tag.
%-"Time Stamp" data type was not tested.
%-"Waveform" fields was not tested.
%-Fixed an error in the 'LV2MatlabDataType' function where LabView data type
% 'tdsTypeSingleFloat' was defined as MATLAB data type 'float64' . Changed
% to 'float32'.
%-Added error trapping.
%-Added feature to count the number of segments for pre-allocation as
% opposed to estimating the number of segments.
%-Added option to save the data in a MAT file.
%-Fixed "invalid field name" error caused by excessive string lengths.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Robert Seltzer - v1.6 2010-09-01
%BorgWarner Morse TEC
%-Tested in MATLAB 2010a.
%-Fixed the "Coversion to cell from char is not possible" error found
% by Francisco Botero in version 1.5.
%-Added capability to process both fragmented or defragmented data.
%-Fixed the "field" error found by Lawrence.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Christian Buxel - V1.7 2010-09-17
%RWTH Aachen
%-Tested in Matlab2007b.
%-Added support for german umlauts (�,�,�,�,�,�,�) in 'propsName'
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Andr� R�egg - V1.7 2010-09-29
%Supercomputing Systems AG
%-Tested in MATLAB 2006a & 2010b
%-Make sure that data can be loaded correctly independently of character
% encoding set in matlab.
%-Fixed error if object consists of several segments with identical segment
% information (if rawdataindex==0, not all segments were loaded)
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Robert Seltzer - v1.7 2010-09-30
%BorgWarner Morse TEC
%-Tested in MATLAB 2010b.
%-Added 'error trapping' to the 'fixcharformatlab' function for group and
% channel names that contain characters that are not 'A' through 'Z',
% 'a' through 'z', 0 through 9 or underscore. The characters are replaced
% with an underscore and a message is written to the Command Window
% explaining to the user what has happened and how to fix it. Only tested
% with a very limited number of "special" characters.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Robert Seltzer - v1.8 2010-10-12
%BorgWarner Morse TEC
%-As a result of an error found by Peter Sulcs when loading data with very
% long channel names, I have re-written the sections of the function that
% creates the channel and property names that are used within the body of
% the function to make them robust against long strings and strings
% containing non-UTF8 characters. The original channel and property
% names (no truncation or character replacement) are now retained and
% included in the output structure. In order to implement this improvement,
% I added a 'Property' field as a structure to the 'ConvertedData' output
% structure.
%-Added a more detailed 'help' description ('doc convertTDMS') of the
% returned structure.
%-List of channel names added as an output parameter of the function.
%-Corrected an error in the time stamp converion. It was off by exactly
% 1 hour.
%-Tested in MATLAB 2010b with a limited number of diverse TDMS files.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Robert Seltzer - v1.8 2010-10-19
%BorgWarner Morse TEC
%-Fixed an error found by Terenzio Girotto with the 'save' routine.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Robert Seltzer - v1.8 2010-10-25
%BorgWarner Morse TEC
%-Fixed an error with channels that contain no data. Previously, if a
% channel contained no data, then it was not passed to the output structure
% even if it did contain properties.
%-Added 'GroupNames' as an optional output variable.
%-Fixed an error with capturing the properties of the Root object
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Philip Top - v1.9 2010-11-09
%John Breneman
%-restructured code as function calls
%-seperated metadata reads from data reads
%-preallocated space for SegInfo with two pass file read
%-preallocated index information and defined segdataoffset for each segment
%-preallocate space for data for speedup in case of fragmented files
%-used matlab input parser instead of nargin switch
%-vectorized timestamp reads for substantial speedup
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Robert Seltzer - v1.9 2010-11-10
%BorgWarner Morse TEC
%-Fixed an error error in the 'offset' calculation for strings
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Philip Top - v1.95 2011-5-10
%Fix Bug with out of order file segments
%Fix some issues with string array reads for newer version files,
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Brad Humphreys - v1.96 2013-10-9
%Fixed problem with error catch messages themselves (interleaved, version,
%big endian) using the TDMSFileName variable which was not available
%(passed into) the getSegInfo function. Added ability to work with
%interleaved and big endian files.
%So function now covers:
% -v1.0-v2.0
% -Interleaved and Decimated Data Formats
% -Big and Little Endian storage
% Does not work with DAQmxRawData
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Brad Humphreys - v1.97 2013-11-13
%Added information to help documentation on how to deal with DAQmxRaw Data.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Brad Humphreys - v1.98 2014-5-27
%Per G. Lohsen's suggestion, added check to verify that first caharters are
% TDMs. If not, errors out and lets user know that the selected file is
% not a TDMS file.
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%Sebastian Schwarzendahl (alias Haench) - v1.99 2014-10-23
% Added support for complex data types
% CSG - tdsTypeComplexSingleFloat=0x08000c
% CDB - tdsTypeComplexDoubleFloat=0x10000d)
% This feature was added in LV2013 (I believ) and produced an error in
% the previous version of this code.
%-------------------------------------------------------------------------
%Initialize outputs
ConvertVer='1.98'; %Version number of this conversion function
ConvertedData=[];
p=inputParser();
p.addRequired('SaveConvertedFile',@(x) islogical(x)||(ismember(x,[0,1])));
p.addOptional('filename','',@(x) iscell(x)||exist(x,'file'));
p.parse(varargin{:});
filename=p.Results.filename;
SaveConvertedFile=p.Results.SaveConvertedFile;
if isempty(filename)
%Prompt the user for the file
[filename,pathname]=uigetfile({'*.tdms','All Files (*.tdms)'},'Choose a TDMS File');
if filename==0
return
end
filename=fullfile(pathname,filename);
end
if iscell(filename)
%For a list of files
infilename=filename;
else
infilename=cellstr(filename);
end
for fnum=1:numel(infilename)
if ~exist(infilename{fnum},'file')
e=errordlg(sprintf('File ''%s'' not found.',infilename{fnum}),'File Not Found');
uiwait(e)
return
end
FileNameLong=infilename{fnum};
[pathstr,name,ext]=fileparts(FileNameLong);
FileNameShort=sprintf('%s%s',name,ext);
FileNameNoExt=name;
FileFolder=pathstr;
if fnum==1
fprintf('\n\n')
end
fprintf('Converting ''%s''...',FileNameShort)
fid=fopen(FileNameLong);
if fid==-1
e=errordlg(sprintf('Could not open ''%s''.',FileNameLong),'File Cannot Be Opened');
uiwait(e)
fprintf('\n\n')
return
end
% Build a database with segment info
[SegInfo,NumOfSeg]=getSegInfo(fid);
% Build a database with channel info
[channelinfo SegInfo]=getChannelInfo(fid,SegInfo,NumOfSeg);
% Add channel count to SegInfo
%SegInfo=addChanCount(SegInfo,channelinfo);
% Get the raw data
ob=getData(fid,channelinfo,SegInfo); %Returns the objects which have data. See postProcess function (appends to all of the objects)
fclose(fid);
%Assign the outputs
ConvertedData(fnum).FileName=FileNameShort;
ConvertedData(fnum).FileFolder=FileFolder;
ConvertedData(fnum).SegTDMSVerNum=SegInfo.vernum;
ConvertedData(fnum).NumOfSegments=NumOfSeg;
[ConvertedData(fnum).Data,CurrGroupNames]=postProcess(ob,channelinfo);
GroupNames(fnum)={CurrGroupNames};
TempChanNames={ConvertedData(fnum).Data.MeasuredData.Name};
TempChanNames(strcmpi(TempChanNames,'Root'))=[];
ChanNames(fnum)={sort(setdiff(TempChanNames',CurrGroupNames))};
if SaveConvertedFile
MATFileNameShort=sprintf('%s.mat',FileNameNoExt);
MATFileNameLong=fullfile(FileFolder,MATFileNameShort);
try
save(MATFileNameLong,'ConvertedData','ConvertVer','ChanNames')
fprintf('\n\nConversion complete (saved in ''%s'').\n\n',MATFileNameShort)
catch exception
fprintf('\n\nConversion complete (could not save ''%s'').\n\t%s: %s\n\n',MATFileNameShort,exception.identifier,...
exception.message)
end
else
fprintf('\n\nConversion complete.\n\n')
end
end
ci=channelinfo;
end
%%
function [SegInfo,NumOfSeg]=getSegInfo(fid)
% Go through the whole file and build a database of the segment lead-in
% information:
%1) Count the total number of segments in the file. (NumOfSeg)
%2) Get all of the "Lead In" Headers from all of the segements and store
% that in the SegInfo variable:
% SegInfo.SegStartPosn: Absolute Starting Position of Segment in file
% SegInfo.MetaStartPosn: Absolute Starting Position of Meta Data in file
% SegInfo.DataStartPosn: Absolute Starting Position of Raw Data in file
% SegInfo.DataLength: Number of bytes of Data in segment
% SegInfo.vernum: LV version number (4712 is v1.0, 4713 is v2.0)
% SegInfo.NumChan: number of channels in this segement. This is
% only instatated in this function to 0. The addChanCount function
% updates this to the actual value later.
%
% SegInfo.SegHasMetaData: There is Meta Data in the Segement
% SegInfo.SegHasRawData: There is Raw Data in the Segment
% SegInfo.SegHasDaqMxRaw: There is DAQmxRaw Data in the Segment
% SegInfo.SegInterleaved: 0: Contigous Data, 1: Interleaved Data
% SegInfo.SegBigEndian: 0: Little, 1:Big (numeric, leadin, raw, and meta)
% SegInfo.SegHasNewObjList: New object list in Segment
%
%3) While doing the above, also include error trapping for incompatibity
%Find the end of the file
fseek(fid,0,'eof');
eoff=ftell(fid);
frewind(fid);
segCnt=0;
CurrPosn=0;
LeadInByteCount=28; %From the National Instruments web page (http://zone.ni.com/devzone/cda/tut/p/id/5696) under
%the 'Lead In' description on page 2: Counted the bytes shown in the table.
while (ftell(fid) ~= eoff)
Ttag=fread(fid,1,'uint8','l');
Dtag=fread(fid,1,'uint8','l');
Stag=fread(fid,1,'uint8','l');
mtag=fread(fid,1,'uint8','l');
if Ttag==84 && Dtag==68 && Stag==83 && mtag==109
%Apparently, this sequence of numbers identifies the start of a new segment.
segCnt=segCnt+1;
%ToC Field
ToC=fread(fid,1,'uint32','l');
kTocBigEndian=bitget(ToC,7);
if kTocBigEndian
kTocEndian='b';
else
kTocEndian='l';
end
%TDMS format version number
vernum=fread(fid,1,'uint32',kTocEndian);
%From the National Instruments web page (http://zone.ni.com/devzone/cda/tut/p/id/5696) under the 'Lead In'
%description on page 2:
%The next eight bytes (64-bit unsigned integer) describe the length of the remaining segment (overall length of the
%segment minus length of the lead in). If further segments are appended to the file, this number can be used to
%locate the starting point of the following segment. If an application encountered a severe problem while writing
%to a TDMS file (crash, power outage), all bytes of this integer can be 0xFF. This can only happen to the last
%segment in a file.
nlen=fread(fid,1,'uint64',kTocEndian);
if (nlen>2^63)
break;
else
segLength=nlen;
end
TotalLength=segLength+LeadInByteCount;
CurrPosn=CurrPosn+TotalLength;
status=fseek(fid,CurrPosn,'bof'); %Move to the beginning position of the next segment
if (status<0)
warning('file glitch');
break;
end
else %TDSm should be the first charaters in a tdms file. If not there, error out to stop hunting.
fclose(fid);
error('Unable to find TDSm tag. This may not be a tdms file, or you forgot to add the .tdms extension to the filename and are reading the wrong file');
end
end
frewind(fid);
CurrPosn=0;
SegInfo.SegStartPosn=zeros(segCnt,1);
SegInfo.MetaStartPosn=zeros(segCnt,1);
SegInfo.DataStartPosn=zeros(segCnt,1);
SegInfo.vernum=zeros(segCnt,1);
SegInfo.DataLength=zeros(segCnt,1);
SegInfo.NumChan=zeros(segCnt,1);
SegInfo.SegHasMetaData=false(segCnt,1);
SegInfo.SegHasRawData=false(segCnt,1);
SegInfo.SegHasDaqMxRaw=false(segCnt,1);
SegInfo.SegInterleaved=false(segCnt,1);
SegInfo.SegBigEndian=false(segCnt,1);
SegInfo.SegHasNewObjList=false(segCnt,1);
segCnt=0;
while (ftell(fid) ~= eoff)
Ttag=fread(fid,1,'uint8','l');
Dtag=fread(fid,1,'uint8','l');
Stag=fread(fid,1,'uint8','l');
mtag=fread(fid,1,'uint8','l');
if Ttag==84 && Dtag==68 && Stag==83 && mtag==109
%Apparently, this sequence of numbers identifies the start of a new segment.
%Leaving the above comment in, because it reflects that state of
%the NI documenation on TDMS files when the contributers first started
%developing this code.
segCnt=segCnt+1;
if segCnt==1
StartPosn=0;
else
StartPosn=CurrPosn;
end
%ToC Field
ToC=fread(fid,1,'uint32','l');
kTocBigEndian=bitget(ToC,7);
kTocMetaData=bitget(ToC,2);
kTocNewObjectList=bitget(ToC,3);
kTocRawData=bitget(ToC,4);
kTocDaqMxRawData=bitget(ToC,8);
kTocInterleavedData=bitget(ToC,6);
if kTocBigEndian
kTocEndian='b';
else
kTocEndian='l';
end
% if kTocInterleavedData
% error([sprintf(['\n Seqment %.0f of the above file has interleaved data which is not supported with this '...
% 'function. '],segCnt),'Interleaved Data Not Supported']);
% fclose(fid);
% end
% if kTocBigEndian
% error(sprintf(['\n Seqment %.0f of the above file uses the big-endian data format which is not supported '...
% 'with this function. '],segCnt),'Big-Endian Data Format Not Supported');
% fclose(fid);
% end
if kTocDaqMxRawData
error(sprintf(['\n Seqment %.0f of the above file contains data in the DAQmxRaw NI datatype format which is not supported '...
'with this function. See help documentation in convertTDMS.m for how to fix this. '],segCnt),'DAQmxRawData Format Not Supported');
fclose(fid);
end
%TDMS format version number
vernum=fread(fid,1,'uint32',kTocEndian);
if ~ismember(vernum,[4712,4713])
error(sprintf(['\n Seqment %.0f of the above file used LabView TDMS file format version %.0f which is not '...
'supported with this function (%s.m).'],segCnt,vernum),...
'TDMS File Format Not Supported');
fclose(fid);
end
%From the National Instruments web page (http://zone.ni.com/devzone/cda/tut/p/id/5696) under the 'Lead In'
%description on page 2:
%The next eight bytes (64-bit unsigned integer) describe the length of the remaining segment (overall length of the
%segment minus length of the lead in). If further segments are appended to the file, this number can be used to
%locate the starting point of the following segment. If an application encountered a severe problem while writing
%to a TDMS file (crash, power outage), all bytes of this integer can be 0xFF. This can only happen to the last
%segment in a file.
segLength=fread(fid,1,'uint64',kTocEndian);
metaLength=fread(fid,1,'uint64',kTocEndian);
if (segLength>2^63)
fseek(fid,0,'eof');
flen=ftell(fid);
segLength=flen-LeadInByteCount-TotalLength;
TotalLength=segLength+LeadInByteCount;
else
TotalLength=segLength+LeadInByteCount;
CurrPosn=CurrPosn+TotalLength;
fseek(fid,CurrPosn,'bof'); %Move to the beginning position of the next segment
end
SegInfo.SegStartPosn(segCnt)=StartPosn;
SegInfo.MetaStartPosn(segCnt)=StartPosn+LeadInByteCount;
SegInfo.DataStartPosn(segCnt)=SegInfo.MetaStartPosn(segCnt)+metaLength;
SegInfo.DataLength(segCnt)=segLength-metaLength;
SegInfo.vernum(segCnt)=vernum;
SegInfo.SegHasMetaData(segCnt)=kTocMetaData;
SegInfo.SegHasRawData(segCnt)=kTocRawData;
SegInfo.SegHasDaqMxRaw(segCnt)=kTocDaqMxRawData;
SegInfo.SegInterleaved(segCnt)=kTocInterleavedData;
SegInfo.SegBigEndian(segCnt)=kTocBigEndian;
SegInfo.SegHasNewObjList(segCnt)=kTocNewObjectList;
end
end
NumOfSeg=segCnt;
end
%%
function [index SegInfo]=getChannelInfo(fid,SegInfo,NumOfSeg)
%Loop through the segements and get all of the object information.
% name: Short name such as Object3
% long_name: The path directory form of the object name
% rawdatacount: number of segments that have raw data for this object
% datastartindex: Absolute position in file of the begining of all raw
% data for the segment. Add the rawdata offset to get to this object's
% start point.
% arrayDim: array dimension of raw data (for now should always be one as
% NI does not support multi-dimension
% nValues: number of values of raw data in segment
% byteSize: number of bytes for string/char data
% index: segment numbers that contain raw data for this object
% rawdataoffset: offset from datastartindex for the begining of raw data
% for this object
% multiplier: As part of LV's optimization, it will appenend raw data writes to the
% TDMS file when meta data does not change, but forget to update the
% nValues. multiplier is the number of times the base pattern
% repeats.
% skip: When multiple raw data writes are appended, this is the number of
% bytes between the end of a channel in one block to the begining of that channel in
% the next block.
% dataType: The data type for this object. See LV2MatlabDataType function.
% datasize: Number of bytes for each raw data entry.
% Property Structures: Structure containg properties.
%Initialize variables for the file conversion
index=struct();
objOrderList={};
for segCnt=1:NumOfSeg
%Go to the segment starting position of the segment
fseek(fid,SegInfo.SegStartPosn(segCnt)+28,'bof');
% +28 bytes: TDSm (4) + Toc (4) + segVer# (4) + segLength (8) + metaLength (8)
%segVersionNum=SegInfo.vernum(segCnt);
kTocMetaData=SegInfo.SegHasMetaData(segCnt);
kTocNewObjectList=SegInfo.SegHasNewObjList(segCnt);
kTocRawData=SegInfo.SegHasRawData(segCnt);
kTocInterleavedData=SegInfo.SegInterleaved(segCnt);
kTocBigEndian=SegInfo.SegBigEndian(segCnt);
if kTocBigEndian
kTocEndian='b';
else
kTocEndian='l';
end
%% Process Meta Data
%If the object list from the last segment should be used....
if (kTocNewObjectList==0)
fnm=fieldnames(index); %Get a list of the objects/channels to loop through
for kk=1:length(fnm)
ccnt=index.(fnm{kk}).rawdatacount;
if (ccnt>0) %If there is raw data info in the previous segement, copy it into the new segement
if (index.(fnm{kk}).index(ccnt)==segCnt-1)
ccnt=ccnt+1;
index.(fnm{kk}).rawdatacount=ccnt;
index.(fnm{kk}).datastartindex(ccnt)=SegInfo.DataStartPosn(segCnt);
index.(fnm{kk}).arrayDim(ccnt)=index.(fnm{kk}).arrayDim(ccnt-1);
index.(fnm{kk}).nValues(ccnt)=index.(fnm{kk}).nValues(ccnt-1);
index.(fnm{kk}).byteSize(ccnt)=index.(fnm{kk}).byteSize(ccnt-1);
index.(fnm{kk}).index(ccnt)=segCnt;
index.(fnm{kk}).rawdataoffset(ccnt)=index.(fnm{kk}).rawdataoffset(ccnt-1);
SegInfo.NumChan(segCnt)=SegInfo.NumChan(segCnt)+1;
end
end
end
end
%If there is Meta data in the segement
if kTocMetaData
numObjInSeg=fread(fid,1,'uint32',kTocEndian);
if (kTocNewObjectList)
objOrderList=cell(numObjInSeg,1);
end
for q=1:numObjInSeg
obLength=fread(fid,1,'uint32',kTocEndian); %Get the length of the objects name
ObjName=convertToText(fread(fid,obLength,'uint8','l'))'; %Get the objects name
if strcmp(ObjName,'/')
long_obname='Root';
else
long_obname=ObjName;
%Delete any apostrophes. If the first character is a slash (forward or backward), delete it too.
long_obname(strfind(long_obname,''''))=[];
if strcmpi(long_obname(1),'/') || strcmpi(long_obname(1),'\')
long_obname(1)=[];
end
end
newob=0;
%Create object's name. Use a generic field name to avoid issues with strings that are too long and/or
%characters that cannot be used in MATLAB variable names. The actual channel name is retained for the final
%output structure.
if exist('ObjNameList','var')
%Check to see if the object already exists
NameIndex=find(strcmpi({ObjNameList.LongName},long_obname)==1,1,'first');
if isempty(NameIndex)
newob=1;
%It does not exist, so create the generic name field name
ObjNameList(end+1).FieldName=sprintf('Object%.0f',numel(ObjNameList)+1);
ObjNameList(end).LongName=long_obname;
NameIndex=numel(ObjNameList);
end
else
%No objects exist, so create the first one using a generic name field name.
ObjNameList.FieldName='Object1';
ObjNameList.LongName=long_obname;
NameIndex=1;
newob=1;
end
%Assign the generic field name
obname=ObjNameList(NameIndex).FieldName;
%Create the 'index' structure
if (~isfield(index,obname))
index.(obname).name=obname;
index.(obname).long_name=long_obname;
index.(obname).rawdatacount=0;
index.(obname).datastartindex=zeros(NumOfSeg,1);
index.(obname).arrayDim=zeros(NumOfSeg,1);
index.(obname).nValues=zeros(NumOfSeg,1);
index.(obname).byteSize=zeros(NumOfSeg,1);
index.(obname).index=zeros(NumOfSeg,1);
index.(obname).rawdataoffset=zeros(NumOfSeg,1);
index.(obname).multiplier=ones(NumOfSeg,1);
index.(obname).skip=zeros(NumOfSeg,1);
end
if (kTocNewObjectList)
objOrderList{q}=obname;
else
if ~ismember(obname,objOrderList)
objOrderList{end+1}=obname;
end
end
%Get the raw data Index
rawdataindex=fread(fid,1,'uint32',kTocEndian);
if rawdataindex==0
if segCnt==0
e=errordlg(sprintf('Seqment %.0f within ''%s'' has ''rawdataindex'' value of 0 (%s.m).',segCnt,...
TDMSFileNameShort,mfilename),'Incorrect ''rawdataindex''');
uiwait(e)
end
if kTocRawData
if (kTocNewObjectList)
ccnt=index.(obname).rawdatacount+1;
else
ccnt=index.(obname).rawdatacount;
end
index.(obname).rawdatacount=ccnt;
index.(obname).datastartindex(ccnt)=SegInfo.DataStartPosn(segCnt);
index.(obname).arrayDim(ccnt)=index.(obname).arrayDim(ccnt-1);
index.(obname).nValues(ccnt)=index.(obname).nValues(ccnt-1);
index.(obname).byteSize(ccnt)=index.(obname).byteSize(ccnt-1);
index.(obname).index(ccnt)=segCnt;
SegInfo.NumChan(segCnt)=SegInfo.NumChan(segCnt)+1;
end
elseif rawdataindex+1==2^32
%Objects raw data index matches previous index - no changes. The root object will always have an
%'FFFFFFFF' entry
if strcmpi(index.(obname).long_name,'Root')
index.(obname).rawdataindex=0;
else
%Need to account for the case where an object (besides the 'root') is added that has no data but reports
%using previous.
if newob
index.(obname).rawdataindex=0;
else
if kTocRawData
if (kTocNewObjectList)
ccnt=index.(obname).rawdatacount+1;
else
ccnt=index.(obname).rawdatacount;
end
index.(obname).rawdatacount=ccnt;
index.(obname).datastartindex(ccnt)=SegInfo.DataStartPosn(segCnt);
index.(obname).arrayDim(ccnt)=index.(obname).arrayDim(ccnt-1);
index.(obname).nValues(ccnt)=index.(obname).nValues(ccnt-1);
index.(obname).byteSize(ccnt)=index.(obname).byteSize(ccnt-1);
index.(obname).index(ccnt)=segCnt;
SegInfo.NumChan(segCnt)=SegInfo.NumChan(segCnt)+1;
end
end
end
else
%Get new object information
if (kTocNewObjectList)
ccnt=index.(obname).rawdatacount+1;
else
ccnt=index.(obname).rawdatacount;
if (ccnt==0)
ccnt=1;
end
end
index.(obname).rawdatacount=ccnt;
index.(obname).datastartindex(ccnt)=SegInfo.DataStartPosn(segCnt);
%index(end).lenOfIndexInfo=fread(fid,1,'uint32');
index.(obname).dataType=fread(fid,1,'uint32',kTocEndian);
if (index.(obname).dataType~=32)
index.(obname).datasize=getDataSize(index.(obname).dataType);
end
index.(obname).arrayDim(ccnt)=fread(fid,1,'uint32',kTocEndian);
index.(obname).nValues(ccnt)=fread(fid,1,'uint64',kTocEndian);
index.(obname).index(ccnt)=segCnt;
SegInfo.NumChan(segCnt)=SegInfo.NumChan(segCnt)+1;
if index.(obname).dataType==32
%Datatype is a string
index.(obname).byteSize(ccnt)=fread(fid,1,'uint64',kTocEndian);
else
index.(obname).byteSize(ccnt)=0;
end
end
%Get the properties
numProps=fread(fid,1,'uint32',kTocEndian);
if numProps>0
if isfield(index.(obname),'PropertyInfo')
PropertyInfo=index.(obname).PropertyInfo;
else
clear PropertyInfo
end
for p=1:numProps
propNameLength=fread(fid,1,'uint32',kTocEndian);
switch 1
case 1
PropName=fread(fid,propNameLength,'*uint8','l')';
PropName=native2unicode(PropName,'UTF-8');
case 2
PropName=fread(fid,propNameLength,'uint8=>char','l')';
otherwise
end
propsDataType=fread(fid,1,'uint32',kTocEndian);
%Create property's name. Use a generic field name to avoid issues with strings that are too long and/or
%characters that cannot be used in MATLAB variable names. The actual property name is retained for the
%final output structure.
if exist('PropertyInfo','var')
%Check to see if the property already exists for this object. Need to get the existing 'PropertyInfo'
%structure for this object. The 'PropertyInfo' structure is not necessarily the same for every
%object in the data file.
PropIndex=find(strcmpi({PropertyInfo.Name},PropName));
if isempty(PropIndex)
%Is does not exist, so create the generic name field name
propExists=false;
PropIndex=numel(PropertyInfo)+1;
propsName=sprintf('Property%.0f',PropIndex);
PropertyInfo(PropIndex).Name=PropName;
PropertyInfo(PropIndex).FieldName=propsName;
else
%Assign the generic field name
propExists=true;
propsName=PropertyInfo(PropIndex).FieldName;
end
else
%No properties exist for this object, so create the first one using a generic name field name.
propExists=false;
PropIndex=p;
propsName=sprintf('Property%.0f',PropIndex);
PropertyInfo(PropIndex).Name=PropName;
PropertyInfo(PropIndex).FieldName=propsName;
end
dataExists=isfield(index.(obname),'data');
%Get the number of samples already found and in the object
if dataExists
nsamps=index.(obname).nsamples+1;
else
nsamps=0;
end
if propsDataType==32 %String data type
PropertyInfo(PropIndex).DataType='String';
propsValueLength=fread(fid,1,'uint32',kTocEndian);
propsValue=convertToText(fread(fid,propsValueLength,'uint8=>char',kTocEndian))';
if propExists
if isfield(index.(obname).(propsName),'cnt')
cnt=index.(obname).(propsName).cnt+1;
else
cnt=1;
end
index.(obname).(propsName).cnt=cnt;
index.(obname).(propsName).value{cnt}=propsValue;
index.(obname).(propsName).samples(cnt)=nsamps;
else
if strcmp(index.(obname).long_name,'Root')
%Header data
index.(obname).(propsName).name=index.(obname).long_name;
index.(obname).(propsName).value={propsValue};
index.(obname).(propsName).cnt=1;
else
index.(obname).(propsName).name=PropertyInfo(PropIndex).Name;
index.(obname).(propsName).datatype=PropertyInfo(PropIndex).DataType;
index.(obname).(propsName).cnt=1;
index.(obname).(propsName).value=cell(nsamps,1); %Pre-allocation
index.(obname).(propsName).samples=zeros(nsamps,1); %Pre-allocation
if iscell(propsValue)
index.(obname).(propsName).value(1)=propsValue;
else
index.(obname).(propsName).value(1)={propsValue};
end
index.(obname).(propsName).samples(1)=nsamps;
end
end
else %Numeric data type
if propsDataType==68 %Timestamp
PropertyInfo(PropIndex).DataType='Time';
%Timestamp data type
if kTocBigEndian
tsec=fread(fid,1,'uint64','b')+fread(fid,1,'uint64','b')/2^64; %time since Jan-1-1904 in seconds
else
tsec=fread(fid,1,'uint64','l')/2^64+fread(fid,1,'uint64','l'); %time since Jan-1-1904 in seconds
end
%tsec=fread(fid,1,'uint64',kTocEndian)/2^64+fread(fid,1,'uint64',kTocEndian); %time since Jan-1-1904 in seconds
%R. Seltzer: Not sure why '5/24' (5 hours) is subtracted from the time value. That's how it was
%coded in the original function I downloaded from MATLAB Central. But I found it to be 1 hour too
%much. So, I changed it to '4/24'.
%propsValue=tsec/86400+695422-5/24; %/864000 convert to days; +695422 days from Jan-0-0000 to Jan-1-1904
propsValue=tsec/86400+695422-4/24; %/864000 convert to days; +695422 days from Jan-0-0000 to Jan-1-1904
else %Numeric
PropertyInfo(PropIndex).DataType='Numeric';
matType=LV2MatlabDataType(propsDataType);
if strcmp(matType,'Undefined')
e=errordlg(sprintf('No MATLAB data type defined for a ''Property Data Type'' value of ''%.0f''.',...
propsDataType),'Undefined Property Data Type');
uiwait(e)
fclose(fid);
return
end
if strcmp(matType,'uint8=>char')
propsValue=convertToText(fread(fid,1,'uint8',kTocEndian));
else
propsValue=fread(fid,1,matType,kTocEndian);
end
end
if propExists
cnt=index.(obname).(propsName).cnt+1;
index.(obname).(propsName).cnt=cnt;
index.(obname).(propsName).value(cnt)=propsValue;
index.(obname).(propsName).samples(cnt)=nsamps;
else
index.(obname).(propsName).name=PropertyInfo(PropIndex).Name;
index.(obname).(propsName).datatype=PropertyInfo(PropIndex).DataType;
index.(obname).(propsName).cnt=1;
index.(obname).(propsName).value=NaN(nsamps,1); %Pre-allocation
index.(obname).(propsName).samples=zeros(nsamps,1); %Pre-allocation
index.(obname).(propsName).value(1)=propsValue;
index.(obname).(propsName).samples(1)=nsamps;
end
end
end %'end' for the 'Property' loop
index.(obname).PropertyInfo=PropertyInfo;
end
end %'end' for the 'Objects' loop
end
%Address Decimation and Interleaving
if (kTocRawData)
singleSegDataSize=0;
rawDataBytes=0;
chanBytes=0;
rawDataOffset=0;
for kk=1:numel(objOrderList)
obname=objOrderList{kk};
ccnt=hasRawDataInSeg(segCnt,index.(obname));
if ccnt>0 %If segement has raw data
index.(obname).rawdataoffset(ccnt)=rawDataOffset;
if index.(obname).dataType==32 %Datatype is a string
rawDataBytes=index.(obname).byteSize(ccnt);
else
rawDataBytes=index.(obname).nValues(ccnt)*index.(obname).datasize;
if kTocInterleavedData
chanBytes=index.(obname).datasize;
else
chanBytes=rawDataBytes;
end
end
rawDataOffset=rawDataOffset+chanBytes;
singleSegDataSize=singleSegDataSize+rawDataBytes;
end
end
%Calculate the offset for each segment. The offset is the amount
%of bytes for one non appened (non optimized) segement
% for kk=1:numel(objOrderList)
% obname=objOrderList{kk};
% if hasRawDataInSeg(segCnt,index.(obname)) %If segement has raw data
% index.(obname).rawdataoffset(ccnt)=rawdataoffset; %This will set the offset correctly for dec data.
%
% %Update the singleSegDataSize value for the next object
% if index.(obname).dataType==32 %Datatype is a string
% singleSegDataSize=singleSegDataSize+index.(obname).byteSize(ccnt);
% else
% singleSegDataSize=singleSegDataSize+index.(obname).nValues(ccnt)*index.(obname).datasize;
% end
%
% if kTocInterleavedData
% rawdataoffset=
% else
% rawdataoffset=singleSegDataSize;
% end
% end
% end
%
%As part of LV's optimization, it will append back to back
%segements into one segement (when nothing changes in meta data).
if (singleSegDataSize~=SegInfo.DataLength(segCnt)) %Multiple appended raw segement (offset did not grow to be greater than the DataLength)
numAppSegs=floor(SegInfo.DataLength(segCnt)/singleSegDataSize); %Number of appened (optimized segements)
for kk=1:numel(objOrderList) %Loop through all the objects
obname=objOrderList{kk};
ccnt=hasRawDataInSeg(segCnt,index.(obname));
if ccnt>0 %If segement has raw data
if kTocInterleavedData
if index.(obname).dataType==32 %Datatype is a string
error('Interleaved string channels are not supported.')
end
index.(obname).multiplier(ccnt)=numAppSegs*index.(obname).nValues(ccnt);
index.(obname).skip(ccnt)=singleSegDataSize/index.(obname).nValues(ccnt)-index.(obname).datasize;
index.(obname).nValues(ccnt)=1;
else %Decimated Data
index.(obname).multiplier(ccnt)=numAppSegs;
if index.(obname).dataType==32 %Datatype is a string
index.(obname).skip(ccnt)=singleSegDataSize-index.(obname).byteSize(ccnt);
else
index.(obname).skip(ccnt)=singleSegDataSize-index.(obname).nValues(ccnt)*index.(obname).datasize;
end
end
end
end
else %If single segement
if kTocInterleavedData
for kk=1:numel(objOrderList) %Loop through all the objects
obname=objOrderList{kk};
ccnt=hasRawDataInSeg(segCnt,index.(obname));
if ccnt>0 %If segement has raw data
if (index.(obname).index(ccnt)==segCnt) %If the object has rawdata in the current segment
if index.(obname).dataType==32 %Datatype is a string
error('Interleaved string channels are not supported.')
end
index.(obname).multiplier(ccnt)=index.(obname).nValues(ccnt);
index.(obname).skip(ccnt)=singleSegDataSize/index.(obname).nValues(ccnt)-index.(obname).datasize;
index.(obname).nValues(ccnt)=1;
end
end
end
end
end
% if (offset~=SegInfo.DataLength(segCnt)) %If the offset grew to be greater than the DataLength (from Meta Data)
% multiplier=floor(SegInfo.DataLength(segCnt)/offset);
% for kk=1:numel(objOrderList) %Loop through the objects
% obname=objOrderList{kk};
% ccnt=index.(obname).rawdatacount;
% if ccnt>0
% index.(obname).multiplier(segCnt)=multiplier;
% if index.(obname).dataType==32 %Datatype is a string
% index.(obname).skip(ccnt)=offset-index.(obname).byteSize(ccnt);
% else
% index.(obname).skip(ccnt)=offset-index.(obname).nValues(ccnt)*index.(obname).datasize;
% end
% end
% end
%
% end
aa=1;
% %Now adjust multiplier and skip if the data is interleaved
% if kTocInterleavedData
% if muliplier>0 %Muliple raw data segements appened
%
% else %Single raw segement
%
% end
% end
end
end
%clean up the index if it has to much data
fnm=fieldnames(index);
for kk=1:numel(fnm)
ccnt=index.(fnm{kk}).rawdatacount+1;
index.(fnm{kk}).datastartindex(ccnt:end)=[];
index.(fnm{kk}).arrayDim(ccnt:end)=[];
index.(fnm{kk}).nValues(ccnt:end)=[];
index.(fnm{kk}).byteSize(ccnt:end)=[];
index.(fnm{kk}).index(ccnt:end)=[];
index.(fnm{kk}).rawdataoffset(ccnt:end)=[];
index.(fnm{kk}).multiplier(ccnt:end)=[];
index.(fnm{kk}).skip(ccnt:end)=[];
end
end
%%
function ob=getData(fid,index,SegInfo)
%Using the file id (fid) and the index database, get the raw data.
%Returns the data in the ob structure. The fields in the structure are the
%generic object names:
% ob.Object3.data - raw data
% ob.Object3.nsamples - number of samples
%
% Note that not all of the objects in the index may not be in ob as it
% contains only those objects which have raw data.
ob=[];
fnm=fieldnames(index); %Get the object names
for kk=1:length(fnm) %Loop through objects
id=index.(fnm{kk});
nsamples=sum(id.nValues.*id.multiplier);
if id.rawdatacount>0 %Only work with channels with raw data
cname=id.name;
ob.(cname).nsamples=0;
%Initialize the data matrix
if id.dataType==32
ob.(cname).data=cell(nsamples,1);
else
ob.(cname).data=zeros(nsamples,1);
end
for rr=1:id.rawdatacount
%Loop through each of the segments and read the raw data
%Move to the raw data start position
fseek(fid,id.datastartindex(rr)+id.rawdataoffset(rr),'bof');
nvals=id.nValues(rr);
segmentNum=index.(cname).index(rr);
segInterleaved=SegInfo.SegInterleaved(segmentNum);
if SegInfo.SegBigEndian(segmentNum)
kTocEndian='b';
else
kTocEndian='l';
end
numChan=SegInfo.NumChan(segmentNum);
if nvals>0 %If there is data in this segement
switch id.dataType
case 32 %String
%From the National Instruments web page (http://zone.ni.com/devzone/cda/tut/p/id/5696) under the
%'Raw Data' description on page 4:
%String type channels are preprocessed for fast random access. All strings are concatenated to a
%contiguous piece of memory. The offset of the first character of each string in this contiguous piece
%of memory is stored to an array of unsigned 32-bit integers. This array of offset values is stored
%first, followed by the concatenated string values. This layout allows client applications to access
%any string value from anywhere in the file by repositioning the file pointer a maximum of three times
%and without reading any data that is not needed by the client.
data=cell(1,nvals*id.multiplier(rr)); %Pre-allocation
for mm=1:id.multiplier(rr)
StrOffsetArray=fread(fid,nvals,'uint32','l');
for dcnt=1:nvals
if dcnt==1
StrLength=StrOffsetArray(dcnt);
else
StrLength=StrOffsetArray(dcnt)-StrOffsetArray(dcnt-1);
end
data{1,dcnt+(mm-1)*nvals}=char(convertToText(fread(fid,StrLength,'uint8=>char','l'))');
end
if (id.multiplier(rr)>1)&&(id.skip(rr)>0)
fseek(fid,id.skip(rr),'cof');
end
end
cnt=nvals*id.multiplier(rr);
case 68 %Timestamp
%data=NaN(1,nvals); %Pre-allocation
data=NaN(1,nvals*id.multiplier(rr));
for mm=1:id.multiplier(rr)
dn=fread(fid,2*nvals,'uint64',kTocEndian);
tsec=dn(1:2:end)/2^64+dn(2:2:end);
data((mm-1)*nvals+1:(mm)*nvals)=tsec/86400+695422-4/24;
fseek(fid,id.skip(rr),'cof');
end
%{
for dcnt=1:nvals
tsec=fread(fid,1,'uint64')/2^64+fread(fid,1,'uint64'); %time since Jan-1-1904 in seconds
%R. Seltzer: Not sure why '5/24' (5 hours) is subtracted from the time value. That's how it was
%coded in the original function I downloaded from MATLAB Central. But I found it to be 1 hour too
%much. So, I changed it to '4/24'.
data(1,dcnt)=tsec/86400+695422-5/24; %/864000 convert to days; +695422 days from Jan-0-0000 to Jan-1-1904
data(1,dcnt)=tsec/86400+695422-4/24; %/864000 convert to days; +695422 days from Jan-0-0000 to Jan-1-1904
end
%}
cnt=nvals*id.multiplier(rr);
otherwise %Numeric
matType=LV2MatlabDataType(id.dataType);
if strcmp(matType,'Undefined') %Bad Data types catch
e=errordlg(sprintf('No MATLAB data type defined for a ''Raw Data Type'' value of ''%.0f''.',...
id.dataType),'Undefined Raw Data Type');
uiwait(e)
fclose(fid);
return
end
if (id.skip(rr)>0)
ntype=sprintf('%d*%s',nvals,matType);
[data,cnt]=fread(fid,nvals*id.multiplier(rr),ntype,id.skip(rr),kTocEndian);
if strcmp(matType,'uint8=>char')
data=convertToText(data);
end
else
% Added by Haench start
if (id.dataType == 524300) || (id.dataType == 1048589) % complex CDB data
[data,cnt]=fread(fid,2*nvals*id.multiplier(rr),matType,kTocEndian);
data= data(1:2:end)+1i*data(2:2:end);
cnt = cnt/2;
else
[data,cnt]=fread(fid,nvals*id.multiplier(rr),matType,kTocEndian);
end
% Haench end
% Original: [data,cnt]=fread(fid,nvals*id.multiplier(rr),matType,kTocEndian);
end
end
%Update the sample counter
if isfield(ob.(cname),'nsamples')
ssamples=ob.(cname).nsamples;
else
ssamples=0;
end
if (cnt>0)
ob.(cname).data(ssamples+1:ssamples+cnt,1)=data;
ob.(cname).nsamples=ssamples+cnt;
end
end
end
end
end
end
function [DataStructure,GroupNames]=postProcess(ob,index)
%Re-organize the 'ob' structure into a more user friendly format for output.
DataStructure.Root=[];
DataStructure.MeasuredData.Name=[];
DataStructure.MeasuredData.Data=[];
obFieldNames=fieldnames(index);
cntData=1;
for i=1:numel(obFieldNames)
cname=obFieldNames{i};
if strcmp(index.(cname).long_name,'Root')
DataStructure.Root.Name=index.(cname).long_name;
%Assign all the 'Property' values
if isfield(index.(cname),'PropertyInfo')
for p=1:numel(index.(cname).PropertyInfo)
cfield=index.(cname).PropertyInfo(p).FieldName;
if isfield(index.(cname).(cfield),'datatype')
DataType=index.(cname).(cfield).datatype;
else
%ASSUME a 'string' data type
DataType='String';
end
DataStructure.Root.Property(p).Name=index.(cname).PropertyInfo(p).Name;
switch DataType
case 'String'
if iscell(index.(cname).(cfield).value)
Value=index.(cname).(cfield).value';
else
Value=cellstr(index.(cname).(cfield).value);
end
case 'Time'
clear Value
if index.(cname).(cfield).cnt==1
if iscell(index.(cname).(cfield).value)
Value=datestr(cell2mat(index.(cname).(cfield).value),'dd-mmm-yyyy HH:MM:SS');
else
Value=datestr(index.(cname).(cfield).value,'dd-mmm-yyyy HH:MM:SS');
end
else
Value=cell(index.(cname).(cfield).cnt,1);
for c=1:index.(cname).(cfield).cnt
if iscell(index.(cname).(cfield).value)
Value(c)={datestr(cell2mat(index.(cname).(cfield).value),'dd-mmm-yyyy HH:MM:SS')};
else
Value(c)={datestr(index.(cname).(cfield).value,'dd-mmm-yyyy HH:MM:SS')};
end
end
end
case 'Numeric'
if isfield(index.(cname).(cfield),'cnt')
Value=NaN(index.(cname).(cfield).cnt,1);
else
if iscell(index.(cname).(cfield).value)
Value=NaN(numel(cell2mat(index.(cname).(cfield).value)),1);
else
Value=NaN(numel(index.(cname).(cfield).value),1);
end
end
for c=1:numel(Value)
if iscell(index.(cname).(cfield).value)
Value(c)=index.(cname).(cfield).value{c};
else
Value(c)=index.(cname).(cfield).value(c);
end
end
otherwise
e=errordlg(sprintf(['No format defined for Data Type ''%s'' in the private function ''postProcess'' '...
'within %s.m.'],index.(cname).(cfield).datatype,mfilename),'Undefined Property Data Type');
uiwait(e)
return
end
if isempty(Value)
DataStructure.Root.Property(p).Value=[];
else
DataStructure.Root.Property(p).Value=Value;
end
end
end
end
DataStructure.MeasuredData(cntData).Name=index.(cname).long_name;
%Should only need the 'ShortName' for debugging the function
%DataStructure.MeasuredData(cntData).ShortName=cname;
if (isfield(ob,cname))
if isfield(ob.(cname),'data')
DataStructure.MeasuredData(cntData).Data=ob.(cname).data;
%The following field is redundant because the information can be obtained from the size of the 'Data' field.
DataStructure.MeasuredData(cntData).Total_Samples=ob.(cname).nsamples;
else
DataStructure.MeasuredData(cntData).Data=[];
DataStructure.MeasuredData(cntData).Total_Samples=0;
end
else
DataStructure.MeasuredData(cntData).Data=[];
DataStructure.MeasuredData(cntData).Total_Samples=0;
end
%Assign all the 'Property' values
if isfield(index.(cname),'PropertyInfo')
for p=1:numel(index.(cname).PropertyInfo)
cfield=index.(cname).PropertyInfo(p).FieldName;
DataStructure.MeasuredData(cntData).Property(p).Name=index.(cname).(cfield).name;
if strcmpi(DataStructure.MeasuredData(cntData).Property(p).Name,'Root')
Value=index.(cname).(cfield).value;
else
switch index.(cname).(cfield).datatype
case 'String'
clear Value
if index.(cname).(cfield).cnt==1
if iscell(index.(cname).(cfield).value)
Value=char(index.(cname).(cfield).value);
else
Value=index.(cname).(cfield).value;
end
else
Value=cell(index.(cname).(cfield).cnt,1);
for c=1:index.(cname).(cfield).cnt
if iscell(index.(cname).(cfield).value)
Value(c)=index.(cname).(cfield).value;
else
Value(c)={index.(cname).(cfield).value};
end
end
end
case 'Time'
clear Value
if index.(cname).(cfield).cnt==1
if iscell(index.(cname).(cfield).value)
Value=datestr(cell2mat(index.(cname).(cfield).value),'dd-mmm-yyyy HH:MM:SS');
else
Value=datestr(index.(cname).(cfield).value,'dd-mmm-yyyy HH:MM:SS');
end
else
Value=cell(index.(cname).(cfield).cnt,1);
for c=1:index.(cname).(cfield).cnt
if iscell(index.(cname).(cfield).value)
Value(c)={datestr(cell2mat(index.(cname).(cfield).value),'dd-mmm-yyyy HH:MM:SS')};
else
Value(c)={datestr(index.(cname).(cfield).value,'dd-mmm-yyyy HH:MM:SS')};
end
end
end
case 'Numeric'
if isfield(index.(cname).(cfield),'cnt')
Value=NaN(index.(cname).(cfield).cnt,1);
else
if iscell(index.(cname).(cfield).value)
Value=NaN(numel(cell2mat(index.(cname).(cfield).value)),1);
else
Value=NaN(numel(index.(cname).(cfield).value),1);
end
end
for c=1:numel(Value)
if iscell(index.(cname).(cfield).value)
Value(c)=index.(cname).(cfield).value{c};
else
Value(c)=index.(cname).(cfield).value(c);
end
end
otherwise
e=errordlg(sprintf(['No format defined for Data Type ''%s'' in the private function ''postProcess'' '...
'within %s.m.'],index.(cname).(cfield).datatype,mfilename),'Undefined Property Data Type');
uiwait(e)
return
end
end
if isempty(Value)
DataStructure.MeasuredData(cntData).Property(p).Value=[];
else
DataStructure.MeasuredData(cntData).Property(p).Value=Value;
end
end
else
DataStructure.MeasuredData(cntData).Property=[];
end
cntData = cntData + 1;
end %'end' for the 'groups/channels' loop
%Extract the Group names
GroupIndices=false(numel(DataStructure.MeasuredData),1);
for d=1:numel(DataStructure.MeasuredData)
if ~strcmpi(DataStructure.MeasuredData(d).Name,'Root')
if (DataStructure.MeasuredData(d).Total_Samples==0)
fs=strfind(DataStructure.MeasuredData(d).Name,'/');
if (isempty(fs))
GroupIndices(d)=true;
end
end
end
end
if any(GroupIndices)
GroupNames=sort({DataStructure.MeasuredData(GroupIndices).Name})';
else
GroupNames=[];
end
end
function SegInfo=addChanCount(SegInfo,channelinfo)
%This function determines the number of channels in each segement by
%looping through the channels in channelinfo. It looks at their index,
%which conatins the segments numbers in which that channel has data.
%It then increments the NumChan field in the segment info database
%(SegInfo).
fnm=fieldnames(channelinfo); %Get the channels
for ccnt=1:numel(fnm) %Loop through the channels
channelName=fnm{ccnt};
ind=channelinfo.(channelName).index; %Get the index
for icnt=1:numel(ind)
segNum=ind(icnt);
SegInfo.NumChan(segNum)=SegInfo.NumChan(segNum)+1;
end
end
end
function sz=getDataSize(LVType)
%Get the number of bytes for each LV data type. See LV2MatlabDataType.
switch(LVType)
case 0
sz=0;
case {1,5,33}
sz=1;
case 68
sz=16;
case {8,10}
sz=8;
case {3,7,9}
sz=4;
case {2,6}
sz=2;
case 32
e=errordlg('Do not call the getDataSize function for strings. Their size is written in the data file','Error');
uiwait(e)
sz=NaN;
case 11
sz=10;
% Added by Haench for tdsTypeComplexSingleFloat=0x08000c,tdsTypeComplexDoubleFloat=0x10000d,
case 524300
sz=8;
case 1048589
sz=16;
% end add haench
otherwise
error('LVData type %d is not defined',LVType)
end
end
function matType=LV2MatlabDataType(LVType)
%Cross Refernce Labview TDMS Data type to MATLAB
switch LVType
case 0 %tdsTypeVoid
matType='';
case 1 %tdsTypeI8
matType='int8';
case 2 %tdsTypeI16
matType='int16';
case 3 %tdsTypeI32
matType='int32';
case 4 %tdsTypeI64
matType='int64';
case 5 %tdsTypeU8
matType='uint8';
case 6 %tdsTypeU16
matType='uint16';
case 7 %tdsTypeU32
matType='uint32';
case 8 %tdsTypeU64
matType='uint64';
case 9 %tdsTypeSingleFloat
matType='single';
case 10 %tdsTypeDoubleFloat
matType='double';
case 11 %tdsTypeExtendedFloat
matType='10*char';
case 25 %tdsTypeSingleFloat with units
matType='Undefined';
case 26 %tdsTypeDoubleFloat with units
matType='Undefined';
case 27 %tdsTypeextendedFloat with units
matType='Undefined';
case 32 %tdsTypeString
matType='uint8=>char';
case 33 %tdsTypeBoolean
matType='bit1';
case 68 %tdsTypeTimeStamp
matType='2*int64';
% Added by Haench for tdsTypeComplexSingleFloat=0x08000c,tdsTypeComplexDoubleFloat=0x10000d,
case 524300
matType='single';
case 1048589
matType='double';
% end add haench
otherwise
matType='Undefined';
end
end
function ccnt=hasRawDataInSeg(segCnt,ob)
%Function to check and see if a channel has raw data in a particulr
%segment. If it does, return the index of index (ccnt). If not present,
%return ccn=0
indexList=ob.index;
ccnt=find(indexList==segCnt);
if isempty(ccnt)
ccnt=0;
else
ccnt=ccnt(1);
end
end
function text=convertToText(bytes)
%Convert numeric bytes to the character encoding localy set in MATLAB (TDMS uses UTF-8)
text=native2unicode(bytes,'UTF-8');
end
e文件
classdef NicoletFile < handle
% NICOLETFILE Reading Nicolet .e files.
%
% GETTING DATA
% You can load data from the .e file using the GETDATA method. The
% inputs to the method are the object, the segment of the data file
% that you want to load data from, the min, and max index you want to
% retrieve, and a vector of channels that you want to retrieve.
%
% Example:
% OUT = GETDATA(OBJ, 1, [1 1000], 1:10) will return the first 1000
% values on the first 10 channels of the first segment of the file.
%
% GETTING Nr OF SAMPLES
% Use the GETNRSAMPLES method to find the number of samples per
% channel in each data segment.
%
% WARNING!
% The .e format allows for changes in the TimeSeries map during the
% recording. This results in multiple TSINFO structures. Depending on
% where these structures are located in the .e file, the appropriate
% TSINFO structure should be used. However, there seems to be a bug in
% the Nicolet .e file writer which sometimes renders the TSINFO structures
% unreadable on disk (verify with hex-edit).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2013 Trustees of the University of Pennsylvania
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Joost Wagenaar, Jan 2015
% Cristian Donos, Dec 2015
% Jan Brogger, Jun 2016
% Callum Stewart, Apr 2018
properties
fileName
patientInfo
segments
eventMarkers
end
properties (Hidden)
sections
index
sigInfo
tsInfos
chInfo
notchFreq
montage
Qi
Qii
allIndexIDs
useTSinfoIdx = 1
end
properties (Constant, Hidden)
LABELSIZE = 32
TSLABELSIZE = 64
UNITSIZE = 16
ITEMNAMESIZE = 64
end
methods
function obj = NicoletFile(filename)
h = fopen(filename,'r','ieee-le');
[folder, name, ext] = fileparts(filename);
assert(strcmp(ext,'.e'), 'File extention must be .e');
if isempty(folder)
filename = fullfile(pwd,filename);
end
obj.fileName = filename;
% Get init
misc1 = fread(h,5, 'uint32'); %#ok<NASGU>
unknown = fread(h,1,'uint32'); %#ok<NASGU>
indexIdx = fread(h,1,'uint32');
% Get TAGS structure and Channel IDS
fseek(h, 172,'bof');
nrTags = fread(h,1, 'uint32');
Tags = struct();
for i = 1:nrTags
Tags(i).tag = deblank(cast(fread(h, 40, 'uint16'),'char')');
Tags(i).index = fread(h,1,'uint32');
switch Tags(i).tag
case 'ExtraDataTags'
Tags(i).IDStr = 'ExtraDataTags';
case 'SegmentStream'
Tags(i).IDStr = 'SegmentStream';
case 'DataStream'
Tags(i).IDStr = 'DataStream';
case 'InfoChangeStream'
Tags(i).IDStr = 'InfoChangeStream';
case 'InfoGuids'
Tags(i).IDStr = 'InfoGuids';
case '{A271CCCB-515D-4590-B6A1-DC170C8D6EE2}'
Tags(i).IDStr = 'TSGUID';
case '{8A19AA48-BEA0-40D5-B89F-667FC578D635}'
Tags(i).IDStr = 'DERIVATIONGUID';
case '{F824D60C-995E-4D94-9578-893C755ECB99}'
Tags(i).IDStr = 'FILTERGUID';
case '{02950361-35BB-4A22-9F0B-C78AAA5DB094}'
Tags(i).IDStr = 'DISPLAYGUID';
case '{8E94EF21-70F5-11D3-8F72-00105A9AFD56}'
Tags(i).IDStr = 'FILEINFOGUID';
case '{E4138BC0-7733-11D3-8685-0050044DAAB1}'
Tags(i).IDStr = 'SRINFOGUID';
case '{C728E565-E5A0-4419-93D2-F6CFC69F3B8F}'
Tags(i).IDStr = 'EVENTTYPEINFOGUID';
case '{D01B34A0-9DBD-11D3-93D3-00500400C148}'
Tags(i).IDStr = 'AUDIOINFOGUID';
case '{BF7C95EF-6C3B-4E70-9E11-779BFFF58EA7}'
Tags(i).IDStr = 'CHANNELGUID';
case '{2DEB82A1-D15F-4770-A4A4-CF03815F52DE}'
Tags(i).IDStr = 'INPUTGUID';
case '{5B036022-2EDC-465F-86EC-C0A4AB1A7A91}'
Tags(i).IDStr = 'INPUTSETTINGSGUID';
case '{99A636F2-51F7-4B9D-9569-C7D45058431A}'
Tags(i).IDStr = 'PHOTICGUID';
case '{55C5E044-5541-4594-9E35-5B3004EF7647}'
Tags(i).IDStr = 'ERRORGUID';
case '{223A3CA0-B5AC-43FB-B0A8-74CF8752BDBE}'
Tags(i).IDStr = 'VIDEOGUID';
case '{0623B545-38BE-4939-B9D0-55F5E241278D}'
Tags(i).IDStr = 'DETECTIONPARAMSGUID';
case '{CE06297D-D9D6-4E4B-8EAC-305EA1243EAB}'
Tags(i).IDStr = 'PAGEGUID';
case '{782B34E8-8E51-4BB9-9701-3227BB882A23}'
Tags(i).IDStr = 'ACCINFOGUID';
case '{3A6E8546-D144-4B55-A2C7-40DF579ED11E}'
Tags(i).IDStr = 'RECCTRLGUID';
case '{D046F2B0-5130-41B1-ABD7-38C12B32FAC3}'
Tags(i).IDStr = 'GUID TRENDINFOGUID';
case '{CBEBA8E6-1CDA-4509-B6C2-6AC2EA7DB8F8}'
Tags(i).IDStr = 'HWINFOGUID';
case '{E11C4CBA-0753-4655-A1E9-2B2309D1545B}'
Tags(i).IDStr = 'VIDEOSYNCGUID';
case '{B9344241-7AC1-42B5-BE9B-B7AFA16CBFA5}'
Tags(i).IDStr = 'SLEEPSCOREINFOGUID';
case '{15B41C32-0294-440E-ADFF-DD8B61C8B5AE}'
Tags(i).IDStr = 'FOURIERSETTINGSGUID';
case '{024FA81F-6A83-43C8-8C82-241A5501F0A1}'
Tags(i).IDStr = 'SPECTRUMGUID';
case '{8032E68A-EA3E-42E8-893E-6E93C59ED515}'
Tags(i).IDStr = 'SIGNALINFOGUID';
case '{30950D98-C39C-4352-AF3E-CB17D5B93DED}'
Tags(i).IDStr = 'SENSORINFOGUID';
case '{F5D39CD3-A340-4172-A1A3-78B2CDBCCB9F}'
Tags(i).IDStr = 'DERIVEDSIGNALINFOGUID';
case '{969FBB89-EE8E-4501-AD40-FB5A448BC4F9}'
Tags(i).IDStr = 'ARTIFACTINFOGUID';
case '{02948284-17EC-4538-A7FA-8E18BD65E167}'
Tags(i).IDStr = 'STUDYINFOGUID';
case '{D0B3FD0B-49D9-4BF0-8929-296DE5A55910}'
Tags(i).IDStr = 'PATIENTINFOGUID';
case '{7842FEF5-A686-459D-8196-769FC0AD99B3}'
Tags(i).IDStr = 'DOCUMENTINFOGUID';
case '{BCDAEE87-2496-4DF4-B07C-8B4E31E3C495}'
Tags(i).IDStr = 'USERSINFOGUID';
case '{B799F680-72A4-11D3-93D3-00500400C148}'
Tags(i).IDStr = 'EVENTGUID';
case '{AF2B3281-7FCE-11D2-B2DE-00104B6FC652}'
Tags(i).IDStr = 'SHORTSAMPLESGUID';
case '{89A091B3-972E-4DA2-9266-261B186302A9}'
Tags(i).IDStr = 'DELAYLINESAMPLESGUID';
case '{291E2381-B3B4-44D1-BB77-8CF5C24420D7}'
Tags(i).IDStr = 'GENERALSAMPLESGUID';
case '{5F11C628-FCCC-4FDD-B429-5EC94CB3AFEB}'
Tags(i).IDStr = 'FILTERSAMPLESGUID';
case '{728087F8-73E1-44D1-8882-C770976478A2}'
Tags(i).IDStr = 'DATEXDATAGUID';
case '{35F356D9-0F1C-4DFE-8286-D3DB3346FD75}'
Tags(i).IDStr = 'TESTINFOGUID';
otherwise
if isstrprop(Tags(i).tag, 'digit')
Tags(i).IDStr = num2str(Tags(i).tag);
else
Tags(i).IDStr = 'UNKNOWN';
end
end
end
obj.sections = Tags;
%% QI index
fseek(h, 172208,'bof');
obj.Qi=struct();
obj.Qi.nrEntries = fread(h,1,'uint32');
obj.Qi.misc1 = fread(h,1,'uint32');
obj.Qi.indexIdx = fread(h,1,'uint32');
obj.Qi.misc3 = fread(h,1,'uint32');
obj.Qi.LQi = fread(h,1,'uint64')';
obj.Qi.firstIdx = fread(h,nrTags,'uint64');
% Don't know what this index is for... Not required to get data and
% can be huge...
% fseek(h, 188664,'bof');
% Qindex = struct();
% for i = 1:obj.Qi.LQi
% % Qindex(i).ftel = ftell(h);
% Qindex(i).index = fread(h,2,'uint16')'; %4
% Qindex(i).misc1 = fread(h,1,'uint32'); %8
% Qindex(i).indexIdx = fread(h,1,'uint32'); %12
% Qindex(i).misc2 = fread(h,3,'uint32')'; %24
% Qindex(i).sectionIdx = fread(h,1,'uint32');%28
% Qindex(i).misc3 = fread(h,1,'uint32'); %32
% Qindex(i).offset = fread(h,1,'uint64'); % 40
% Qindex(i).blockL = fread(h,1,'uint32');%44
% Qindex(i).dataL = fread(h,1,'uint32')';%48
% end
% obj.Qi.index = Qindex;
%% Get Main Index:
% Index consists of multiple blocks, after each block is the pointer
% to the next block. Total number of entries is in obj.Qi.nrEntries
Index = struct();
curIdx = 0;
nextIndexPointer = indexIdx;
fprintf('Parsing index ');
curIdx2 = 1;
while curIdx < obj.Qi.nrEntries
if mod(curIdx2,20)
fprintf('.');
else
fprintf('\n.');
end
fseek(h, nextIndexPointer, 'bof');
nrIdx = fread(h,1, 'uint64');
Index(curIdx + nrIdx).sectionIdx = 0; % Preallocate next set of indices
var = fread(h,3*nrIdx, 'uint64');
for i = 1: nrIdx
% Index(curIdx + i).sectionIdx = fread(h,1, 'uint64');
% Index(curIdx + i).offset = fread(h,1, 'uint64');
% Index(curIdx + i).blockL = fread(h,1, 'uint32');
% Index(curIdx + i).sectionL = fread(h,1, 'uint32');
Index(curIdx + i).sectionIdx = var(3*(i-1)+1);
Index(curIdx + i).offset = var(3*(i-1)+2);
Index(curIdx + i).blockL = mod(var(3*(i-1)+3),2^32);
Index(curIdx + i).sectionL = round(var(3*(i-1)+3)/2^32);
end
nextIndexPointer = fread(h,1, 'uint64');
curIdx = curIdx + i;
curIdx2=curIdx2+1;
end
fprintf('done\n');
obj.index = Index;
obj.allIndexIDs = [obj.index.sectionIdx];
%---READ DYNAMIC PACKETS---%
dynamicPackets = struct();
indexIdx = Tags(find(strcmp({Tags.IDStr},'InfoChangeStream'),1)).index;
offset = Index(indexIdx).offset;
nrDynamicPackets = Index(indexIdx).sectionL / 48;
fseek(h, offset, 'bof');
%Read first only the dynamic packets structure without actual data
for i = 1: nrDynamicPackets
dynamicPackets(i).offset = offset+i*48;
guidmixed = fread(h,16, 'uint8')';
guidnonmixed = [guidmixed(04), guidmixed(03), guidmixed(02), guidmixed(01), ...
guidmixed(06), guidmixed(05), guidmixed(08), guidmixed(07), ...
guidmixed(09), guidmixed(10), guidmixed(11), guidmixed(12), ...
guidmixed(13), guidmixed(14), guidmixed(15), guidmixed(16)];
dynamicPackets(i).guid = num2str(guidnonmixed, '%02X');
dynamicPackets(i).guidAsStr = sprintf('{%02X%02X%02X%02X-%02X%02X-%02X%02X-%02X%02X-%02X%02X%02X%02X%02X%02X}', guidnonmixed);
dynamicPackets(i).date = datenum(1899,12,31) + fread(h,1,'double');
dynamicPackets(i).datefrac = fread(h,1,'double');
dynamicPackets(i).internalOffsetStart = fread(h,1, 'uint64')';
dynamicPackets(i).packetSize = fread(h,1, 'uint64')';
dynamicPackets(i).data = zeros(0, 1,'uint8');
switch dynamicPackets(i).guid
case 'BF7C95EF6C3B4E709E11779BFFF58EA7'
dynamicPackets(i).IDStr = 'CHANNELGUID';
case '8A19AA48BEA040D5B89F667FC578D635'
dynamicPackets(i).IDStr = 'DERIVATIONGUID';
case 'F824D60C995E4D949578893C755ECB99'
dynamicPackets(i).IDStr = 'FILTERGUID';
case '0295036135BB4A229F0BC78AAA5DB094'
dynamicPackets(i).IDStr = 'DISPLAYGUID';
case '782B34E88E514BB997013227BB882A23'
dynamicPackets(i).IDStr = 'ACCINFOGUID';
case 'A271CCCB515D4590B6A1DC170C8D6EE2'
dynamicPackets(i).IDStr = 'TSGUID';
case 'D01B34A09DBD11D393D300500400C148'
dynamicPackets(i).IDStr = 'AUDIOINFOGUID';
otherwise
dynamicPackets(i).IDStr = 'UNKNOWN';
end
end
%Then read the actual data from the pointers above
for i = 1: nrDynamicPackets
%Look up the GUID of this dynamic packet in the Tags
% to find the section index
infoIdx = Tags(find(strcmp({Tags.tag},dynamicPackets(i).guidAsStr),1)).index;
%Matching index segments
indexInstances = Index([Index.sectionIdx] == infoIdx);
%Then, treat all these sections as one contiguous memory block
% and grab this packet across these instances
internalOffset = 0;
remainingDataToRead = dynamicPackets(i).packetSize;
%disp(['Target packet ' dynamicPackets(i).IDStr ' : ' num2str(dynamicPackets(i).internalOffsetStart) ' to ' num2str(dynamicPackets(i).internalOffsetStart+dynamicPackets(i).packetSize) ' target read length ' num2str(remainingDataToRead)]);
currentTargetStart = dynamicPackets(i).internalOffsetStart;
for j = 1: size(indexInstances,2)
currentInstance = indexInstances(j);
%hitInThisSegment = '';
if (internalOffset <= currentTargetStart) && (internalOffset+currentInstance.sectionL) >= currentTargetStart
startAt = currentTargetStart;
stopAt = min(startAt+remainingDataToRead, internalOffset+currentInstance.sectionL);
readLength = stopAt-startAt;
filePosStart = currentInstance.offset+startAt-internalOffset;
fseek(h,filePosStart, 'bof');
dataPart = fread(h,readLength,'uint8=>uint8');
dynamicPackets(i).data = cat(1, dynamicPackets(i).data, dataPart);
%hitInThisSegment = ['HIT at ' num2str(startAt) ' to ' num2str(stopAt)];
%if (readLength < remainingDataToRead)
% hitInThisSegment = [hitInThisSegment ' (partial ' num2str(readLength) ' )'];
%else
% hitInThisSegment = [hitInThisSegment ' (finished - this segment contributed ' num2str(readLength) ' )'];
%end
%hitInThisSegment = [hitInThisSegment ' abs file pos ' num2str(filePosStart) ' - ' num2str(filePosStart+readLength)];
remainingDataToRead = remainingDataToRead-readLength;
currentTargetStart = currentTargetStart + readLength;
end
%disp([' Index ' num2str(j) ' Offset: ' num2str(internalOffset) ' to ' num2str(internalOffset+currentInstance.sectionL) ' ' num2str(hitInThisSegment)]);
internalOffset = internalOffset + currentInstance.sectionL;
end
end
%% Get PatientGUID
info = struct();
infoProps = { 'patientID', 'firstName','middleName','lastName',...
'altID','mothersMaidenName','DOB','DOD','street','sexID','phone',...
'notes','dominance','siteID','suffix','prefix','degree','apartment',...
'city','state','country','language','height','weight','race','religion',...
'maritalStatus'};
infoIdxStruct = Tags(find(strcmp({Tags.IDStr},'PATIENTINFOGUID'),1));
if ~isempty(infoIdxStruct)
infoIdx = infoIdxStruct.index;
indexInstance = Index(find([Index.sectionIdx]==infoIdx,1));
fseek(h, indexInstance.offset,'bof');
guid = fread(h, 16, 'uint8'); %#ok<NASGU>
lSection = fread(h, 1, 'uint64'); %#ok<NASGU>
% reserved = fread(h, 3, 'uint16'); %#ok<NASGU>
nrValues = fread(h,1,'uint64');
nrBstr = fread(h,1,'uint64');
for i = 1:nrValues
id = fread(h,1,'uint64');
switch id
case {7,8}
unix_time = (fread(h,1, 'double')*(3600*24)) - 2209161600;% 2208988800; %8
obj.segments(i).dateStr = datestr(unix_time/86400 + datenum(1970,1,1));
value = datevec( obj.segments(i).dateStr );
value = value([3 2 1]);
case {23,24}
value = fread(h,1,'double');
otherwise
value = 0;
end
info.(infoProps{id}) = value;
end
strSetup = fread(h,nrBstr*2,'uint64');
for i=1:2:(nrBstr*2)
id = strSetup(i);
value = deblank(cast(fread(h, strSetup(i+1) + 1, 'uint16'),'char')');
info.(infoProps{id}) = value;
end
end
obj.patientInfo = info;
%% Get INFOGUID
infoIdx = Tags(find(strcmp({Tags.IDStr},'InfoGuids'),1)).index;
indexInstance = Index(find([Index.sectionIdx]==infoIdx,1));
fseek(h, indexInstance.offset,'bof');
% Ignoring, is list of GUIDS in file.
%% Get SignalInfo (SIGNALINFOGUID): One per file
SIG_struct = struct();
sensorIdx = Tags(find(strcmp({Tags.IDStr},'SIGNALINFOGUID'),1)).index;
indexInstance = Index(find([Index.sectionIdx]==sensorIdx,1));
fseek(h, indexInstance.offset,'bof');
SIG_struct.guid = fread(h, 16, 'uint8');
SIG_struct.name = fread(h, obj.ITEMNAMESIZE, '*char');
unkown = fread(h, 152, '*char'); %#ok<NASGU>
fseek(h, 512, 'cof');
nrIdx = fread(h,1, 'uint16'); %783
misc1 = fread(h,3, 'uint16'); %#ok<NASGU>
obj.sigInfo = struct();
for i = 1: nrIdx
obj.sigInfo(i).sensorName = deblank(cast(fread(h, obj.LABELSIZE, 'uint16'),'char')');
obj.sigInfo(i).transducer = deblank(cast(fread(h, obj.UNITSIZE, 'uint16'),'char')');
obj.sigInfo(i).guid = fread(h, 16, '*uint8');
obj.sigInfo(i).bBiPolar = logical(fread(h, 1 ,'uint32'));
obj.sigInfo(i).bAC = logical(fread(h, 1 ,'uint32'));
obj.sigInfo(i).bHighFilter = logical(fread(h, 1 ,'uint32'));
obj.sigInfo(i).color = fread(h, 1 ,'uint32');
reserved = fread(h, 256, '*char'); %#ok<NASGU>
end
%% Get CHANNELINFO (CHANNELGUID)
CH_struct = struct();
sensorIdx = Tags(find(strcmp({Tags.IDStr},'CHANNELGUID'),1)).index;
indexInstance = Index(find([Index.sectionIdx]==sensorIdx,1));
fseek(h, indexInstance.offset,'bof');
CH_struct.guid = fread(h, 16, 'uint8');
CH_struct.name = fread(h, obj.ITEMNAMESIZE, '*char');
fseek(h, 152, 'cof');
CH_struct.reserved = fread(h, 16, 'uint8');
CH_struct.deviceID = fread(h, 16, 'uint8');
fseek(h, 488, 'cof');
nrIdx = fread(h,2, 'uint32'); %783
obj.chInfo = struct();
for i = 1: nrIdx(2)
obj.chInfo(i).sensor = deblank(cast(fread(h, obj.LABELSIZE, 'uint16'),'char')');
obj.chInfo(i).samplingRate = fread(h,1,'double');
obj.chInfo(i).bOn = logical(fread(h, 1 ,'uint32'));
obj.chInfo(i).lInputID = fread(h, 1 ,'uint32');
obj.chInfo(i).lInputSettingID = fread(h,1,'uint32');
obj.chInfo(i).reserved = fread(h,4,'char');
fseek(h, 128, 'cof');
end
curIdx = 0;
for i = 1: length(obj.chInfo)
if obj.chInfo(i).bOn
obj.chInfo(i).indexID = curIdx;
curIdx = curIdx+1;
else
obj.chInfo(i).indexID = -1;
end
end
%% Get TS info (TSGUID):(One per segment, last used if no new for segment)
%% To simplify things, we only read the first TSINFO.
tsPackets = dynamicPackets(strcmp({dynamicPackets.IDStr},'TSGUID'));
if isempty(tsPackets)
warning(['No TSINFO found']);
else
obj.tsInfos = {};
for j = 1:length(tsPackets)
tsPacket = tsPackets(j);
elems = typecast(tsPacket.data(753:756),'uint32');
alloc = typecast(tsPacket.data(757:760),'uint32');
offset = 761;
tsInfo = struct();
for i = 1:elems
internalOffset = 0;
tsInfo(i).label = deblank(char(typecast(tsPacket.data(offset:(offset+obj.TSLABELSIZE-1))','uint16')));
internalOffset = internalOffset + obj.TSLABELSIZE*2;
tsInfo(i).activeSensor = deblank(char(typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+obj.LABELSIZE))','uint16')));
internalOffset = internalOffset + obj.TSLABELSIZE;
tsInfo(i).refSensor = deblank(char(typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','uint16')));
internalOffset = internalOffset + 8;
internalOffset = internalOffset + 56;
tsInfo(i).dLowCut = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
internalOffset = internalOffset + 8;
tsInfo(i).dHighCut = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
internalOffset = internalOffset + 8;
tsInfo(i).dSamplingRate = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
internalOffset = internalOffset + 8;
tsInfo(i).dResolution = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
internalOffset = internalOffset + 8;
tsInfo(i).bMark = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+2))','uint16');
internalOffset = internalOffset + 2;
tsInfo(i).bNotch = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+2))','uint16');
internalOffset = internalOffset + 2;
tsInfo(i).dEegOffset = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
offset = offset + 552;
%disp([num2str(i) ' : ' TSInfo(i).label ' : ' TSInfo(i).activeSensor ' : ' TSInfo(i).refSensor ' : ' num2str(TSInfo(i).samplingRate)]);
end
obj.tsInfos{j} = tsInfo;
end
end
% -- -- --
%% Get Segment Start Times
segmentIdx = Tags(find(strcmp({Tags.IDStr}, 'SegmentStream'),1)).index;
indexIdx = find([Index.sectionIdx] == segmentIdx, 1);
segmentInstance = Index(indexIdx);
nrSegments = segmentInstance.sectionL/152;
fseek(h, segmentInstance.offset,'bof');
obj.segments = struct();
for i = 1: nrSegments
dateOLE = fread(h,1, 'double');
obj.segments(i).dateOLE = dateOLE;
unix_time = (dateOLE*(3600*24)) - 2209161600;% 2208988800; %8
obj.segments(i).dateStr = datestr(unix_time/86400 + datenum(1970,1,1));
datev = datevec( obj.segments(i).dateStr );
obj.segments(i).startDate = datev(1:3);
obj.segments(i).startTime = datev(4:6);
fseek(h, 8 , 'cof'); %16
obj.segments(i).duration = fread(h,1, 'double');%24
fseek(h, 128 , 'cof'); %152
end
% Get nrValues per segment and channel
for iSeg = 1:length(obj.segments)
if iSeg > length(obj.tsInfos)
obj.tsInfos{iSeg} = obj.tsInfos{iSeg-1}
end
% Add Channel Names to segments
obj.segments(iSeg).chName = {obj.tsInfos{iSeg}.label};
obj.segments(iSeg).refName = {obj.tsInfos{iSeg}.refSensor};
obj.segments(iSeg).samplingRate = [obj.tsInfos{iSeg}.dSamplingRate];
obj.segments(iSeg).scale = [obj.tsInfos{iSeg}.dResolution];
end
%% Get events - Andrei Barborica, Dec 2015
% Find sequence of events, that are stored in the section tagged 'Events'
idxSection = find(strcmp('Events',{Tags.tag}));
indexIdx = find([obj.index.sectionIdx] == obj.sections(idxSection).index);
offset = obj.index(indexIdx).offset;
ePktLen = 272; % Event packet length, see EVENTPACKET definition
eMrkLen = 240; % Event marker length, see EVENTMARKER definition
evtPktGUID = hex2dec({'80', 'F6', '99', 'B7', 'A4', '72', 'D3', '11', '93', 'D3', '00', '50', '04', '00', 'C1', '48'}); % GUID for event packet header
HCEVENT_ANNOTATION = '{A5A95612-A7F8-11CF-831A-0800091B5BDA}';
HCEVENT_SEIZURE = '{A5A95646-A7F8-11CF-831A-0800091B5BDA}';
HCEVENT_FORMATCHANGE = '{08784382-C765-11D3-90CE-00104B6F4F70}';
HCEVENT_PHOTIC = '{6FF394DA-D1B8-46DA-B78F-866C67CF02AF}';
HCEVENT_POSTHYPERVENT = '{481DFC97-013C-4BC5-A203-871B0375A519}';
HCEVENT_REVIEWPROGRESS = '{725798BF-CD1C-4909-B793-6C7864C27AB7}';
HCEVENT_EXAMSTART = '{96315D79-5C24-4A65-B334-E31A95088D55}';
HCEVENT_HYPERVENTILATION = '{A5A95608-A7F8-11CF-831A-0800091B5BDA}';
HCEVENT_IMPEDANCE = '{A5A95617-A7F8-11CF-831A-0800091B5BDA}';
HCEVENT_AMPLIFIERDISCONNECT = '{A71A6DB5-4150-48BF-B462-1C40521EBD6F}';
HCEVENT_AMPLIFIERRECONNECT = '{6387C7C8-6F98-4886-9AF4-FA750ED300DE}';
HCEVENT_PAUSED = '{71EECE80-EBC4-41C7-BF26-E56911426FB4}';
DAYSECS = 86400.0; % From nrvdate.h
fseek(h,offset,'bof');
pktGUID = fread(h,16,'uint8');
pktLen = fread(h,1,'uint64');
obj.eventMarkers = struct();
i = 0; % Event counter
while (pktGUID == evtPktGUID)
i = i + 1;
% Please refer to EVENTMARKER structure in the Nervus file documentation
fseek(h,8,'cof'); % Skip eventID, not used
evtDate = fread(h,1,'double');
evtDateFraction = fread(h,1,'double');
obj.eventMarkers(i).dateOLE = evtDate;
obj.eventMarkers(i).dateFraction = evtDateFraction;
evtPOSIXTime = evtDate*DAYSECS + evtDateFraction - 2209161600;% 2208988800; %8
obj.eventMarkers(i).dateStr = datestr(evtPOSIXTime/DAYSECS + datenum(1970,1,1),'dd-mmmm-yyyy HH:MM:SS.FFF'); % Save fractions of seconds, as well
obj.eventMarkers(i).duration = fread(h,1,'double');
fseek(h,48,'cof');
evtUser = fread(h,12,'uint16');
obj.eventMarkers(i).user = deblank(char(evtUser).');
evtTextLen = fread(h,1,'uint64');
evtGUID = fread(h,16,'uint8');
obj.eventMarkers(i).GUID = sprintf('{%.2X%.2X%.2X%.2X-%.2X%.2X-%.2X%.2X-%.2X%.2X-%.2X%.2X%.2X%.2X%.2X%.2X}',evtGUID([4 3 2 1 6 5 8 7 9:16]));
fseek(h,16,'cof'); % Skip Reserved4 array
evtLabel = fread(h,32,'uint16'); % LABELSIZE = 32;
evtLabel = deblank(char(evtLabel).'); % Not used
eventMarkers(i).label = evtLabel;
%disp(sprintf('Offset: %s, TypeGUID:%s, User:%s, Label:%s',dec2hex(offset),evtGUID,evtUser,evtLabel));
% Only a subset of all event types are dealt with
switch obj.eventMarkers(i).GUID
case HCEVENT_SEIZURE
obj.eventMarkers(i).IDStr = 'Seizure';
%disp(' Seizure event');
case HCEVENT_ANNOTATION
obj.eventMarkers(i).IDStr = 'Annotation';
fseek(h,32,'cof'); % Skip Reserved5 array
evtAnnotation = fread(h,evtTextLen,'uint16');
obj.eventMarkers(i).annotation = deblank(char(evtAnnotation).');
%disp(sprintf(' Annotation:%s',evtAnnotation));
case HCEVENT_FORMATCHANGE
obj.eventMarkers(i).IDStr = 'Format change';
case HCEVENT_PHOTIC
obj.eventMarkers(i).IDStr = 'Photic';
case HCEVENT_POSTHYPERVENT
obj.eventMarkers(i).IDStr = 'Posthyperventilation';
case HCEVENT_REVIEWPROGRESS
obj.eventMarkers(i).IDStr = 'Review progress';
case HCEVENT_EXAMSTART
obj.eventMarkers(i).IDStr = 'Exam start';
case HCEVENT_HYPERVENTILATION
obj.eventMarkers(i).IDStr = 'Hyperventilation';
case HCEVENT_IMPEDANCE
obj.eventMarkers(i).IDStr = 'Impedance';
case HCEVENT_AMPLIFIERDISCONNECT
obj.eventMarkers(i).IDStr = 'Amplifier Disconnect';
case HCEVENT_AMPLIFIERRECONNECT
obj.eventMarkers(i).IDStr = 'Amplifier Reconnect';
case HCEVENT_PAUSED
obj.eventMarkers(i).IDStr = 'Recording Paused';
otherwise
obj.eventMarkers(i).IDStr = 'UNKNOWN';
end
% Next packet
offset = offset + pktLen;
fseek(h,offset,'bof');
pktGUID = fread(h,16,'uint8');
pktLen = fread(h,1,'uint64');
end
%% Get montage - Andrei Barborica, Dec 2015
% Derivation (montage)
mtgIdx = Tags(find(strcmp({Tags.IDStr},'DERIVATIONGUID'),1)).index;
indexIdx = find([obj.index.sectionIdx]==mtgIdx,1);
fseek(h,obj.index(indexIdx(1)).offset + 40,'bof'); % Beginning of current montage name
mtgName = deblank(char(fread(h,32,'uint16')).');
fseek(h,640,'cof'); % Number of traces in the montage
numDerivations = fread(h,1,'uint32');
numDerivations2 = fread(h,1,'uint32');
obj.montage = struct();
for i = 1:numDerivations
obj.montage(i).derivationName = deblank(char(fread(h,64,'uint16')).');
obj.montage(i).signalName1 = deblank(char(fread(h,32,'uint16')).');
obj.montage(i).signalName2 = deblank(char(fread(h,32,'uint16')).');
fseek(h,264,'cof'); % Skip additional info
end
% Display properties
dispIdx = Tags(find(strcmp({Tags.IDStr},'DISPLAYGUID'),1)).index;
indexIdx = find([obj.index.sectionIdx]==dispIdx,1);
fseek(h,obj.index(indexIdx(1)).offset + 40,'bof'); % Beginning of current montage name
displayName = deblank(char(fread(h,32,'uint16')).');
fseek(h,640,'cof'); % Number of traces in the montage
numTraces = fread(h,1,'uint32');
numTraces2 = fread(h,1,'uint32');
if (numTraces == numDerivations)
for i = 1:numTraces
fseek(h,32,'cof');
obj.montage(i).color = fread(h,1,'uint32'); % Use typecast(uint32(montage(i).color),'uint8') to convert to RGB array
fseek(h,136-4,'cof');
end
else
disp('Could not match montage derivations with display color table');
end
% Close File
fclose(h);
end
function out = getNrSamples(obj, segment)
% GETNRSAMPLES Returns the number of samples per channel in segment.
%
% OUT = GETNRSAMPLES(OBJ, SEGMENT) returns a 1xn array of values
% indicating the number of samples for each of the channels in the
% associated SEGMENT, where SEGMENT is the index of the
% OBJ.segments array.
assert(length(obj.segments)>= segment, ...
'Incorrect SEGMENT argument; must be integer representing segment index.');
out = obj.segments(segment).samplingRate .* obj.segments(segment).duration;
end
function cSumSegs = getCSumSegs(obj, chI)
% GETCSUMSEGS Returns the cumulative sum of a channels segments
%
% CSUMSEGS = GETCSUMSEGS(OBJ, CHI) returns a 1xn array of the
% cumulative sum of the given channel, where chI is the channel
% index. Different TsInfos can have different numbers of channels
% and different sampling rates.
samplingRates = zeros(1, length(obj.tsInfos));
for ts_i = 1:length(obj.tsInfos)
if length(obj.tsInfos{ts_i}) < chI
continue
else
samplingRates(ts_i) = obj.tsInfos{ts_i}(chI).dSamplingRate;
end
end
cSumSegs = [0 cumsum(samplingRates.*[obj.segments.duration])];
end
function out = getdata(obj, segment, range, chIdx)
% GETDATA Returns data from Nicolet file.
%
% OUT = GETDATA(OBJ, SEGMENT, RANGE, CHIDX) returns data in an nxm array of
% doubles where n is the number of datapoints and m is the number
% of channels. RANGE is a 1x2 array with the [StartIndex EndIndex]
% and CHIDX is a vector of channel indeces.
% Assert range is 1x2 vector
if isempty([obj.tsInfos{segment}.label])
warning('Segment %d has an empty tsInfo. Skipping', segment)
out = 0;
return
end
assert(length(range) == 2, 'Range is [firstIndex lastIndex]');
assert(length(segment) == 1, 'Segment must be single value.');
% Reopen .e file.
h = fopen(obj.fileName,'r','ieee-le');
% Find sectionID for channels
lChIdx = length(chIdx);
sectionIdx = zeros(lChIdx,1);
for i = 1:lChIdx
tmp = find(strcmp(num2str(chIdx(i)-1),{obj.sections.tag}),1);
sectionIdx(i) = obj.sections(tmp).index;
end
% Iterate over all requested channels and populate array.
out = zeros(range(2) - range(1) + 1, lChIdx);
for i = 1 : lChIdx
% Get cumulative sum segments.
cSumSegments = obj.getCSumSegs(chIdx(i));
% Get sampling rate for current channel
curSF = obj.segments(segment).samplingRate(chIdx(i));
mult = obj.segments(segment).scale(chIdx(i));
% Find all sections
allSectionIdx = obj.allIndexIDs == sectionIdx(i);
allSections = find(allSectionIdx);
% Find relevant sections
sectionLengths = [obj.index(allSections).sectionL]./2;
cSectionLengths = [0 cumsum(sectionLengths)];
skipValues = cSumSegments(segment);
firstSectionForSegment = find(cSectionLengths > skipValues, 1) - 1 ;
lastSectionForSegment = firstSectionForSegment + ...
find(cSectionLengths > curSF*obj.segments(segment).duration,1) - 2 ;
if isempty(lastSectionForSegment)
lastSectionForSegment = length(cSectionLengths);
end
offsetSectionLengths = cSectionLengths - cSectionLengths(firstSectionForSegment);
firstSection = find(offsetSectionLengths < range(1) ,1,'last');
lastSection = find(offsetSectionLengths >= range(2),1)-1;
if isempty(lastSection)
lastSection = length(offsetSectionLengths);
end
if lastSection > lastSectionForSegment
error('Index out of range for current section: %i > %i, on channel: %i', ...
range(2), cSectionLengths(lastSectionForSegment+1), chIdx(i));
end
useSections = allSections(firstSection: lastSection) ;
useSectionL = sectionLengths(firstSection: lastSection) ;
% First Partial Segment
curIdx = 1;
curSec = obj.index(useSections(1));
fseek(h, curSec.offset,'bof');
firstOffset = range(1) - offsetSectionLengths(firstSection);
lastOffset = min([range(2) useSectionL(1)]);
lsec = lastOffset-firstOffset + 1;
fseek(h, (firstOffset-1) * 2,'cof');
out(1 : lsec,i) = fread(h, lsec, 'int16') * mult;
curIdx = curIdx + lsec;
if length(useSections) > 1
% Full Segments
for j = 2: (length(useSections)-1)
curSec = obj.index(useSections(j));
fseek(h, curSec.offset,'bof');
out(curIdx : (curIdx + useSectionL(j) - 1),i) = ...
fread(h, useSectionL(j), 'int16') * mult;
curIdx = curIdx + useSectionL(j);
end
% Final Partial Segment
curSec = obj.index(useSections(end));
fseek(h, curSec.offset,'bof');
out(curIdx : end,i) = fread(h, length(out)-curIdx + 1, 'int16') * mult;
end
end
% Close the .e file.
fclose(h);
end
function out = getdataQ(obj, segment, range, chIdx)
% GETDATAQ Returns data from Nicolet file. This is a "QUICK" version of getdata,
% that uses more memory but operates faster on large datasets by reading
% a single block of data from disk that contains all data of interest.
%
% OUT = GETDATAQ(OBJ, SEGMENT, RANGE, CHIDX) returns data in an nxm array of
% doubles where n is the number of datapoints and m is the number
% of channels. RANGE is a 1x2 array with the [StartIndex EndIndex]
% and CHIDX is a vector of channel indeces.
%
% Andrei Barborica, Dec 2015
%
if isempty([obj.tsInfos{segment}.label])
warning('Segment %d has an empty tsInfo. Skipping', segment)
out = 0;
return
end
% Assert range is 1x2 vector
assert(length(range) == 2, 'Range is [firstIndex lastIndex]');
assert(length(segment) == 1, 'Segment must be single value.');
% Reopen .e file.
h = fopen(obj.fileName,'r','ieee-le');
% Find sectionID for channels
lChIdx = length(chIdx);
sectionIdx = zeros(lChIdx,1);
for i = 1:lChIdx
tmp = find(strcmp(num2str(chIdx(i)-1),{obj.sections.tag}),1);
sectionIdx(i) = obj.sections(tmp).index;
end
usedIndexEntries = zeros(size([obj.index.offset]));
% Iterate over all requested channels and populate array.
out = zeros(range(2) - range(1) + 1, lChIdx);
for i = 1 : lChIdx
% Get cumulative sum segments.
cSumSegments = obj.getCSumSegs(chIdx(i));
% Get sampling rate for current channel
curSF = obj.segments(segment).samplingRate(chIdx(i));
mult = obj.segments(segment).scale(chIdx(i));
% Find all sections
allSectionIdx = obj.allIndexIDs == sectionIdx(i);
allSections = find(allSectionIdx);
% Find relevant sections
sectionLengths = [obj.index(allSections).sectionL]./2;
cSectionLengths = [0 cumsum(sectionLengths)];
skipValues = cSumSegments(segment);
firstSectionForSegment = find(cSectionLengths > skipValues, 1) - 1 ;
lastSectionForSegment = firstSectionForSegment + ...
find(cSectionLengths > curSF*obj.segments(segment).duration,1) - 2 ;
if isempty(lastSectionForSegment)
lastSectionForSegment = length(cSectionLengths);
end
offsetSectionLengths = cSectionLengths - cSectionLengths(firstSectionForSegment);
firstSection = find(offsetSectionLengths < range(1) ,1,'last');
lastSection = find(offsetSectionLengths >= range(2),1)-1;
if isempty(lastSection)
lastSection = length(offsetSectionLengths);
end
if lastSection > lastSectionForSegment
error('Index out of range for current section: %i > %i, on channel: %i', ...
range(2), cSectionLengths(lastSectionForSegment+1), chIdx(i));
end
useSections = allSections(firstSection: lastSection) ;
useSectionL = sectionLengths(firstSection: lastSection) ;
% First Partial Segment
usedIndexEntries(useSections(1)) = 1;
if length(useSections) > 1
% Full Segments
for j = 2: (length(useSections)-1)
usedIndexEntries(useSections(j)) = 1;
end
% Final Partial Segment
usedIndexEntries(useSections(end)) = 1;
end
end
% Read a big chunk of the file, containing data of interest.
ix = find(usedIndexEntries);
fseek(h, obj.index(ix(1)).offset,'bof');
dsize = obj.index(ix(end)).offset - obj.index(ix(1)).offset + obj.index(ix(end)).sectionL;
tmp = fread(h,dsize/2,'int16').';
% Close the .e file.
fclose(h);
baseOffset = obj.index(ix(1)).offset;
% Extract specified channels
for i = 1 : lChIdx
% Get cumulative sum segments.
cSumSegments = obj.getCSumSegs(chIdx(i));
% Get sampling rate for current channel
curSF = obj.segments(segment).samplingRate(chIdx(i));
mult = obj.segments(segment).scale(chIdx(i));
% Find all sections
allSectionIdx = obj.allIndexIDs == sectionIdx(i);
allSections = find(allSectionIdx);
% Find relevant sections
sectionLengths = [obj.index(allSections).sectionL]./2;
cSectionLengths = [0 cumsum(sectionLengths)];
skipValues = cSumSegments(segment);
firstSectionForSegment = find(cSectionLengths > skipValues, 1) - 1 ;
lastSectionForSegment = firstSectionForSegment + ...
find(cSectionLengths > curSF*obj.segments(segment).duration,1) - 2 ;
if isempty(lastSectionForSegment)
lastSectionForSegment = length(cSectionLengths);
end
offsetSectionLengths = cSectionLengths - cSectionLengths(firstSectionForSegment);
firstSection = find(offsetSectionLengths < range(1) ,1,'last');
lastSection = find(offsetSectionLengths >= range(2),1)-1;
if isempty(lastSection)
lastSection = length(offsetSectionLengths);
end
if lastSection > lastSectionForSegment
error('Index out of range for current section: %i > %i, on channel: %i', ...
range(2), cSectionLengths(lastSectionForSegment+1), chIdx(i));
end
useSections = allSections(firstSection: lastSection) ;
useSectionL = sectionLengths(firstSection: lastSection) ;
% First Partial Segment
curIdx = 1;
curSec = obj.index(useSections(1));
%fseek(h, curSec.offset,'bof');
firstOffset = range(1) - offsetSectionLengths(firstSection);
lastOffset = min([range(2) useSectionL(1)]);
lsec = lastOffset-firstOffset + 1;
out(1 : lsec,i) = tmp( (curSec.offset - baseOffset)/2 + (firstOffset-1) + (1:lsec) ) * mult;
curIdx = curIdx + lsec;
if length(useSections) > 1
% Full Segments
for j = 2: (length(useSections)-1)
curSec = obj.index(useSections(j));
out(curIdx : (curIdx + useSectionL(j) - 1),i) = ...
tmp( (curSec.offset - baseOffset)/2 + (1:useSectionL(j)) ) * mult;
curIdx = curIdx + useSectionL(j);
end
% Final Partial Segment
curSec = obj.index(useSections(end));
out(curIdx : end,i) = tmp( (curSec.offset - baseOffset)/2 + (1:(length(out)-curIdx + 1)) ) * mult; % length(out) ??????
end
end
end
function labels = getlabels(obj,str)
% Returns annotations containing specified string
%
% Cristian Donos, Dec 2015
%
labels=[]; counter = 1;
for i = 1:length(obj.eventMarkers)
if strfind(lower(obj.eventMarkers(i).annotation),lower(str))
labels{counter,1} = obj.eventMarkers(i).annotation; % annotation string
labels{counter,2} = i; % annotation index in obj.eventMarkers
% identify segment
time_vector = [];
for j = 1:length(obj.segments)
time_vector = [time_vector etime(datevec(obj.eventMarkers(i).dateStr),datevec(obj.segments(j).dateStr))];
end
labels{counter,3}= find(time_vector==min(time_vector(time_vector>0))); % annotation part of this segment
labels{counter,4}= min(time_vector(time_vector>0)); % annotation offset in seconds, relative to its segment start
counter = counter+1;
end
end
end
end
end
obj = NicoletFile('janbrogger.e');
% save data
out = getdata(obj, 1, [1 20000], 1:16);
plot(out)
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
