edfopen.m

该程序可以读取欧洲格式.edf和.rec的数据

function [EDF,H1,h2]=edfopen(arg1,arg2,arg3)
% [EDF]=edfopen(FILENAME,'r')
% Opens an EDF File for reading 
% Defines an EDF-struct
%
% EDF.FILE.FID = -1 indicates that file has not been opened properly.
%
%	For use under Octave define the function UPPER.M
% function s=upper(s); s=toupper(s); end;


% [EDF]=edfopen(EDF,'w')
% Opens an EDF File for writing
% EDF is a struct as defined in EDFOPEN(filename,'r')
%
% The following elements have to be defined:
% EDF.FileName     Filename
% EDF.NS           number of channels
% EDF.Dur
% EDF.NRec         -1 if unknown
% EDF.SPR
% EDF.DigMax
% EDF.DigMin
%
% optional can be defined
% !!! ATTENTION, the elements have to have the SAME TYPE and SIZE as defined in the EDF-struct
% EDF.PID
% EDF.RID
% EDF.T0
% EDF.Label
% EDF.Transducer
% EDF.PhysDim
% EDF.PreFilt
% EDF.sPhysMin
% EDF.sPhysMax
%

INFO='(C) 1997-1998 by Alois Schloegl, 13.Nov.1998';
%	[e-mail]
%	12.11.1998

% 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.

if nargin<2 arg2='r';end;
if nargin<3 arg3='ieee-le'; end;

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];

%%%%%%% ============= READ ===========%%%%%%%%%%%%
if arg2=='r'

[fid,MESSAGE]=fopen(FILENAME,'r',arg3);          
EDF.FILE.FID=fid;
if fid<0 
        fprintf(2,'Error EDFOPEN: %s %s\n',MESSAGE,FILENAME);  
        H1=MESSAGE; H2=FILENAME;
	return;
end;
EDF.FILE.OPEN = 1;
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 %%%

H1=setstr(fread(EDF.FILE.FID,256,'char')');     %
%H1=setstr(h1);

EDF.VERSION=H1(1:8);                     % 8 Byte  Versionsnummer 
%if 0 fprintf(2,'LOADEDF: WARNING  Version EDF Format %i',ver); end;
EDF.PID = deblank(H1(9:88));                  % 80 Byte local patient identification
EDF.RID = deblank(H1(89:168));                % 80 Byte local recording identification
%EDF.H.StartDate = H1(169:176);         % 8 Byte		
%EDF.H.StartTime = H1(177:184);         % 8 Byte		
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])) ];

% Y2K compatibility until year 2090
if EDF.VERSION(1)=='0'
        if EDF.T0(1)<91
                EDF.T0(1)=2000+EDF.T0(1);
        else
                EDF.T0(1)=1900+EDF.T0(1);
        end;
else ;
        % in a future version, this is hopefully not needed   
end;

EDF.HeadLen = str2num(H1(185:192));  % 8 Byte  Length of Header
% reserved = H1(193:236);	         % 44 Byte		
EDF.NRec = str2num(H1(237:244));     % 8 Byte  # of data records
EDF.Dur = str2num(H1(245:252));      % 8 Byte  # duration of data record in sec
EDF.NS = str2num(H1(253:256));       % 8 Byte  # of signals

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

h2=zeros(EDF.NS,256);
H2=fread(EDF.FILE.FID,EDF.NS*256,'char');
for k=1:length(H2idx);
        h2(:,idx1(k)+1:idx1(k+1))=reshape(H2(idx2(k)+1:idx2(k+1)),H2idx(k),EDF.NS)';
end;
h2=setstr(h2);

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.sPhysmin   =         h2(:,idx1(4)+1:idx1(5));
EDF.sPhysmax   =         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));

