sdfread.m

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

% sdfread() -  Reads selected seconds of an EDF File (European Data Format for Biosignals)
%              The EDF data and header format is specified in [1].
% Usage:
%   [S,EDF] = sdfread(EDF [,SecRead [,StartSec]] )
%
% Inputs:
%   EDF      = Opened EDF file header structure returned by sdfopen();  
%   SecRead  = Number of seconds to read. Inf reads the whole file. Default is 1 (sec). 
%   StartSec = Starting position in seconds. If not provided, the data is read
%              continously from the EDF file. No file pointer repositioning is performed.
% Example:
%          >> EDF     = sdfopen(Filename,'r',CHAN,ReRefMx,TSR,OFCHK); % see sdfopen() 
% % then
%          >> [S,EDF] = sdfread(EDF, SecRead, StartSec) % Read SecRead seconds of data 
%                                                       % beginning at StartSec, else
%          >> [S,EDF] = sdfread(EDF, Inf);              % read the whole file data, or
%          >> [S,EDF] = sdfread(EDF, EDF.NRec*EDF.Dur); % equivalent...
%
% Note:
% 1) Versions>0.75 request SecRead and StartSec in seconds. 
%    Previously, Versions<0.76, the units were Records. 
% 2) In this version, modified for eeglab(), functions sdfseek(), sdftell(), 
%    sdfeof() and dt() have been included in the code of sdfread(). -AD
%
% Author: (C) 1997-2002 by Alois Schloegl, 15 Jun 2002 #0.85, (Header reworked for 
%         EEGLAB format and subfunctions inserted in the code, Arnaud Delorme and 
%         Scott Makeig, 27 Dec 2002)
%
% See also: fread, SDFOPEN, SDFWRITE, SDFCLOSE, SDFREWIND, [SDFSEEK, SDFTELL, SDFEOF]

%	Copyright (c) 1997-2002 by Alois Schloegl
%	[e-mail]	
%	Version 0.85
%	24. Jun. 2002

% References:
% [1] Bob Kemp, Alpo Värri, Agostinho C. Rosa, Kim D. Nielsen and John Gade.
%     A simple format for exchange of digitized polygraphic recordings.
%     Electroencephalography and Clinical Neurophysiology, 82 (1992) 391-393.
% See also: http://www.medfac.leidenuniv.nl/neurology/knf/kemp/edf/edf_spec.htm

% .changelog
% 13022002  V0.84    minor bug fixed
% 24062002  V0.85    check of input arguments corrected

function [S,EDF,OFCHK] = sdfread(EDF,NoS,StartPos)
global SDF_TB_TOGGLE_CHECK

if isfield(EDF.AS,'TOGGLE');
        if EDF.AS.TOGGLE~=SDF_TB_TOGGLE_CHECK
                fprintf(2,'Warning SDFREAD: [s,EDF]=sdfread(EDF, ...) \nYou forgot to pass EDF in %i call(s) of SDFREAD\n',SDF_TB_TOGGLE_CHECK-EDF.AS.TOGGLE);
        end;
else
        EDF.AS.TOGGLE=0;
        SDF_TB_TOGGLE_CHECK=0;
end;
SDF_TB_TOGGLE_CHECK = SDF_TB_TOGGLE_CHECK+1;
EDF.AS.TOGGLE = EDF.AS.TOGGLE+1;

OptiMEM=1;             % if you get low on memory (so your CPU starts swapping), 
                       % set this to one
OptiSPEED=~OptiMEM;
MAX_BLOCK_NUMBER=16;	% max. # of blocks in case of OptiMEM

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

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

