edfload.m

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

function [S,count]=edfload(EDF_Filename,arg2,Mode,arg4)
% Loads selected Records of an EDF File (European Data Format for Biosignals)
% [S]=edfload(EDF_Filename,[NoR STARTPOS],Mode,CHAN)
% all elements of CHAN must be non-negative
%         In case CHAN=0 all channels are loaded
%         size(S)=[max(samples_per_chan), number_of_channels]
%
% [S]=edfload(EDF_Filename,[NoR STARTPOS],Mode,ReRefMx)
% the channels are re-referenced with ReRefMx. Usually Mode=0 has to be used.
% any element in ReRefMx must be negative, otherwise arg4 is used as channel selector CHAN
%         the channels are ordered column-wise
%         size(S)=[max(samples_per_chan), size(ReRefMx,2)]
%
% [S]=edfload(EDF_Filename,[NoR STARTPOS],Mode) 
% one block per column (EDF raw form)
%          size(S)=[samples_per_record x number_of_Records]
%
% INPUT:
%   EDF_Filename must be a valid EDF-File 
%   NoR      read next NoR units, units can be seconds or blocks 
%   STARTPOS [otional] goto STARTPOS before reading, units can be seconds or blocks 
%   Mode	0 	autoscaling,      read NoR blocks
%       	1* 	No autoscaling,   read NoR blocks
%        	2 	autoscaling,      read NoR seconds
%       	3 	No autoscaling,   read NoR seconds
%
% OUTPUT:
%   S           data
%
% This program is the Matlab4.x compatible version of the EDF3 toolbox
% 
   
%	Version 3.07
%	05.11.1998
%	Copyright (c) 1997-98 by Alois Schloegl
%	[e-mail]	
                      
% 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.

fid=fopen(EDF_Filename,'r','ieee-le');
if fid<0 
        fprintf(2,['Error LOADEDF: File ' FILENAME ' not found\n']);  
        return;
end;

H1=setstr(fread(fid,256,'char')');      %
VER=H1(1:8);                            %       8 Byte          Versionsnummer 
% if 0 fprintf(2,'LOADEDF: WARNING  Version EDF Format %i',ver); end;
pid = H1(9:88);                 %       80 Byte         local patient identification
rid = H1(89:168);               %       80 Byte         local recording identification

%startdate = H1(169:176);        %       8 Byte          
%starttime = H1(177:184);        %       8 Byte          
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 VER(1)=='0'
        if T0(1)<91
                T0(1)=2000+T0(1);
        else
                T0(1)=1900+T0(1);
        end;
else ;
        % in a future version, this is hopefully not needed   
end;

headlen = str2num(H1(185:192)); %       8 Byte          Length of Header
% reserved = H1(193:236);       %       44 Byte         
nrec = str2num(H1(237:244));    %       8 Byte          # of data records
dur = str2num(H1(245:252));     %       8 Byte          # duration of data record in sec
ns = str2num(H1(253:256));      %       8 Byte          # of signals

Label = setstr(fread(fid,[16,ns],'char')');     %       
transducer = setstr(fread(fid,[80,ns],'char')');        %       
PhysDim = setstr(fread(fid,[8,ns],'char')');    %       

PHYSMIN=setstr(fread(fid,[8,ns],'char')');      %       
PHYSMAX=setstr(fread(fid,[8,ns],'char')');      %       

physmin=zeros(ns,1);
physmax=zeros(ns,1);
for k=1:ns,
        tmp=str2num(PHYSMIN(k,:));
        if isempty(tmp) physmin(k)=0;   %       
        else physmin(k)=tmp; end;       %       

        tmp=str2num(PHYSMAX(k,:));
        if isempty(tmp) physmax(k)=-1;  %       
        else physmax(k)=tmp; end;       %       
end;

digmin = str2num(setstr(fread(fid,[8,ns],'char')'));    %       
digmax = str2num(setstr(fread(fid,[8,ns],'char')'));    %       
prefilt= setstr(fread(fid,[80,ns],'char')');    %       
spr = fread(fid,[8,ns],'char')';        %       samples per data record
SPR = str2num(setstr(spr));     %       samples per data record

%fread(fid,[32,ns],'char')';    %       reserved
fseek(fid,32*ns,0);

CAL = (physmax-physmin)./(digmax-digmin);
OFF = physmin-CAL.*digmin;
%OFF = physmax-CAL.*digmax;
tmp = find(CAL<0);
CAL(tmp) = ones(size(tmp));
OFF(tmp) = ones(size(tmp));

SAMPLERATE = SPR/dur;

spb = sum(SPR); % Samples per Block
bi=[0; cumsum(SPR)];

if nargin<4 
        chan=1:ns; 
else
        if arg4==0 
                chan=1:ns; 
                %ReRefMx=eye(ns);
        elseif any(arg4<0)
                ReRefMx=arg4;
                [chan]=find(any(arg4,2))';
        else
                chan=reshape(arg4,1,prod(size(arg4)));
                ReRefMx=sparse(arg4,arg4,1,ns,ns);
        end;
        if any(SPR(chan)~=SPR(chan(1))) fprintf(2,'Warning EDFLOAD: channels do not have the same sampling rate\n');end;
end;
if nargin<3 Mode=1; end;
if nargin<2 arg2=1;end;

% Position file pointer and calculate number of Records
if fix(Mode./2)==1
        if ~all(~rem(arg2,dur))
                fprintf(2,'Warning EDFLOAD: arg2 does not fit to the EDF blocklength of %i s\n',dur);
        end;
        arg2=floor(arg2/dur);
end;

if length(arg2)>1
        [status]=fseek(fid,(headlen+arg2(2)*spb*2),'bof');
end;
Records=arg2(1);

[s, count]=fread(fid,[spb, Records],'int16');
count=floor(count/spb);

if count~=Records fprintf(2,'Warning EDFLOAD: only % i blocks were read instead of %i\n',count,Records); end;

bi=[0;cumsum(SPR)]; 
if ~rem(Mode,2)         % Autocalib 
        for k=chan,
                s(bi(k)+1:bi(k+1),:)=s(bi(k)+1:bi(k+1),:)*CAL(k)+OFF(k);
        end;
end;

if nargin==4
        maxspr=max(SPR(chan));
        S=zeros(maxspr*Records,length(chan)); 

        if ~all(SPR(chan)==SPR(chan(1)))
                for k=1:length(chan),
                        K=chan(k);
                        tmp=reshape(s(bi(K)+1:bi(K+1),:),SPR(K)*Records,1);
                        S(:,k)=reshape(tmp(:,ones(1,maxspr/SPR(K)))',maxspr*Records,1);
                end;
        else
                for k=1:length(chan),
                        K=chan(k);
                        S(:,k)=reshape(s(bi(K)+1:bi(K+1),:),maxspr*Records,1);
                end;
        end;
        if any(arg4<0)
                S=S*ReRefMx(chan,:);
        end;
else
        S=s;
end;