if 0
%fseek(EDF.FILE.FID,256,'bof');
EDF.Label = setstr(fread(EDF.FILE.FID,[16,EDF.NS],'char')');		
EDF.Transducer = setstr(fread(EDF.FILE.FID,[80,EDF.NS],'char')');	
EDF.PhysDim = setstr(fread(EDF.FILE.FID,[8,EDF.NS],'uchar')');
EDF.sPhysmin= setstr(fread(EDF.FILE.FID,[8,EDF.NS],'char')');	
EDF.sPhysmax= setstr(fread(EDF.FILE.FID,[8,EDF.NS],'char')');	
EDF.DigMin = str2num(setstr(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));	%	
EDF.DigMax = str2num(setstr(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));	%	
EDF.PreFilt= setstr(fread(EDF.FILE.FID,[80,EDF.NS],'char')');	%	
tmp = fread(EDF.FILE.FID,[8,EDF.NS],'char')';	%	samples per data record
EDF.SPR = str2num(setstr(tmp));	%	samples per data record
fseek(EDF.FILE.FID,32*EDF.NS,0);
end;

EDF.PhysMin = str2num(EDF.sPhysmin);	
EDF.PhysMax = str2num(EDF.sPhysmax);	


% check validity of DigMin and DigMax
if (length(EDF.DigMin) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Digital Minimum\n');
        EDF.DigMin = -(2^15)*ones(EDF.NS,1);
end
if (length(EDF.DigMax) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Digital Maximum\n');
        EDF.DigMax = (2^15-1)*ones(EDF.NS,1);
end
if (any(EDF.DigMin >= EDF.DigMax))
        fprintf(2,'Warning OPENEDF: Digital Minimum larger than Maximum\n');
end  
% check validity of PhysMin and PhysMax
if (length(EDF.PhysMin) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Physical Minimum\n');
        EDF.PhysMin = EDF.DigMin;
end
if (length(EDF.PhysMax) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Physical Maximum\n');
        EDF.PhysMax = EDF.DigMax;
end
if (any(EDF.PhysMin >= EDF.PhysMax))
        fprintf(2,'Warning OPENEDF: Physical Minimum larger than Maximum\n');
        EDF.PhysMin = EDF.DigMin;
        EDF.PhysMax = EDF.DigMax;
end  

EDF.Cal = (EDF.PhysMax-EDF.PhysMin)./ (EDF.DigMax-EDF.DigMin);
EDF.Off = EDF.PhysMin - EDF.Cal .* EDF.DigMin;
tmp = find(EDF.Cal < 0);
EDF.Cal(tmp) = ones(size(tmp));
EDF.Off(tmp) = zeros(size(tmp));

EDF.Calib=[EDF.Off';(diag(EDF.Cal))];
%EDF.Calib=sparse(diag([1; EDF.Cal]));
%EDF.Calib(1,2:EDF.NS+1)=EDF.Off';

EDF.SampleRate = EDF.SPR / EDF.Dur;

%EDF.FILE.POS = ftell(EDF.FILE.FID);
if EDF.NRec == -1   % unknown record size, determine correct NRec
        fseek(EDF.FILE.FID, 0, 'eof');
        endpos = ftell(EDF.FILE.FID);
        EDF.NRec = floor((endpos - EDF.HeadLen) / (sum(EDF.SPR) * 2));
        fseek(EDF.FILE.FID, EDF.HeadLen, 'bof');
        H1(237:244)=sprintf('%-8i',EDF.NRec);      % write number of records
end; 

EDF.Chan_Select=(EDF.SPR==max(EDF.SPR));
for k=1:EDF.NS;
	if EDF.Chan_Select(k);
	    EDF.ChanTyp(k)='N';
	else
	    EDF.ChanTyp(k)=' ';
	end;         
	if findstr(upper(EDF.Label(k,:)),'ECG')
	    EDF.ChanTyp(k)='C';
	elseif findstr(upper(EDF.Label(k,:)),'EKG')
	    EDF.ChanTyp(k)='C';
	elseif findstr(upper(EDF.Label(k,:)),'EEG')
	    EDF.ChanTyp(k)='E';
	elseif findstr(upper(EDF.Label(k,:)),'EOG')
	    EDF.ChanTyp(k)='O';
	elseif findstr(upper(EDF.Label(k,:)),'EMG')
	    EDF.ChanTyp(k)='M';
	end;
end;

bi=[0;cumsum(EDF.SPR)]; 

idx=[];idx2=[];
for k=1:EDF.NS, 
	idx2=[idx2, (k-1)*max(EDF.SPR)+(1:EDF.SPR(k))];
end;
maxspr=max(EDF.SPR);
idx3=zeros(EDF.NS*maxspr,1);
for k=1:EDF.NS, idx3(maxspr*(k-1)+(1:maxspr))=bi(k)+ceil((1:maxspr)'/maxspr*EDF.SPR(k));end;

EDF.AS.spb = sum(EDF.SPR);	% Samples per Block
EDF.AS.bi=bi;
%EDF.AS.IDX2=idx2;
%EDF.AS.IDX3=idx3;

EDF.AS.Ver=3.09;
EDF.AS.Date=fix(clock);
EDF.AS.Info=INFO;
if exist('OCTAVE_VERSION')
        EDF.AS.Method='edfopen';
else 
        EDF.AS.Method=mfilename;
end;
return;
        
%%%%%%% ============= WRITE ===========%%%%%%%%%%%%        

elseif arg2=='w'
        fprintf(2,'error EDFOPEN: write mode not possible.\n'); 
        H1=[];H2=[];
        return;

end;