sdfopen.m

含有多种ICA算法的eeglab工具箱

% sdfopen() - Opens EDF/GDF/SDF files for reading and writing. The EDF format is specified 
%             in [1], GDF in [2]. SDF is the SIESTA convention (see [3], pp.8-9); SDF uses 
%             the EDF format.
% Usage:
%       >> EDF     = sdfopen(Filename,'r' [,CHAN]);  % prepare to read CHAN channels
%       >> EDF     = sdfopen(Filename,'r' [,CHAN [,MODE [,TSR]]]); % all arguments
% % then
%       >> [S,EDF] = sdfread(EDF, SecRead, SecStart); % reads the opened EDF data file.
%       >> [S,EDF] = sdfread(EDF, Inf); % reads the whole opened EDF data file.
%
% Inputs: 
%      CHAN - [int vector] Specifies the channel(s) to read. Else, a re-referencing 
%                  matrix. [Default|0 -> read all channels].
%      MODE  -
%
%   'UCAL' [Default mode] Indicates that no calibration (re-scaling) to physical dim. 
%                  is performed. Outputs are 16-bit integers. 
% 
%   'SIESTA' - Indicates that channels #1-#6 are re-referenced to (M1+M2)/2          
%                 (Note: 'UCAL' overrides 'SIESTA')
%   'AFIR'   - Indicates that Adaptive FIR filtering is used for ECG removal.
%		       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. 
%   'SIESTA+AFIR' - Does both.
%   'RAW'       - One column represents one EDF-block
%   'Notch50Hz' - Implements a simple FIR-notch filter at 50 Hz
%   'RECG'      - Implements ECG minimization with regression analysis
%   'TECG'      - Implements ECG minimization with template removal (test status)
%   'HPF0.0Hz'  - Implements a high-pass filter (with zero at z=+1, i.e. a differentiator).
%                 In this case, a notch-filter and/or sub-sampling is recommended. 
%   'TAUx.yS'   - Compensates time-constant of x.y seconds
%   'EOG[hvr]'  - Produces HEOG, VEOG and/or REOG output (CHAN not considered)
%   'OVERFLOW'  - Performs overflow detection
%   'Units_Blocks' - Requests the EDF-field arguments to SDFREAD in blocks 
%                    [default is seconds]
%      TSR - [optional] The target (re)sampling rate. Currently, only downsampling 
%            from 256 Hz or 200 Hz to 100 Hz is supported.  The details are described 
%            in the appendix of [4].
%
% Outputs: 
%         EDF - data structure read from the input file header.
%           EDF.ErrNo   ~= 0  Indicates that an error occurred 
%              1: First 8 bytes are not '0       ', violating the EDF spec.
%              2: Invalid date (unable to guess correct date)
%              4: Incorrect date information (later than current date) 
%             16: Incorrect filesize: Header information does not match actual size
%           EDF.FILE.FID = -1 indicates that file has not been opened
%                           (for compatibility with former versions).
%           EDF.WarnNo  ~=0 Indicates EDF structure problems
%              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  
%
% Example: 
%      To open an EDF/SDF file for writing:
%      >> [EDF] = sdfopen(EDF,'w') % or equivalently
%      >> [EDF] = sdfopen(EDF.FileName,'w',EDF.Dur,EDF.SampleRate);
%
% Note: Fields EDF.FileName, EDF.NS, EDF.Dur and EDF.EDF.SampleRate must be defined.
% 
% Author: (C) 1997-2002 by Alois Schloegl, 15 Jun 2002 #0.85, (Header reworked for 
%         EEGLAB format, Arnaud Delorme and Scott Makeig, 27 Dec 2002)
%
% See also: fopen, SDFREAD, SDFWRITE, SDFCLOSE, SDFSEEK, SDFREWIND, SDFTELL, SDFEOF