if isfield(EDF,'SIE')   % SIESTA; Layer 4
        Mode_RAW=EDF.SIE.RAW; %any(Mode==[[4:7 12:15] [4:7 12:15]+16]);
	%Mode_CAL=EDF.SIE.RR; %any(Mode==[[2 3 6 7 10 11 14 15] [2 3 6 7 10 11 14 15]+16]);
        Mode_REF=0; %any(Mode==[[8:15] [8:15]+16]);
	Mode_SEC=0; %any(Mode==[[1:2:15] [1:2:15]+16]);
        Mode_RS100=EDF.SIE.RS; %bitand(Mode,2^4)>0;
        
        EDF.AS.startrec=EDF.AS.startrec+EDF.AS.numrec;
        EDF.AS.numrec=0;
        
        % prepare input arguments
        if ~EDF.SIE.TimeUnits_Seconds;        % units of time in records      
                NoR=NoS;%chan;
                if ~(nargin<3) 
                        [EDF]=sdfseek(EDF,StartPos,'bof'); 
                end;
                if isinf(NoR), 
                        NoR=EDF.NRec-EDF.FILE.POS; 
                end;
        else         % units of time in seconds
                % Prepare input arguments, transform NoS and StartPos into units of samples 
                if (nargin>2) %nargin==3 or larger
                        %tmp = floor(StartPos/EDF.Dur);
                        if StartPos+NoS > EDF.Dur*EDF.NRec;
                                fprintf(2,'Warning EDFREAD: %s has only %i seconds\n',EDF.FileName, EDF.Dur*EDF.NRec);
                                StartPos=min(StartPos*EDF.AS.MAXSPR/EDF.Dur,EDF.AS.MAXSPR*EDF.NRec); % transformation of seconds to samples + Overflowcheck
                                NoS=EDF.AS.MAXSPR*EDF.NRec-StartPos;
                        else
                                StartPos = StartPos*EDF.AS.MAXSPR/EDF.Dur;% EDF.Block.number(4);%/EDF.AS.MAXSPR*EDF.Dur;
                                NoS=NoS*EDF.AS.MAXSPR/EDF.Dur;
                        end;
                        NoR = ceil((StartPos+NoS)/EDF.AS.MAXSPR)-floor(StartPos/EDF.AS.MAXSPR);
                        %EDF.AS.A_Block_number=0;
                elseif nargin==2 %nargin==2
                        StartPos = EDF.Block.number(4); % EDF.Block.number(4);%/EDF.AS.MAXSPR*EDF.Dur;
                        NoS = NoS * EDF.AS.MAXSPR/EDF.Dur; % 
                else  % nargin<2
                        %NoR = 1, %ceil(NoS/EDF.Dur);
                        NoS = EDF.AS.MAXSPR/EDF.Dur; %EDF.Dur; % default one second
                        StartPos = EDF.Block.number(4);%/EDF.AS.MAXSPR*EDF.Dur; % EDF.FILE.POS*EDF.Dur;
                        %StartPos = EDF.FILE.POS*EDF.Dur,
                end;
                if isinf(NoS), 
                        NoS=EDF.NRec*EDF.AS.MAXSPR-EDF.Block.number(4); 
                end;
                
                % Q? whether repositioning is required, otherwise read continously, 
                if floor(StartPos/EDF.AS.MAXSPR)~=EDF.FILE.POS; % if start not next block 
                        if floor(StartPos/EDF.AS.MAXSPR)~=floor(EDF.Block.number(1)/EDF.AS.MAXSPR); % if start not same as Block.data
                                [EDF]=sdfseek(EDF,floor(StartPos/EDF.AS.MAXSPR),'bof');
                        end;
                end;
                % Q whether last bufferblock is needed or not
                if (StartPos >= EDF.Block.number(2)) | (StartPos < EDF.Block.number(1)) | diff(EDF.Block.number(1:2))==0
                        EDF.Block.number=[0 0 0 0 ];
                        EDF.Block.data=[];
                        NoR=ceil((StartPos+NoS)/EDF.AS.MAXSPR)-floor(StartPos/EDF.AS.MAXSPR);
                else
                        NoR=ceil((StartPos+NoS-EDF.Block.number(2))/EDF.AS.MAXSPR);
                end;
                
        end;
        
        
        %clear chan;
	InChanSelect=EDF.SIE.InChanSelect;
        Mode_CHANSAME = ~all(EDF.SPR(InChanSelect)==EDF.SPR(InChanSelect(1)));
        clear Mode;
        
