edfread源码

  1 function [hdr, record] = edfread(fname, varargin)
  2 % Read European Data Format file into MATLAB
  3 %
  4 % [hdr, record] = edfread(fname)
  5 %         Reads data from ALL RECORDS of file fname ('*.edf'). Header
  6 %         information is returned in structure hdr, and the signals
  7 %         (waveforms) are returned in structure record, with waveforms
  8 %         associated with the records returned as fields titled 'data' of
  9 %         structure record.
 10 %
 11 % [...] = edfread(fname, 'assignToVariables', assignToVariables)
 12 %         Triggers writing of individual output variables, as defined by
 13 %         field 'labels', into the caller workspace.
 14 %
 15 % [...] = edfread(...,'desiredSignals',desiredSignals)
 16 %         Allows user to specify the names (or position numbers) of the
 17 %         subset of signals to be read. |desiredSignals| may be either a
 18 %         string, a cell array of comma-separated strings, or a vector of
 19 %         numbers. (Default behavior is to read all signals.)
 20 %         E.g.:
 21 %         data = edfread(mydata.edf,'desiredSignals','Thoracic');
 22 %         data = edfread(mydata.edf,'desiredSignals',{'Thoracic1','Abdominal'});
 23 %         or
 24 %         data = edfread(mydata.edf,'desiredSignals',[2,4,6:13]);
 25 %
 26 % FORMAT SPEC: Source: http://www.edfplus.info/specs/edf.html SEE ALSO:
 27 % http://www.dpmi.tu-graz.ac.at/~schloegl/matlab/eeg/edf_spec.htm
 28 %
 29 % The first 256 bytes of the header record specify the version number of
 30 % this format, local patient and recording identification, time information
 31 % about the recording, the number of data records and finally the number of
 32 % signals (ns) in each data record. Then for each signal another 256 bytes
 33 % follow in the header record, each specifying the type of signal (e.g.
 34 % EEG, body temperature, etc.), amplitude calibration and the number of
 35 % samples in each data record (from which the sampling frequency can be
 36 % derived since the duration of a data record is also known). In this way,
 37 % the format allows for different gains and sampling frequencies for each
 38 % signal. The header record contains 256 + (ns * 256) bytes.
 39 %
 40 % Following the header record, each of the subsequent data records contains
 41 % 'duration' seconds of 'ns' signals, with each signal being represented by
 42 % the specified (in the header) number of samples. In order to reduce data
 43 % size and adapt to commonly used software for acquisition, processing and
 44 % graphical display of polygraphic signals, each sample value is
 45 % represented as a 2-byte integer in 2's complement format. Figure 1 shows
 46 % the detailed format of each data record.
 47 %
 48 % DATA SOURCE: Signals of various types (including the sample signal used
 49 % below) are available from PHYSIONET: http://www.physionet.org/
 50 %
 51 %
 52 % % EXAMPLE 1:
 53 % % Read all waveforms/data associated with file 'ecgca998.edf':
 54 %
 55 % [header, recorddata] = edfRead('ecgca998.edf');
 56 %
 57 % % EXAMPLE 2:
 58 % % Read records 3 and 5, associated with file 'ecgca998.edf':
 59 %
 60 % header = edfRead('ecgca998.edf','AssignToVariables',true);
 61 % % Header file specifies data labels 'label_1'...'label_n'; these are
 62 % % created as variables in the caller workspace.
 63 %
 64 % Coded 8/27/09 by Brett Shoelson, PhD
 65 % brett.shoelson@mathworks.com
 66 % Copyright 2009 - 2012 MathWorks, Inc.
 67 %
 68 % Modifications:
 69 % 5/6/13 Fixed a problem with a poorly subscripted variable. (Under certain
 70 % conditions, data were being improperly written to the 'records' variable.
 71 % Thanks to Hisham El Moaqet for reporting the problem and for sharing a
 72 % file that helped me track it down.)
 73 % 
 74 % 5/22/13 Enabled import of a user-selected subset of signals. Thanks to
 75 % Farid and Cindy for pointing out the deficiency. Also fixed the import of
 76 % signals that had "bad" characters (spaces, etc) in their names.
 77 
 78 % HEADER RECORD
 79 % 8 ascii : version of this data format (0)
 80 % 80 ascii : local patient identification
 81 % 80 ascii : local recording identification
 82 % 8 ascii : startdate of recording (dd.mm.yy)
 83 % 8 ascii : starttime of recording (hh.mm.ss)
 84 % 8 ascii : number of bytes in header record
 85 % 44 ascii : reserved
 86 % 8 ascii : number of data records (-1 if unknown)
 87 % 8 ascii : duration of a data record, in seconds
 88 % 4 ascii : number of signals (ns) in data record
 89 % ns * 16 ascii : ns * label (e.g. EEG FpzCz or Body temp)
 90 % ns * 80 ascii : ns * transducer type (e.g. AgAgCl electrode)
 91 % ns * 8 ascii : ns * physical dimension (e.g. uV or degreeC)
 92 % ns * 8 ascii : ns * physical minimum (e.g. -500 or 34)
 93 % ns * 8 ascii : ns * physical maximum (e.g. 500 or 40)
 94 % ns * 8 ascii : ns * digital minimum (e.g. -2048)
 95 % ns * 8 ascii : ns * digital maximum (e.g. 2047)
 96 % ns * 80 ascii : ns * prefiltering (e.g. HP:0.1Hz LP:75Hz)
 97 % ns * 8 ascii : ns * nr of samples in each data record
 98 % ns * 32 ascii : ns * reserved
 99 
