gdfread.m

matlab下实现

function [S,count]=gdfread(EDF,chan,Mode,NoR,StartPos)
% Loads selected Records of an EDF File (European Data Format for Biosignals)
% [S]=gdfread(EDF_Struct, CHAN, Mode, NoR, StartPos)
%
% Functions of EDFREAD: 
% - Channel selection,           CHAN is list of selected channels, CHAN=0 [default] selects all channels
% - Segment selection            NoR reads NoR units
%                                StartPos sets the begin of reading, if not given reading starts at current position 
%                                      units can be blocks or seconds depending on Mode
% - {R}ereferencing                Mode = 'R' (8), CHAN is the referencing Matrix 
% - Ra{W} vs. Channel order        Mode = 'W' (4), S is Raw EDF format with one block per column
%                                                  size(S)=[samples_per_record x number_of_Records]
%                                                  channel selection and Mode='S' not possible, 
% - {A}utoscaling (Calibration)      Mode = 'A' (2), defines whether Scaling should be applied or Not [default]. 
% - units {S}econds or blocks      Mode = 'S' (1), defines whether units of NoR & StartPos are seconds or blocks 
%                                                  'S' not recommended
%
% Mode can be string e.g. 'RAS' or integer e.g. 8+2+1; both are the same.
% default Mode is '' (0+0+0+0), 'RW' not possible
%
% OUTPUT:
%   S           data
%   

%	Copyright (c) 1997-99 by Alois Schloegl
%	[e-mail]	
%	Version 0.43
%	11. Feb. 1999

% 04 Feb 1999 V0.42 speed improved for reading EDF files
% 11 Feb 1999 V0.43 bug fixed, count for EDF files corrected

% 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

OptiMEM=1;             % if you get low of memory (your CPU starts swapping), set this to one
OptiSPEED=~OptiMEM;

GDF=strcmp(EDF.VERSION(1:3),'GDF');

if nargin<1
        fprintf(2,'error EDFREAD: missing input arguments\n');
elseif nargin<2
        chan=0;
elseif nargin<3
        Mode=0;
elseif nargin<4 
        NoR=EDF.NRec;
elseif nargin<5
        StartPos=0;
end;

if chan==0        
        chan=1:EDF.NS;
end;

if ischar(Mode)         %%% not OCTAVE-compatible
        arg3=upper(Mode);
        Mode = any(arg3=='R')*8 + any(arg3=='W')*4 + (any(arg3=='A') & ~any(arg3=='N'))*2 + any(arg3=='S'); 
end;

Mode_RAW=any(Mode==[4:7 12:15]);
Mode_CAL=any(Mode==[2 3 6 7 10 11 14 15]);
Mode_REF=any(Mode==[8:15]);
Mode_SEC=any(Mode==[1:2:15]);

% Position file pointer and calculate number of Records
if Mode_SEC
        if ~all(~rem([NoR StartPos],EDF.Dur))
                fprintf(2,'Warning EDFREAD: NoR &/ StartPos do not fit to blocklength of EDF File of %i s.\n',EDF.Dur);
        end;
        StartPos=StartPos/EDF.Dur;
        NoR=NoR/EDF.Dur;
end;
        
if ~(nargin<5) 
        [status,POS]=gdfseek(EDF,StartPos,'bof'); 
end;

if Mode_REF
        ReRefMx=chan;
        chan=find(any(chan,2))';
else
        chan=reshape(chan,1,prod(size(chan)));
        ReRefMx=eye(length(chan));
end;

Mode_CHANSAME = ~all(EDF.SPR(chan)==EDF.SPR(chan(1)));

%if any(EDF.SPR(chan)~=EDF.SPR(chan(1))) fprintf(2,'Warning EDFREAD: channels do not have the same sampling rate\n');end;

bi=[0;cumsum(EDF.SPR)];
Records=NoR;  % -1 if bodged to miss the last line in caseit is corrupted
maxspr=max(EDF.SPR(chan));
count=0;

if GDF
        if Mode_RAW;
                S=zeros(sum(EDF.SPR(chan)),Records);
                bi=[0;cumsum(EDF.SPR(chan))];
        else
                S=zeros(maxspr*Records,1+length(chan));
                S(:,1)=1;
        end;        
        
        l=0;
        while (l<Records) & ~feof(EDF.FILE.FID), 