else % Layer 3
	%if chan==0        
                InChanSelect=1:EDF.NS;
        %else
        %        InChanSelect=chan;
	%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] [4:7 12:15]+16]);
        EDF.SIE.RR=any(Mode==[[2 3 6 7 10 11 14 15] [2 3 6 7 10 11 14 15]+16]);
        Mode_REF=any(Mode==[[8:15] [8:15]+16]);
        Mode_SEC=any(Mode==[[1:2:15] [1:2:15]+16]);
        Mode_RS100=0; %bitand(Mode,2^4)>0;
        EDF.SIE.FILT=0;
        
	% Position file pointer and calculate number of Records
	if Mode_SEC
    		if ~all(~rem([NoR StartPos],EDF.Dur))
            		fprintf(2,'Warning EDFREAD: NoR and/or 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) 
    	        EDF=sdfseek(EDF,StartPos,'bof'); 
        else
                EDF.AS.startrec=EDF.AS.startrec+EDF.AS.numrec;
                EDF.AS.numrec=0;
	end;

	if Mode_REF
    		ReRefMx=chan;
                InChanSelect=find(any(ReRefMx,1))';
	else 
                InChanSelect=reshape(chan,1,prod(size(chan)));
                ReRefMx=eye(EDF.NS);
                ReRefMx=ReRefMx(:,chan);
	end;

        tmp = EDF.SIE.RR + Mode_REF*2;
        if tmp==0
                EDF.Calib = [zeros(1,EDF.NS);eye(EDF.NS)];
                EDF.Calib = EDF.Calib(:,chan);
                %[zeros(1,EDF.NS); eye(EDF.NS)];
        elseif tmp==1
                %EDF.Calib=EDF.Calib;
                EDF.Calib = EDF.Calib(:,chan);
                %S=[ones(size(S,1),1) S]*EDF.Calib([1 chan+1],chan);
        elseif tmp==2
                EDF.Calib=[zeros(1,length(chan)); ReRefMx];
                %S=S*ReRefMx(chan,:);
        elseif tmp==3
                EDF.Calib=(EDF.Calib*ReRefMx);
                %S=[ones(size(S,1),1) S]*(EDF.Calib([1 chan+1],chan)*ReRefMx(chan,:));
        end;
        
	if exist('OCTAVE_VERSION')~=1
        	InChanSelect=find(any(EDF.Calib(2:EDF.NS+1,:),2))';
        else
                InChanSelect=find(any(EDF.Calib(2:EDF.NS+1,:)'));
        end;
        
        Mode_CHANSAME = ~all(EDF.SPR(InChanSelect)==EDF.SPR(InChanSelect(1)));
	%if any(EDF.SPR(chan)~=EDF.SPR(chan(1))) fprintf(2,'Warning EDFREAD: channels do not have the same sampling rate\n');end;
        clear chan;
end;
        
bi=[0;cumsum(EDF.SPR)];
Records=NoR;
maxspr=EDF.AS.MAXSPR; %max(EDF.SPR(EDF.SIE.ChanSelect));
count=0;

if Mode_RAW;
        S=zeros(sum(EDF.SPR(InChanSelect)),Records);
        bi=[0;cumsum(EDF.SPR(InChanSelect))];
else
        S=zeros(maxspr*Records,length(InChanSelect));
end;        

if ~EDF.AS.SAMECHANTYP; % all(EDF.GDFTYP(:)~=EDF.GDFTYP(1))
        idx0=0;
        count=0;
        while (count<Records*EDF.AS.spb) & ~sdfeof(EDF), 
                %   for K=1:length(InChanSelect), K=InChanSelect(k);
                for k=1:EDF.NS,
                        %if GDF-format should be used 
                        %datatyp=dt(EDF.GDFTYP(k));
                        %if exist('OCTAVE_VERSION')
                        %        tmp1=(datatyp~='');
                        %else
                        %        tmp1=~isempty(datatyp);
                        %end;
                        %if tmp1
                        %        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),datatyp);
                        %else 
                        %        fprintf(2,'Error SDFREAD: Invalid SDF channel type in %s at channel %i',EDF.FileName,k);
                        %end;
                        
                        [tmp,cnt]=fread(EDF.FILE.FID,EDF.SPR(k),int16);
                        if any(InChanSelect==k), K=find(InChanSelect==k);
                                if Mode_RAW;
                                        S(bi(K)+1:bi(K+1),count+1)=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(idx0+(1:maxspr),K)=tmp; %RS%
                                end;  
                        end;      
                        count=count+cnt;
                end; 
                idx0 = idx0 + maxspr;
        end;
        S=S(1:idx0,:);
        count=count/EDF.AS.spb;
        
        %AFIR%
        if EDF.SIE.AFIR