sdfopen.m

matlab数字信号处理工具箱

function [EDF,H1,h2]=sdfopen(arg1,arg2,arg3,arg4,arg5,arg6)
% Opens EDF/GDF/SDF files for reading and writing. 
% The EDF format is specified in [1], GDF in [2].
% SDF is the SIESTA convention about polygraphic recordings (see [3], pp.8-9). 
% SDF used the EDF format.
%
% EDF=sdfopen(Filename,'r' [,CHAN [,Mode [,TSR]]]);
% [S,EDF]=sdfread(EDF, NoR, StartPos)
%
% CHAN defines the selected Channels or a re-referencing matrix
%
% Mode=='SIESTA' indicates that #1 - #6 is rereferenced to (M1+M2)/2           
% Mode=='AFIR' indicates that Adaptive FIR filtering is used for ECG removing
%			Implements Adaptive FIR filtering for ECG removal in EDF/GDF-tb.
% 			based on the Algorithm of Mikko Koivuluoma <[e-mail]>
%                 A delay of EDF.AFIR.delay number of samples has to be considered. 
% Mode=='SIESTA+AFIR' for both
% Mode=='UCAL' [default]
%                 indicates that no calibration (re-scaling) to Physical Dim. is used. 
%                 The output are 16bit interger numbers. 'UCAL' overrides 'SIESTA'
% Mode=='RAW' One column represents one EDF-block
% Mode=='Notch50Hz' Implements a simple FIR-notch filter for 50Hz
% Mode=='RECG' Implements ECG minimization with regression analysis
% Mode=='TECG' Implements ECG minimization with template removal (Test status)
% Mode=='HPF0.0Hz'  Implements a high-pass filter (with the zero at z=+1, i.e. a differentiator 
%                     In this case a Notch-filter and/or sub-sampling is recommended. 
% Mode=='TAUx.yS' compensates time-constant of x.y seconds
% Mode=='EOG[hvr]' produces HEOG, VEOG and/or REOG output (CHAN not considered)
% 
% Mode=='OVERFLOW' overflow detection
% Mode=='FailingElectrodeDetector' using optimized threshold based AUC for FED & OFC
% Mode=='Units_Blocks' requests the arguments in SDFREAD in Blocks [default is seconds]
%
% TSR [optional] is the target sampling rate
%         Currently only downsampling from 256 and 200 to 100Hz is supported.  
%         The details are described in the appendix of [4].
%
% EDF.ErrNo~=0  indicates that an error occured.
% For compatibility to former versions, 
%    EDF.FILE.FID = -1 indicates that file has not been opened.
%
% Opening of an EDF/SDF File for writing
% [EDF]=sdfopen(EDF,'w') is equal to  
% [EDF]=sdfopen(EDF.FileName,'w',EDF.Dur,EDF.SampleRate);
% At least EDF.FileName, EDF.NS, EDF.Dur and EDF.EDF.SampleRate must be defined
% 
% 
% See also: fopen, SDFREAD, SWRITE, SCLOSE, SSEEK, SREWIND, STELL, SEOF


% References: 
% [1] Bob Kemp, Alpo V鋜ri, Agostinho C. Rosa, Kim D. Nielsen and John Gade
%     A simple format for exchange of digitized polygraphic recordings.
%     Electroencephalography and Clinical Neurophysiology, 82 (1992) 391-393.
% see also: http://www.medfac.leidenuniv.nl/neurology/knf/kemp/edf/edf_spec.htm
%
% [2] Alois Schl鰃l, Oliver Filz, Herbert Ramoser, Gert Pfurtscheller
%     GDF - A GENERAL DATAFORMAT FOR BIOSIGNALS
%     Technical Report, Department for Medical Informatics, Universtity of Technology, Graz (2004)
% see also: http://www.dpmi.tu-graz.ac.at/~schloegl/matlab/eeg/gdf4/TR_GDF.PDF
%
% [3] The SIESTA recording protocol. 
% see http://www.ai.univie.ac.at/siesta/protocol.html
% and http://www.ai.univie.ac.at/siesta/protocol.rtf 
%
% [4] Alois Schl鰃l
%     The electroencephalogram and the adaptive autoregressive model: theory and applications. 
%     (ISBN 3-8265-7640-3) Shaker Verlag, Aachen, Germany.
% see also: "http://www.shaker.de/Online-Gesamtkatalog/Details.idc?ISBN=3-8265-7640-3"


% Testing state
%
% (1) reads header information and writes the header; can be used to check SDFOPEN or for correcting headerinformation
% EDF=sdfopen(EDF,'r+b'); EDF=sdfclose(EDF); 
% 
% (2a) Minimal requirements for opening an EDF-File
% EDF.FileName='file.edf'; % define filename
% EDF.NS = 5; % fix number of channels
% EDF=sdfopen(EDF,'wb');
%     write file
%     define header somewhen before 
% EDF=sdfclose(EDF); % writes the corrected header information
% 
% (2b) Minimal requirements for opening an EDF-File
% EDF=sdfopen('file.edf','wb',N); % N number of channels
%      .. do anything, e.g define header
% EDF=sdfopen(EDF,'w+b'); % writes Header information
%
%
% 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.

% EDF.WarnNo   1: ascii(0) in 1st header
%              2: ascii(0) in 2nd header
%              4: invalid SPR
%              4: invalid samples_per_record-values
%              8: date not in EDF-format (tries to guess correct date, see also E2)
%             16: invalid number_of-channels-value
%             32: invalid value of the EDF header length
%             64: invalid value of block_duration
%            128: Polarity of #7 probably inverted  
%
% EDF.ErrNo    1: first 8 bytes different to '0       ', this violates the EDF spec.
%              2: Invalid date (unable to guess correct date)
%              4: Incorrect date information (later than actual date) 
%             16: incorrect filesize, Header information does not match actual size

%	$Revision: 1.22 $
%	$Id: sdfopen.m,v 1.22 2004/08/19 00:36:40 schloegl Exp $
INFO='(C) 1997-2002, 2004 by Alois Schloegl, 17 Aug 2004 #0.87';
%	[e-mail]

if nargin<2, 
        arg2='rb'; 
elseif ~any(arg2=='b');
        arg2= [arg2,'b']; % force binary open. 
end;

if 1,exist('OCTAVE_VERSION');
        EDF.AS.Method='sdfopen';
        mfilename=EDF.AS.Method;
else 
        EDF.AS.Method=mfilename;
end;


EDF.AS.Method=[EDF.AS.Method '-' arg2];
EDF.AS.Date=fix(clock);
EDF.AS.Info=INFO;
EDF.AS.Ver = 0.86;

if isstruct(arg1) 
        EDF=arg1; 
        FILENAME=EDF.FileName;
else
        FILENAME=arg1;
end;

H1idx=[8 80 80 8 8 8 44 8 8 4];
H2idx=[16 80 8 8 8 8 8 80 8 32];

%%%%% Define Valid Data types %%%%%%
%GDFTYPES=[0 1 2 3 4 5 6 7 16 17 255+(1:64) 511+(1:64)];
GDFTYPES=[0 1 2 3 4 5 6 7 16 17 255+[1 12 22 24] 511+[1 12 22 24]];

%%%%% Define Size for each data type %%%%%
GDFTYP_BYTE=zeros(1,512+64);
GDFTYP_BYTE(256+(1:64))=(1:64)/8;
GDFTYP_BYTE(512+(1:64))=(1:64)/8;
GDFTYP_BYTE(1:18)=[1 1 1 2 2 4 4 8 8 4 8 0 0 0 0 0 4 8]';

%EDF.GDFTYP.TEXT={'char','int8','uint8','int16','uint16','int32','uint32','int64','uint64','float32','float64'};
%GDFTYP_BYTE=[1 1 1 2 2 4 4 8 8 4 8 0 0 0 0 0 4 8]';
%GDFTYPES=[0 1 2 3 4 5 6 7 16 17];

EDF.ErrNo = 0;

%%%%%%% ============= READ ===========%%%%%%%%%%%%
if any(arg2=='r'), 
        
[EDF.FILE.FID,MESSAGE]=fopen(FILENAME,arg2,'ieee-le');          
%EDF.FILE.FID=fid;

EDF.FILE.stderr=2;
if EDF.FILE.FID<0 
        %fprintf(EDF.FILE.stderr,'Error SDFOPEN: %s %s\n',MESSAGE,FILENAME);  
        H1=MESSAGE; H2=FILENAME;
        EDF.ErrNo = [32,EDF.ErrNo];
	return;
end;
EDF.FILE.OPEN = 1 + any(arg2=='+'); 
EDF.FileName = FILENAME;

PPos=min([max(find(FILENAME=='.')) length(FILENAME)+1]);
SPos=max([0 find(FILENAME==filesep)]);
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 filesep EDF.FILE.Name '.' EDF.FILE.Ext];
%EDF.datatype = 'EDF';

