sdfopen.m

matlab数字信号处理工具箱

                % 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 = str2double(H1(185:192));           % 8 Bytes  Length of Header
        EDF.reserved1=H1(193:236);              % 44 Bytes reserved   
        EDF.NRec    = str2double(H1(237:244));     % 8 Bytes  # of data records
        EDF.Dur     = str2double(H1(245:252));     % 8 Bytes  # duration of data record in sec
        EDF.NS      = str2double(H1(253:256));     % 4 Bytes  # of signals
	EDF.AS.H1   = H1;	                     % for debugging the EDF Header
end;

if strcmp(EDF.reserved1(1:4),'EDF+'),	% EDF+ specific header information 
	[EDF.Patient.Id,   tmp] = strtok(EDF.PID,' ');
	[EDF.Patient.Sex,  tmp] = strtok(tmp,' ');
	[EDF.Patient.Birthday, tmp] = strtok(tmp,' ');
	[EDF.Patient.Name, tmp] = strtok(tmp,' ');

	[chk, tmp] = strtok(EDF.RID,' ');
	if ~strcmp(chk,'Startdate')
		fprinf(EDF.FILE.stderr,'Warning SDFOPEN: EDF+ header is corrupted.\n');
	end;
	[EDF.Date2, tmp] = strtok(tmp,' ');
	[EDF.ID.Investigation, tmp] = strtok(tmp,' ');
	[EDF.ID.Investigator,  tmp] = strtok(tmp,' ');
	[EDF.ID.Equiment, tmp] = strtok(tmp,' ');
end;

if any(size(EDF.NS)~=1) %%%%% 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;
% Octave assumes EDF.NS is a matrix instead of a scalare. Therefore, we need
% Otherwise, eye(EDF.NS) will be executed as eye(size(EDF.NS)).
EDF.NS = EDF.NS(1);     

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 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;
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 ~strcmp(EDF.VERSION(1:3),'GDF'),
        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    = str2double(h2(:,idx1(4)+1:idx1(5)));
        EDF.PhysMax    = str2double(h2(:,idx1(5)+1:idx1(6)));
        EDF.DigMin     = str2double(h2(:,idx1(6)+1:idx1(7)));
        EDF.DigMax     = str2double(h2(:,idx1(7)+1:idx1(8)));
        EDF.PreFilt    =         h2(:,idx1(8)+1:idx1(9));
        EDF.SPR        = str2double(h2(:,idx1(9)+1:idx1(10)));
        %EDF.reserved  =       h2(:,idx1(10)+1:idx1(11));
	if ~all(abs(EDF.VERSION)==[255,abs('BIOSEMI')]),
		EDF.GDFTYP     = 3*ones(1,EDF.NS);	%	datatype
	else
		EDF.GDFTYP     = (255+24)*ones(1,EDF.NS);	%	datatype
	end;
	
        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');	
	tmp            =         fread(EDF.FILE.FID,[2*EDF.NS,1],'int32');	
        EDF.DigMin     = tmp((1:EDF.NS)*2-1);
        tmp            =         fread(EDF.FILE.FID,[2*EDF.NS,1],'int32');	
        EDF.DigMax     = tmp((1:EDF.NS)*2-1);
        
        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.Filter.LowPass = repmat(nan,1,EDF.NS);
EDF.Filter.HighPass = repmat(nan,1,EDF.NS);
for k=1:EDF.NS,
	tmp = EDF.PreFilt(k,:);

	ixh=findstr(tmp,'HP');
	ixl=findstr(tmp,'LP');
	ix=findstr(lower(tmp),'hz');
	%tmp(tmp==':')=' ';
try;
	if any(tmp==';')
		[tok1,tmp] = strtok(tmp,';');
		[tok2,tmp] = strtok(tmp,';');
	else
		[tok1,tmp] = strtok(tmp,' ');
		[tok2,tmp] = strtok(tmp,' ');
	end;
	[T1, F1 ] = strtok(tok1,': ');
	[T2, F2 ] = strtok(tok2,': ');

	[F1 ] = strtok(F1,': ');
	[F2 ] = strtok(F2,': ');
	
	if strcmp(F1,'DC'), F1='0'; end;
	if strcmp(F2,'DC'), F2='0'; end;
	
	tmp = findstr(lower(F1),'hz');
	if ~isempty(tmp), F1=F1(1:tmp-1);end;
	tmp = findstr(lower(F2),'hz');
	if ~isempty(tmp), F2=F2(1:tmp-1);end;

	if strcmp(T1,'LP'), 
		EDF.Filter.LowPass(k) =str2double(F1);
	elseif strcmp(T1,'HP'), 
		EDF.Filter.HighPass(k)=str2double(F1);
	end;
	if strcmp(T2,'LP'), 
		EDF.Filter.LowPass(k) =str2double(F2);
	elseif strcmp(T2,'HP'), 
		EDF.Filter.HighPass(k)=str2double(F2);
	end;