100 % DATA RECORD
101 % nr of samples[1] * integer : first signal in the data record
102 % nr of samples[2] * integer : second signal
103 % ..
104 % ..
105 % nr of samples[ns] * integer : last signal
106 
107 if nargin > 5
108     error('EDFREAD: Too many input arguments.');
109 end
110 
111 if ~nargin
112     error('EDFREAD: Requires at least one input argument (filename to read).');
113 end
114 
115 [fid,msg] = fopen(fname,'r');
116 if fid == -1
117     error(msg)
118 end
119 
120 assignToVariables = false; %Default
121 targetSignals = []; %Default
122 for ii = 1:2:numel(varargin)
123     switch lower(varargin{ii})
124         case 'assigntovariables'
125             assignToVariables = varargin{ii+1};
126         case 'targetsignals'
127             targetSignals = varargin{ii+1};
128         otherwise
129             error('EDFREAD: Unrecognized parameter-value pair specified. Valid values are ''assignToVariables'' and ''targetSignals''.')
130     end
131 end
132 
133 
134 % HEADER
135 hdr.ver        = str2double(char(fread(fid,8)'));
136 hdr.patientID  = fread(fid,80,'*char')';
137 hdr.recordID   = fread(fid,80,'*char')';
138 hdr.startdate  = fread(fid,8,'*char')';% (dd.mm.yy)
139 % hdr.startdate  = datestr(datenum(fread(fid,8,'*char')','dd.mm.yy'), 29); %'yyyy-mm-dd' (ISO 8601)
140 hdr.starttime  = fread(fid,8,'*char')';% (hh.mm.ss)
141 % hdr.starttime  = datestr(datenum(fread(fid,8,'*char')','hh.mm.ss'), 13); %'HH:MM:SS' (ISO 8601)
142 hdr.bytes      = str2double(fread(fid,8,'*char')');
143 reserved       = fread(fid,44);
144 hdr.records    = str2double(fread(fid,8,'*char')');
145 hdr.duration   = str2double(fread(fid,8,'*char')');
146 % Number of signals
147 hdr.ns    = str2double(fread(fid,4,'*char')');
148 for ii = 1:hdr.ns
149     hdr.label{ii} = regexprep(fread(fid,16,'*char')','W','');
150 end
151 
152 if isempty(targetSignals)
153     targetSignals = 1:numel(hdr.label);
154 elseif iscell(targetSignals)||ischar(targetSignals)
155     targetSignals = find(ismember(hdr.label,regexprep(targetSignals,'W','')));
156 end
157 if isempty(targetSignals)
158     error('EDFREAD: The signal(s) you requested were not detected.')
159 end
160 
161 for ii = 1:hdr.ns
162     hdr.transducer{ii} = fread(fid,80,'*char')';
163 end
164 % Physical dimension
165 for ii = 1:hdr.ns
166     hdr.units{ii} = fread(fid,8,'*char')';
167 end
168 % Physical minimum
169 for ii = 1:hdr.ns
170     hdr.physicalMin(ii) = str2double(fread(fid,8,'*char')');
171 end
172 % Physical maximum
173 for ii = 1:hdr.ns
174     hdr.physicalMax(ii) = str2double(fread(fid,8,'*char')');
175 end
176 % Digital minimum
177 for ii = 1:hdr.ns
178     hdr.digitalMin(ii) = str2double(fread(fid,8,'*char')');
179 end
180 % Digital maximum
181 for ii = 1:hdr.ns
182     hdr.digitalMax(ii) = str2double(fread(fid,8,'*char')');
183 end
184 for ii = 1:hdr.ns
185     hdr.prefilter{ii} = fread(fid,80,'*char')';
186 end
187 for ii = 1:hdr.ns
188     hdr.samples(ii) = str2double(fread(fid,8,'*char')');
189 end
190 for ii = 1:hdr.ns
191     reserved    = fread(fid,32,'*char')';
192 end
193 hdr.label = hdr.label(targetSignals);
194 hdr.label = regexprep(hdr.label,'W','');
195 hdr.units = regexprep(hdr.units,'W','');
196 disp('Step 1 of 2: Reading requested records. (This may take a few minutes.)...');
197 if nargout > 1 || assignToVariables
198     % Scale data (linear scaling)
199     scalefac = (hdr.physicalMax - hdr.physicalMin)./(hdr.digitalMax - hdr.digitalMin);
200     dc = hdr.physicalMax - scalefac .* hdr.digitalMax;
201     
202     % RECORD DATA REQUESTED
203     tmpdata = struct;
204     for recnum = 1:hdr.records
205         for ii = 1:hdr.ns
206             % Read or skip the appropriate number of data points
207             if ismember(ii,targetSignals)
208                 % Use a cell array for DATA because number of samples may vary
209                 % from sample to sample
210                 tmpdata(recnum).data{ii} = fread(fid,hdr.samples(ii),'int16') * scalefac(ii) + dc(ii);
211             else
212                 fseek(fid,hdr.samples(ii)*2,0);
213             end
214         end
215     end
216     hdr.units = hdr.units(targetSignals);
217     hdr.physicalMin = hdr.physicalMin(targetSignals);
218     hdr.physicalMax = hdr.physicalMax(targetSignals);
219     hdr.digitalMin = hdr.digitalMin(targetSignals);
220     hdr.digitalMax = hdr.digitalMax(targetSignals);
221     hdr.prefilter = hdr.prefilter(targetSignals);
222     hdr.transducer = hdr.transducer(targetSignals);
223     
224     record = zeros(numel(hdr.label), hdr.samples(1)*hdr.records);
225     % NOTE: 5/6/13 Modified for loop below to change instances of hdr.samples to
226     % hdr.samples(ii). I think this underscored a problem with the reader.
227     
228     disp('Step 2 of 2: Parsing data...');
229     recnum = 1;
230     for ii = 1:hdr.ns
231         if ismember(ii,targetSignals)
232             ctr = 1;
233             for jj = 1:hdr.records
234                 try
235                     record(recnum, ctr : ctr + hdr.samples(ii) - 1) = tmpdata(jj).data{ii};
236                 end
237                 ctr = ctr + hdr.samples(ii);
238             end
239             recnum = recnum + 1;
240         end
241     end
242     hdr.ns = numel(hdr.label);
243     hdr.samples = hdr.samples(targetSignals);
244 
245     if assignToVariables
246         for ii = 1:numel(hdr.label)
247             try
248                 eval(['assignin(''caller'',''',hdr.label{ii},''',record(ii,:))'])
249             end
250         end
251         record = [];
252     end
253 end
254 fclose(fid);
原文地址:https://www.cnblogs.com/myohao/p/matlab.html