%%% Read Fixed Header %%%
[tmp,count]=fread(EDF.FILE.FID,184,'uchar');     %
if count<184,
        EDF.ErrNo = [64,EDF.ErrNo];
        return;
end;
H1=setstr(tmp');     %


EDF.VERSION=H1(1:8);                     % 8 Byte  Versionsnummer 
if ~(strcmp(EDF.VERSION,'0       ') | all(abs(EDF.VERSION)==[255,abs('BIOSEMI')]) | strcmp(EDF.VERSION(1:3),'GDF'))
        EDF.ErrNo = [1,EDF.ErrNo];
	if ~strcmp(EDF.VERSION(1:3),'   '); % if not a scoring file, 
%	    return; 
	end;
end;
EDF.PID = deblank(H1(9:88));                  % 80 Byte local patient identification
EDF.RID = deblank(H1(89:168));                % 80 Byte local recording identification

IsGDF=strcmp(EDF.VERSION(1:3),'GDF');

if strcmp(EDF.VERSION(1:3),'GDF'),
        if 1, % if strcmp(EDF.VERSION(4:8),' 0.12'); % in future versions the date format might change. 
      		EDF.T0(1,1) = str2double( H1(168 + [ 1:4]));
		EDF.T0(1,2) = str2double( H1(168 + [ 5 6]));
        	EDF.T0(1,3) = str2double( H1(168 + [ 7 8]));
        	EDF.T0(1,4) = str2double( H1(168 + [ 9 10]));
        	EDF.T0(1,5) = str2double( H1(168 + [11 12]));
        	EDF.T0(1,6) = str2double( H1(168 + [13:16]))/100;
     	end; 
     
	if str2double(EDF.VERSION(4:8))<0.12
                tmp = setstr(fread(EDF.FILE.FID,8,'uchar')');    % 8 Byte  Length of Header
                EDF.HeadLen = str2double(tmp);    % 8 Byte  Length of Header
        else
		%EDF.HeadLen = fread(EDF.FILE.FID,1,'int64');    % 8 Byte  Length of Header
		EDF.HeadLen = fread(EDF.FILE.FID,1,'int32');    % 8 Byte  Length of Header
		tmp         = fread(EDF.FILE.FID,1,'int32');    % 8 Byte  Length of Header
        end;
        EDF.reserved1 = fread(EDF.FILE.FID,8+8+8+20,'uchar');     % 44 Byte reserved
        
        %EDF.NRec = fread(EDF.FILE.FID,1,'int64');     % 8 Byte  # of data records
        EDF.NRec = fread(EDF.FILE.FID,1,'int32');     % 8 Byte  # of data records
                   fread(EDF.FILE.FID,1,'int32');     % 8 Byte  # of data records
        if strcmp(EDF.VERSION(4:8),' 0.10')
                EDF.Dur =  fread(EDF.FILE.FID,1,'float64');    % 8 Byte  # duration of data record in sec
        else
                tmp =  fread(EDF.FILE.FID,2,'uint32');    % 8 Byte  # duration of data record in sec
                EDF.Dur =  tmp(1)./tmp(2);
        end;
        EDF.NS =   fread(EDF.FILE.FID,1,'uint32');     % 4 Byte  # of signals
else 
        if exist('OCTAVE_VERSION')>=5
                tmp=(find((toascii(H1)<32) | (toascii(H1)>126))); 	%%% snytax for OCTAVE
        else
                tmp=(find((H1<32) | (H1>126))); 		%%% syntax for Matlab
        end;
        if ~isempty(tmp) %%%%% not EDF because filled out with ASCII(0) - should be spaces
                %H1(tmp)=32; 
                EDF.ErrNo=[1025,EDF.ErrNo];
        end;
        
        EDF.T0 = zeros(1,6);
        ErrT0=0;
        tmp = str2double( H1(168 + [ 7  8]));
        if ~isnan(tmp), EDF.T0(1) = tmp; else ErrT0 = 1; end;
        tmp = str2double( H1(168 + [ 4  5]));
        if ~isnan(tmp), EDF.T0(2) = tmp; else ErrT0 = 1; end;
        tmp = str2double( H1(168 + [ 1  2]));
        if ~isnan(tmp), EDF.T0(3) = tmp; else ErrT0 = 1; end;
        tmp = str2double( H1(168 + [ 9 10]));
        if ~isnan(tmp), EDF.T0(4) = tmp; else ErrT0 = 1; end;
        tmp = str2double( H1(168 + [12 13]));
        if ~isnan(tmp), EDF.T0(5) = tmp; else ErrT0 = 1; end;
        tmp = str2double( H1(168 + [15 16]));
        if ~isnan(tmp), EDF.T0(6) = tmp; else ErrT0 = 1; end;
        
	if any(EDF.T0~=fix(EDF.T0)); ErrT0=1; end;

        if ErrT0,
                ErrT0=0;
                EDF.ErrNo = [1032,EDF.ErrNo];
                
                tmp = H1(168 + [1:8]);
                for k = [3 2 1],
                        %fprintf(1,'\n zz%szz \n',tmp);
                        [tmp1,tmp] = strtok(tmp,' :./-');
			tmp1 = str2double([tmp1,' ']);
			
                        if isempty(tmp1)
                                ErrT0 = ErrT0 | 1;
                        else
                                EDF.T0(k)  = tmp1;
                        end;
                end;
                tmp = H1(168 + [9:16]);
                for k = [4 5 6],
                        [tmp1,tmp] = strtok(tmp,' :./-');
                        tmp1=str2double([tmp1,' ']);
                        if isempty(tmp1)
                                ErrT0 = ErrT0 | 1;
                        else
                                EDF.T0(k)  = tmp1;
                        end;
                end;
                if ErrT0
                        EDF.ErrNo = [2,EDF.ErrNo];
                end;
        else