sdfopen.m

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

        EDF.Dur     = str2num(H1(245:252));     % 8 Bytes  # duration of data record in sec
        EDF.NS      = str2num(H1(253:256));     % 4 Bytes  # of signals
	EDF.AS.H1   = H1;	                     % for debugging the EDF Header
end;

if isempty(EDF.NS) %%%%% not EDF because filled out with ASCII(0) - should be spaces
        fprintf(EDF.FILE.stderr, 'Warning SDFOPEN: invalid NS-value in header of %s\n',EDF.FileName);
        EDF.ErrNo=[1040,EDF.ErrNo];
        EDF.NS=1;
end;

if isempty(EDF.HeadLen) %%%%% not EDF because filled out with ASCII(0) - should be spaces
        EDF.ErrNo=[1056,EDF.ErrNo];
        EDF.HeadLen=256*(1+EDF.NS);
end;


if any(~isempty(EDF.reserved1)) %%%%% not EDF because filled out with ASCII(0) - should be spaces
        fprintf(2, 'Warning SDFOPEN: 44bytes-reserved-field of fixed header is not empty in %s\n',EDF.FILE.Name);
        %EDF.ErrNo=[1057,EDF.ErrNo];
end;


if isempty(EDF.NRec) %%%%% not EDF because filled out with ASCII(0) - should be spaces
        EDF.ErrNo=[1027,EDF.ErrNo];
        EDF.NRec = -1;
end;

if isempty(EDF.Dur) %%%%% not EDF because filled out with ASCII(0) - should be spaces
        EDF.ErrNo=[1088,EDF.ErrNo];
        EDF.Dur=30;
else 
	if ~IsGDF,
		if (EDF.Dur>1) & (EDF.Dur~=fix(EDF.Dur)), % Blockduration must be in seconds or subsecond range 
		        EDF.ErrNo=[1088,EDF.ErrNo];
		end;
	end;	
end;

if  any(EDF.T0>[2084 12 31 24 59 59]) | any(EDF.T0<[1985 1 1 0 0 0])
        EDF.ErrNo = [4, EDF.ErrNo];
end;

%%% Read variable Header %%%
if ~IsGDF
        idx1=cumsum([0 H2idx]);
        idx2=EDF.NS*idx1;

        h2=zeros(EDF.NS,256);
        [H2,count]=fread(EDF.FILE.FID,EDF.NS*256,'uchar');
        if count < EDF.NS*256 
	        EDF.ErrNo=[8,EDF.ErrNo];
                return; 
        end;
                
        %tmp=find((H2<32) | (H2>126)); % would confirm 
        tmp = find((H2<32) | ((H2>126) & (H2~=255) & (H2~=181)& (H2~=230))); 
        if ~isempty(tmp) %%%%% not EDF because filled out with ASCII(0) - should be spaces
                H2(tmp) = 32; 
	        EDF.ErrNo = [1026,EDF.ErrNo];
        end;
        
        for k=1:length(H2idx);
                %disp([k size(H2) idx2(k) idx2(k+1) H2idx(k)]);
                h2(:,idx1(k)+1:idx1(k+1))=reshape(H2(idx2(k)+1:idx2(k+1)),H2idx(k),EDF.NS)';
        end;
        %size(h2),
        h2=setstr(h2);
        %(h2(:,idx1(9)+1:idx1(10))),
        %abs(h2(:,idx1(9)+1:idx1(10))),
        
        EDF.Label      =         h2(:,idx1(1)+1:idx1(2));
        EDF.Transducer =         h2(:,idx1(2)+1:idx1(3));
        EDF.PhysDim    =         h2(:,idx1(3)+1:idx1(4));
        EDF.PhysMin    = str2num(h2(:,idx1(4)+1:idx1(5)));
        EDF.PhysMax    = str2num(h2(:,idx1(5)+1:idx1(6)));
        EDF.DigMin     = str2num(h2(:,idx1(6)+1:idx1(7)));
        EDF.DigMax     = str2num(h2(:,idx1(7)+1:idx1(8)));
        EDF.PreFilt    =         h2(:,idx1(8)+1:idx1(9));
        EDF.SPR        = str2num(h2(:,idx1(9)+1:idx1(10)));
        %EDF.reserved  =       h2(:,idx1(10)+1:idx1(11));
        EDF.GDFTYP     = 3*ones(1,EDF.NS);	%	datatype

        if isempty(EDF.SPR), 
                fprintf(EDF.FILE.stderr, 'Warning SDFOPEN: invalid SPR-value in header of %s\n',EDF.FileName);
                EDF.SPR=ones(EDF.NS,1);
	        EDF.ErrNo=[1028,EDF.ErrNo];
        end;