% 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 (1999)
% See also: http://www-dpmi.tu-graz.ac.at/~schloegl/matlab/eeg/gdf4/tr_gdf.ps
%
% [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+'); 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,'w');
%           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','w',N); % N number of channels
%            .. do anything, e.g define header
%       EDF=sdfopen(EDF,'w+'); % 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.

function [EDF,H1,h2]=sdfopen(arg1,arg2,arg3,arg4,arg5,arg6)

    INFO='(C) 1997-2002 by Alois Schloegl, 15 Jun 2002 #0.85';
%	[e-mail]
%	Version 0.85

if nargin<2, arg2='r';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.82;

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 (strcmp(arg2,'r') | strcmp(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 GDFOPEN: %s %s\n',MESSAGE,FILENAME);  
        H1=MESSAGE; H2=FILENAME;
        EDF.ErrNo = [32,EDF.ErrNo];
	return;
end;
if (arg2=='r') 
        EDF.FILE.OPEN = 1;
elseif (arg2=='r+') 
        EDF.FILE.OPEN = 2;
end;
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];

%%% 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       ') | 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 IsGDF
        %EDF.T0=[str2num(H1(168+[1:4])) str2num(H1(168+[5 6])) str2num(H1(168+[7 8])) str2num(H1(168+[9 10])) str2num(H1(168+[11 12])) str2num(H1(168+[13:16]))/100 ];
        if 1, % if strcmp(EDF.VERSION(4:8),' 0.12'); 
	% Note: in future versions the date format might change. 
      		EDF.T0(1) = str2num( H1(168 + [ 1:4]));
		EDF.T0(2) = str2num( H1(168 + [ 5 6]));
        	EDF.T0(3) = str2num( H1(168 + [ 7 8]));
        	EDF.T0(4) = str2num( H1(168 + [ 9 10]));
        	EDF.T0(5) = str2num( H1(168 + [11 12]));
        	EDF.T0(6) = str2num( H1(168 + [13:16]))/100;
     	end; 
     
	if str2num(EDF.VERSION(4:8))<0.12
                tmp = setstr(fread(EDF.FILE.FID,8,'uchar')');    % 8 Byte  Length of Header
                EDF.HeadLen = str2num(tmp);    % 8 Byte  Length of Header
        else
                EDF.HeadLen = fread(EDF.FILE.FID,1,'int64');    % 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
        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=[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])) ];
        
        EDF.T0 = zeros(1,6);
        ErrT0=0;
        tmp = str2num( H1(168 + [ 7  8]));
        if ~isempty(tmp), EDF.T0(1) = tmp; else ErrT0 = 1; end;
        tmp = str2num( H1(168 + [ 4  5]));
        if ~isempty(tmp), EDF.T0(2) = tmp; else ErrT0 = 1; end;
        tmp = str2num( H1(168 + [ 1  2]));
        if ~isempty(tmp), EDF.T0(3) = tmp; else ErrT0 = 1; end;
        tmp = str2num( H1(168 + [ 9 10]));
        if ~isempty(tmp), EDF.T0(4) = tmp; else ErrT0 = 1; end;
        tmp = str2num( H1(168 + [12 13]));
        if ~isempty(tmp), EDF.T0(5) = tmp; else ErrT0 = 1; end;
        tmp = str2num( H1(168 + [15 16]));
        if ~isempty(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 = str2num([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=str2num([tmp1,' ']);
                        if isempty(tmp1)
                                ErrT0 = ErrT0 | 1;
                        else
                                EDF.T0(k)  = tmp1;
                        end;
                end;
                if ErrT0
                        EDF.ErrNo = [2,EDF.ErrNo];
                end;
        else
                % Y2K compatibility until year 2084
                if EDF.T0(1) < 85    % for biomedical data recorded in the 1950's and converted to EDF
                        EDF.T0(1) = 2000+EDF.T0(1);
                elseif EDF.T0(1) < 100
                        EDF.T0(1) = 1900+EDF.T0(1);
                %else % already corrected, do not change
                end;
        end;     
        H1(185:256)=setstr(fread(EDF.FILE.FID,256-184,'uchar')');     %
        EDF.HeadLen = str2num(H1(185:192));           % 8 Bytes  Length of Header
        EDF.reserved1=H1(193:136);              % 44 Bytes reserved   
        EDF.NRec    = str2num(H1(237:244));     % 8 Bytes  # of data records