%               for K=1:length(chan), K=chan(k);
                for k=1:EDF.NS,
                                if EDF.GDFTYP(k)==0
                                       [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'uchar');
                                elseif EDF.GDFTYP(k)==1
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'int8');
                                elseif EDF.GDFTYP(k)==2
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'uint8');
                                elseif EDF.GDFTYP(k)==3
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'int16');
                                elseif EDF.GDFTYP(k)==4
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'uint16');
                                elseif EDF.GDFTYP(k)==5
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'int32');
                                elseif EDF.GDFTYP(k)==6
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'uint32');
                                elseif EDF.GDFTYP(k)==7
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'int64');
                                elseif 0; EDF.GDFTYP(k)==8
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'uint64');
                                elseif EDF.GDFTYP(k)==16
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'float32');
                                elseif EDF.GDFTYP(k)==17
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'float64');
                                elseif 0;EDF.GDFTYP(k)>255 & EDF.GDFTYP(k)< 256+64
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),['bit' int2str(EDF.GDFTYP(k))]);
                                elseif 0;EDF.GDFTYP(k)>511 & EDF.GDFTYP(k)< 511+64
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),['ubit' int2str(EDF.GDFTYP(k))]);
                                elseif EDF.GDFTYP(k)==256
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'bit1');
                                elseif EDF.GDFTYP(k)==512
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'ubit1');
                                elseif EDF.GDFTYP(k)==255+12
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'bit12');
                                elseif EDF.GDFTYP(k)==511+12
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'ubit12');
                                elseif EDF.GDFTYP(k)==255+22
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'bit22');
                                elseif EDF.GDFTYP(k)==511+22
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'ubit22');
                                elseif EDF.GDFTYP(k)==255+24
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'bit24');
                                elseif EDF.GDFTYP(k)==511+24
                                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),'ubit24');
                                else 
                                        fprintf(2,'Error GDFREAD: Invalid GDF channel type in %s at channel %i',EDF.FileName,k);
                                end;
                        if any(chan==k), K=find(chan==k);
                                if Mode_RAW;
                                        S(bi(K)+1:bi(K+1),l)=tmp;
                                else 
                                        if Mode_CHANSAME
                                                tmp=reshape(tmp(:,ones(1,maxspr/EDF.SPR(k)))',maxspr,1);
                                        end;
                                        %disp([size(tmp) size(S) k K l maxspr])
                                        S(l*maxspr+(1:maxspr),K+1)=tmp;
                                end;  
                        end;      
                        count=count+cnt;
                end; 
        l=l+1;
        end;
        count=count/EDF.AS.spb;

        if Mode_RAW 
                if Mode_CAL          % Autocalib 
                        for k=1:EDF.NS,
                                S(bi(k)+1:bi(k+1),:)=S(bi(k)+1:bi(k+1),:)*EDF.Cal(k)+EDF.Off(k);
                        end;
                end;
        else
                tmp = Mode_CAL + Mode_REF*2;
                if tmp==0
                        S=S(:,1+(1:length(chan)));
                elseif tmp==1
                        S=S*EDF.Calib([1 chan+1],chan);
                elseif tmp==2
                        S=S*[zeros(1,size(ReRefMx,2)); ReRefMx(chan,:)];
                elseif tmp==3
                        S=S*(EDF.Calib([1 chan+1],chan)*ReRefMx(chan,:));
                end;
        end;
        
else %%%%%%%%% if ~GDF

if Mode_RAW;
        [S, count]=fread(EDF.FILE.FID,[EDF.AS.spb, Records],'int16');
        %count = floor(count/EDF.AS.spb);
        count = (count/EDF.AS.spb);
        if count<1; fprintf(2,'Warning EDFREAD: only %.1f blocks were read instead  of %.1f\n',count,Records); end;

        if Mode_CAL          % Autocalib 
                for k=chan,
                        S(bi(k)+1:bi(k+1),:)=S(bi(k)+1:bi(k+1),:)*EDF.Cal(k)+EDF.Off(k);
                end;
        end;
else
        S=zeros(maxspr*Records,length(chan));
        if all(EDF.SPR(chan)==EDF.SPR(chan(1)))
                if ~OptiMEM % but OptiSPEED
                        [s, count]=fread(EDF.FILE.FID,[EDF.AS.spb,Records],'int16');
                        count = (count/EDF.AS.spb);
                        for k=1:length(chan), K=chan(k);
                                S(:,k)=reshape(s(bi(K)+1:bi(K+1),:),maxspr*Records,1);
                        end;
                else %OptiMEM % but ~OptiSPEED
                        for l=1:Records;
                                [s, count]=fread(EDF.FILE.FID,[EDF.AS.spb,1],'int16');
                                count = (count/EDF.AS.spb);
                                for k=1:length(chan), K=chan(k);
                                        S((l-1)*maxspr+(1:maxspr),k)=s(bi(K)+1:bi(K+1));
                                end;
                        end;
                end;
        else
                if ~OptiMEM % but OptiSPEED
                        [s, count]=fread(EDF.FILE.FID,[EDF.AS.spb,Records],'int16');
                        count = (count/EDF.AS.spb);
                        for k=1:length(chan), K=chan(k);
                                tmp=reshape(s(bi(K)+1:bi(K+1),:),EDF.SPR(K)*Records,1);
                                S(:,k)=reshape(tmp(:,ones(1,maxspr/EDF.SPR(K)))',maxspr*Records,1);
                        end;
                else %OptiMEM % but ~OptiSPEED
                        for l=1:Records;
                                [s, count]=fread(EDF.FILE.FID,[EDF.AS.spb,1],'int16');
                                count = (count/EDF.AS.spb);
                                for k=1:length(chan), K=chan(k);
                                        tmp=reshape(s(bi(K)+1:bi(K+1)),EDF.SPR(K),1);
                                        S((l-1)*maxspr+(1:maxspr),k)=reshape(tmp(:,ones(1,maxspr/EDF.SPR(K)))',maxspr,1);
                                end;
                        end;
                end;
        end;
        
        tmp = Mode_CAL + Mode_REF*2;
        if tmp==0
        elseif tmp==1
                S=[ones(size(S,1),1) S]*EDF.Calib([1 chan+1],chan);
        elseif tmp==2
                S=S*ReRefMx(chan,:);
        elseif tmp==3
                S=[ones(size(S,1),1) S]*(EDF.Calib([1 chan+1],chan)*ReRefMx(chan,:));
        end;
end;

end; % GDF