else
        fseek(EDF.FILE.FID,256,'bof');
        EDF.Label      =  setstr(fread(EDF.FILE.FID,[16,EDF.NS],'uchar')');		
        EDF.Transducer =  setstr(fread(EDF.FILE.FID,[80,EDF.NS],'uchar')');	
        EDF.PhysDim    =  setstr(fread(EDF.FILE.FID,[ 8,EDF.NS],'uchar')');
%       EDF.AS.GDF.TEXT = EDF.GDFTYP.TEXT;
        EDF.PhysMin    =         fread(EDF.FILE.FID,[EDF.NS,1],'float64');	
        EDF.PhysMax    =         fread(EDF.FILE.FID,[EDF.NS,1],'float64');	
        EDF.DigMin     =         fread(EDF.FILE.FID,[EDF.NS,1],'int64');	
        EDF.DigMax     =         fread(EDF.FILE.FID,[EDF.NS,1],'int64');	
        
        EDF.PreFilt    =  setstr(fread(EDF.FILE.FID,[80,EDF.NS],'uchar')');	%	
        EDF.SPR        =         fread(EDF.FILE.FID,[ 1,EDF.NS],'uint32')';	%	samples per data record
        EDF.GDFTYP     =         fread(EDF.FILE.FID,[ 1,EDF.NS],'uint32');	%	datatype
        %                        fread(EDF.FILE.FID,[32,EDF.NS],'uchar')';	%	datatype
end;

%		EDF=gdfcheck(EDF,1);
if any(EDF.PhysMax==EDF.PhysMin), EDF.ErrNo=[1029,EDF.ErrNo]; end;	
if any(EDF.DigMax ==EDF.DigMin ), EDF.ErrNo=[1030,EDF.ErrNo]; end;	
EDF.Cal = (EDF.PhysMax-EDF.PhysMin)./(EDF.DigMax-EDF.DigMin);
EDF.Off = EDF.PhysMin - EDF.Cal .* EDF.DigMin;
EDF.Calib=[EDF.Off';(diag(EDF.Cal))];
EDF.SampleRate = EDF.SPR / EDF.Dur;

EDF.AS.spb = sum(EDF.SPR);	% Samples per Block
EDF.AS.bi = [0;cumsum(EDF.SPR)]; 
EDF.AS.BPR  = ceil(EDF.SPR.*GDFTYP_BYTE(EDF.GDFTYP+1)'); 
EDF.AS.SAMECHANTYP = all(EDF.AS.BPR == (EDF.SPR.*GDFTYP_BYTE(EDF.GDFTYP+1)')) & all(EDF.GDFTYP(:)==EDF.GDFTYP(1));
EDF.AS.GDFbi = [0;cumsum(ceil(EDF.SPR.*GDFTYP_BYTE(EDF.GDFTYP+1)'))]; 
EDF.AS.bpb = sum(ceil(EDF.SPR.*GDFTYP_BYTE(EDF.GDFTYP+1)'));	% Bytes per Block
EDF.AS.startrec = 0;
EDF.AS.numrec = 0;
EDF.FILE.POS = 0;

status = fseek(EDF.FILE.FID, 0, 'eof');
EDF.AS.endpos = ftell(EDF.FILE.FID);
fseek(EDF.FILE.FID, EDF.HeadLen, 'bof');

%[status EDF.AS.endpos EDF.HeadLen EDF.AS.bpb EDF.NRec]
if EDF.NRec == -1   % unknown record size, determine correct NRec
        EDF.NRec = floor((EDF.AS.endpos - EDF.HeadLen) / EDF.AS.bpb);
elseif  EDF.NRec ~= ((EDF.AS.endpos - EDF.HeadLen) / EDF.AS.bpb);
        EDF.ErrNo=[16,EDF.ErrNo];
%        EDF.NRec = floor((EDF.AS.endpos - EDF.HeadLen) / EDF.AS.bpb);
end; 

% if Channelselect, ReReferenzing and Resampling
% Overflowcheck, Adaptive FIR
% Layer 4 

%if nargin<3 %%%%%          Channel Selection 
if nargin <3 %else
        arg3=0;
end;

EDF.SIE.ChanSelect = 1:EDF.NS;
EDF.SIE.InChanSelect = 1:EDF.NS;

EDF.SIE.RR=1;
EDF.SIE.RS=0; %exist('arg5')==1; if EDF.SIE.RS, EDF.SIE.RS==(arg5>0); end;
EDF.SIE.TH=0; %(exist('arg6')==1);
EDF.SIE.RAW=0;
EDF.SIE.REGC=0;
EDF.SIE.TECG=0;
EDF.SIE.AFIR=0;
EDF.SIE.FILT=0;
EDF.SIE.TimeUnits_Seconds=1;
EDF.AS.MAXSPR=max(EDF.SPR(EDF.SIE.ChanSelect)); % Layer 3 defines EDF.AS.MAXSPR in GDFREAD

EDF.SIE.ReRefMx=eye(EDF.NS);
EDF.SIE.REG=eye(EDF.NS);
EDF.SIE.ReRefMx = EDF.SIE.ReRefMx(:,EDF.SIE.ChanSelect);
EDF.Calib=EDF.Calib*EDF.SIE.REG*EDF.SIE.ReRefMx; 

%if nargin>2 %%%%%          Channel Selection 
        EDF.SIE.REG=eye(EDF.NS);
        if arg3==0
                EDF.SIE.ChanSelect = 1:EDF.NS;
                EDF.SIE.InChanSelect = 1:EDF.NS;
                EDF.SIE.ReRefMx = eye(EDF.NS);
        else
                [nr,nc]=size(arg3);
                if all([nr,nc]>1), % Re-referencing
                        EDF.SIE.ReRefMx = [[arg3 zeros(size(arg3,1),EDF.NS-size(arg3,2))]; zeros(EDF.NS-size(arg3,1),EDF.NS)];
                        EDF.SIE.InChanSelect = find(any(EDF.SIE.ReRefMx'));
                        EDF.SIE.ChanSelect = find(any(EDF.SIE.ReRefMx));
                        if nargin>3
                        %        fprintf(EDF.FILE.stderr,'Error SDFOPEN: Rereferenzing does not work correctly with %s (more than 3 input arguments)\n',arg4);
                        end;
                else
                        EDF.SIE.ChanSelect = arg3; %full(sparse(1,arg3(:),1,1,EDF.NS));
                        EDF.SIE.InChanSelect = arg3;
                        EDF.SIE.ReRefMx = eye(EDF.NS);
                end;
        end;
        
        if (exist('sedfchk')==2),  
                EDF=sedfchk(EDF); % corrects incorrect Header information. 
        end;
        
	%EDF.SIE.CS=1;
	EDF.SIE.RR=1;
        EDF.SIE.RS=exist('arg5')==1; if EDF.SIE.RS, EDF.SIE.RS==(arg5>0); end;
        EDF.SIE.TH=0; %(exist('arg6')==1);
        EDF.SIE.RAW=0;
        EDF.SIE.REGC=0;
        EDF.SIE.TECG=0;
        EDF.SIE.AFIR=0;
        EDF.SIE.FILT=0;
        EDF.AS.MAXSPR=max(EDF.SPR(EDF.SIE.ChanSelect)); % Layer 3 defines EDF.AS.MAXSPR in GDFREAD
        
        EDF.SIE.REG=eye(EDF.NS);

        %elseif nargin>3   %%%%% RE-REFERENCING
        if nargin>3   
                if ~isempty(findstr(upper(arg4),'SIESTA'))
                        EDF.SIE.ReRefMx = eye(EDF.NS);% speye(EDF.NS);
                        EDF.SIE.ReRefMx(7,1:6)=[1 1 1 -1 -1 -1]/2;
                        EDF.SIE.TH=1;
	                % calculate (In)ChanSelect based on ReRefMatrix
    		        EDF.SIE.InChanSelect = find(any(EDF.SIE.ReRefMx'));
            		EDF.SIE.ChanSelect = find(any(EDF.SIE.ReRefMx));
                end;
                if ~isempty(findstr(upper(arg4),'EOG'))
                        tmp=findstr(upper(arg4),'EOG');
                        tmp=lower(strtok(arg4(tmp:length(arg4)),' +'));
                        EDF.SIE.ReRefMx = sparse(EDF.NS,EDF.NS);% speye(EDF.NS);
                        if any(tmp=='h'); EDF.SIE.ReRefMx(8:9,1)=[ 1 -1 ]';  end;
                        if any(tmp=='v'); EDF.SIE.ReRefMx(1:9,2)=[ 1 0 0 1 0 0 1 -1 -1]';end;
                        if any(tmp=='r'); EDF.SIE.ReRefMx(3:9,3)=[ 1 0 0 1 2 -2 -2]'/2;  end;
                        %EDF.SIE.TH=1;
	                % calculate (In)ChanSelect based on ReRefMatrix
    		        EDF.SIE.InChanSelect = find(any(EDF.SIE.ReRefMx'));
            		EDF.SIE.ChanSelect = find(any(EDF.SIE.ReRefMx));
                end;
                
                if ~isempty(findstr(upper(arg4),'UCAL'))
                        %EDF.SIE.ReRefMx = speye(EDF.NS); % OVERRIDES 'SIESTA' and 'REGRESS_ECG'
                        EDF.SIE.RR = 0;
                end;
                if ~isempty(findstr(upper(arg4),'RAW'))
                        EDF.SIE.RAW = 1;
                end;
                if ~isempty(findstr(upper(arg4),'OVERFLOW'))
                        EDF.SIE.TH = 1;
                end;
                if ~isempty(findstr(upper(arg4),'FailingElectrodeDetector'))
                        EDF.SIE.FED = 1;
			EDF.SIE.TH = 2; 
                end;
		
                if ~isempty(findstr(upper(arg4),'ECG')), % identify ECG channel for some ECG artifact processing method   
                        if ~isempty(findstr(upper(arg4),'SIESTA'))
                                channel1=12;
                                %channel2=1:9;
                                M=zeros(EDF.NS,1);M(channel1)=1;
                        end
                        if isfield(EDF,'ChanTyp')
                                M=upper(char(EDF.ChanTyp(channel1)))=='C';
                                channel1=find(M);
                                M=(upper(char(EDF.ChanTyp))=='E' | upper(char(EDF.ChanTyp))=='O' );
                                %channel2=find(M);
                        else
                                channel1=12;
                                %channel2=1:9;
                                M=zeros(EDF.NS,1);M(channel1)=1;
                        end;
                        
                end;
                if ~isempty(findstr(upper(arg4),'RECG'))
                        EDF.SIE.REGC = 1;
                        if all(EDF.SIE.InChanSelect~=channel1)
                                EDF.SIE.InChanSelect=[EDF.SIE.InChanSelect channel1];        
                        end;
                end;
                
                if ~isempty(findstr(upper(arg4),'TECG'))
                        fprintf(EDF.FILE.stderr,'SDFOPEN: status of TECG Mode: alpha test passed\n');    
                        %%%% TECG - ToDo
                        % - optimize window
                        if exist([lower(EDF.FILE.Name) 'ECG.mat']);
                                if exist('OCTAVE_VERSION')==5
                                        load(file_in_loadpath([lower(EDF.FILE.Name) 'ECG.mat']));
                                else
                                        load([lower(EDF.FILE.Name) 'ECG.mat']);
                                end;
                                if isstruct(QRS)
                                        EDF.SIE.TECG = 1;
					%%%  EDF.AS.MAXSPR=size(QRS.Templates,1)/3;
                                        EDF.AS.MAXSPR=max(EDF.AS.MAXSPR,EDF.SPR(channel1));
				else
                                        fprintf(EDF.FILE.stderr,'WARNING SDFOPEN: %s invalid for TECG\n',[ lower(EDF.FILE.Name) 'ECG.mat']);
                                end;
                        else
                                fprintf(EDF.FILE.stderr,'WARNING SDFOPEN: %s not found\t (needed for Mode=TECG)\n',[lower(EDF.FILE.Name) 'ECG.mat']);
                        end;
                end;
                
                if ~isempty(findstr(upper(arg4),'AFIR')) 
                        % Implements Adaptive FIR filtering for ECG removal in EDF/GDF-tb.
                        % based on the Algorithm of Mikko Koivuluoma <[e-mail]>
                                
                        % channel2 determines which channels should be corrected with AFIR 
                        
                        %EDF = sdf_afir(EDF,12,channel2);
                        channel2=EDF.SIE.ChanSelect;	
                        fprintf(EDF.FILE.stderr,'Warning SDFOPEN: option AFIR still buggy\n');    
                        if isempty(find(channel2))
                                EDF.SIE.AFIR=0;
                        else
                                EDF.SIE.AFIR=1;
                                EDF.AFIR.channel2=channel2;
                                
                                EDF.AFIR.alfa=0.01;
                                EDF.AFIR.gamma=1e-32;
                                
                                EDF.AFIR.delay  = ceil(0.05*EDF.AS.MAXSPR/EDF.Dur); 
                                EDF.AFIR.nord = EDF.AFIR.delay+EDF.AS.MAXSPR/EDF.Dur; 
                                
                                EDF.AFIR.nC = length(EDF.AFIR.channel2);
                                EDF.AFIR.w = zeros(EDF.AFIR.nC, max(EDF.AFIR.nord));
                                EDF.AFIR.x = zeros(1, EDF.AFIR.nord);
                                EDF.AFIR.d = zeros(EDF.AFIR.delay, EDF.AFIR.nC);
                                
                                channel1=12;
                                
                                if isfield(EDF,'ChanTyp')
                                        if upper(char(EDF.ChanTyp(channel1)))=='C';
                                                EDF.AFIR.channel1 = channel1;
                                        else
                                                EDF.AFIR.channel1 = find(EDF.ChanTyp=='C');
                                                fprintf(EDF.FILE.stderr,'Warning %s: #%i is not an ECG channel, %i used instead\n' ,filename,channel1,EDF.AFIR.channel1);