catch
	fprintf(1,'Cannot interpret: %s\n',EDF.PreFilt(k,:));
end;
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.EDF.SampleRate = EDF.SPR / EDF.Dur;
EDF.AS.MAXSPR=EDF.SPR(1);
for k=2:EDF.NS,
        EDF.AS.MAXSPR = lcm(EDF.AS.MAXSPR,EDF.SPR(k));
end;
EDF.SampleRate = EDF.AS.MAXSPR/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.AS.EVENTTABLEPOS = -1;


% EDF+: 
tmp = strmatch('EDF Annotations',EDF.Label);
if strcmp(EDF.TYPE,'EDF') & (length(tmp)==1),
        EDF.EDF.Annotations = tmp;
        EDF.Cal(EDF.EDF.Annotations) = 1;
        EDF.Off(EDF.EDF.Annotations) = 0;
        
        fseek(EDF.FILE.FID,EDF.HeadLen+EDF.AS.bi(EDF.EDF.Annotations)*2,'bof');
        t = fread(EDF.FILE.FID,EDF.SPR(EDF.EDF.Annotations),'uchar',EDF.AS.bpb-EDF.SPR(EDF.EDF.Annotations)*2);
        t = char(t)';
        lt = length(t);
        EDF.EDF.ANNONS = t;
        N = 0; 
        ix = 1; 
	t = [t,' '];
        while ix < length(t),
                while (ix<=lt) & (t(ix)==0), ix = ix+1; end;
                ix1 = ix; 
                while (ix<=lt) & (t(ix)~=0), ix = ix+1; end;
                ix2 = ix; 
                if (ix < lt),
                        v = t(ix1:ix2-1);
                        [s1,v]=strtok(v,20);
                        s1(s1==21) = 0;
                        tmp=str2double(char(s1));
                        
                        [s2,v]=strtok(v,20);
                        [s3,v]=strtok(v,20);
                        
                        N = N+1;
                        EDF.EVENT.POS(N,1) = tmp(1);
                        if length(tmp)>1,
                                EDF.EVENT.DUR(N,1) = tmp(2);
                        end;
                        EDF.EVENT.TeegType{N,1} = char(s2);
                        EDF.EVENT.TeegDesc{N,1} = char(s3);
                end;
        end;
        EDF.EVENT.TYP(1:N,1) = 0;
        EDF.EVENT.N = N; 

elseif strcmp(EDF.TYPE,'EDF') & (length(EDF.FILE.Name)==8) & any(lower(EDF.FILE.Name(1))=='bchmnpsu') 
        if strcmp(lower(EDF.FILE.Name([3,6:8])),'001a'),
                % load scoring of ADB database if available 

                fid2 = fopen(fullfile(EDF.FILE.Path, [EDF.FILE.Name(1:7),'.txt']),'r');
                if fid2<0,
                        fid2 = fopen(fullfile(EDF.FILE.Path,[lower(EDF.FILE.Name(1:7)),'.txt']),'r');
                end
                if fid2<0,
                        fid2 = fopen(fullfile(EDF.FILE.Path,[EDF.FILE.Name(1:7),'.TXT']),'r');
                end
                if fid2>0,
                        tmp = fread(fid2,inf,'char');
                        fclose(fid2);
                        [ma,status] = str2double(char(tmp'));
                        if ~any(isnan(status(:))),
                                %%% TODO: include ADB2EVENT here

				% code from MAK2BIN.M (C) 1998-2004 A. Schl鰃l 
				ERG = zeros(size(ma));

				%%%% one artifact %%%%
				for k=0:9,
				        if exist('OCTAVE_VERSION')==5
				                ERG = ERG+(ma==k)*2^k;
				        else
				                ERG(ma==k) = 2^k;
				        end;
				end;

				%%%% more than one artifact %%%%
				[i,j] = find(ma>9);
				L='123456789';
				for k=1:length(i),
				        b=int2str(ma(i(k),j(k)));
				        erg=0;
				        for l=1:9,
				                if any(b==L(l)), erg=erg+2^l; end;
				        end;        
				        ERG(i(k),j(k)) = erg;
				end;
				
				N   = 0;
				POS = [];
				TYP = [];
				DUR = [];
				CHN = [];
				cc  = zeros(1,10);
				for k = 1:9,
				        for c = 1:7;%size(ERG,2),
				                tmp = [0;~~(bitand(ERG(:,c),2^k));0];
				 
				                cc(k+1) = cc(k+1) + sum(tmp);
				                pos = find(diff(tmp)>0);
				                pos2 = find(diff(tmp)<0);
				                n   = length(pos);
				                
				                POS = [POS; pos(:)];
				                TYP = [TYP; repmat(k,n,1)];
				                CHN = [CHN; repmat(c,n,1)];
				                DUR = [DUR; pos2(:)-pos(:)];
				                N   = N + n;
				        end;
				end;