sdfread.m

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

                %EDF.AFIR.xin=S(:,find(EDF.AFIR.channel1==EDF.SIE.InChanSelect)); 
                EDF.AFIR.xin=S(:,EDF.AFIR.channel1); 
        end;
        
        % Overflow Check for SDF
        if nargout>2, %
                OFCHK=S(:,InChanSelect);
                if EDF.SIE.TH
                        for k=1:length(InChanSelect),K=InChanSelect(k);
                                OFCHK(:,k)=(S(:,k) < EDF.SIE.THRESHOLD(K,1)) + (S(:,k)>= EDF.SIE.THRESHOLD(K,2));
                        end;
                else 
                        fprintf(2,'Error SDFREAD (#SDF): No Overflow Threshold defined\n');
                end;
                OFCHK=OFCHK*abs(EDF.Calib([InChanSelect+1],:)~=0); %RS%
        end;
        if EDF.SIE.TH
                for k=1:length(InChanSelect),K=InChanSelect(k);
                        tmp=find(((S(:,k) < EDF.SIE.THRESHOLD(K,1)) + (S(:,k)>= EDF.SIE.THRESHOLD(K,2)))>0);
        %                S(tmp,k)=NaN;
                        %S(find(((S(:,k) < EDF.SIE.THRESHOLD(K,1)) + (S(:,k)>= EDF.SIE.THRESHOLD(K,2)))>0),k)=NaN;
                end;
        end;
        
else %%%%%%%%% if all channels of same type
        
datatyp=dt(EDF.GDFTYP(1));        
if Mode_RAW;
        [S, count]=fread(EDF.FILE.FID,[EDF.AS.spb, Records],datatyp);
        
        %count = floor(count/EDF.AS.spb);
        count = (count/EDF.AS.spb);
        if count<1; fprintf(2,'Warning EDFREAD: only %3.1f blocks were read instead  of %3.1f\n',count,Records); end;

        % Overflow Check for EDF-RAW
        %if nargout>2, %
                OFCHK = S;
                if EDF.SIE.TH>1
                        for k=1:length(InChanSelect),K=InChanSelect(k);
		                [y1,EDF.Block.z1{k}]=filter([1 -1],1,tmp,EDF.Block.z1{k});
            			[y2,EDF.Block.z2{k}]=filter(ones(1,EDF.SPR(K)/EDF.Dur)/(EDF.SPR(K)/EDF.Dur),1,y1==0,EDF.Block.z2{k});
            			%y2=y2>0.079*(EDF.SPR(k)/EDF.Dur);
    				[y3,EDF.Block.z3{k}]=filter(ones(1,EDF.SPR(K)/EDF.Dur)/(EDF.SPR(K)/EDF.Dur),1,(tmp>=EDF.SIE.THRESHOLD(K,2)) | (tmp<=EDF.SIE.THRESHOLD(K,1)),EDF.Block.z3{k});
            			%y3=y3>0.094*(EDF.SPR(k)/EDF.Dur);
				OFCHK(bi(k)+1:bi(k+1),:) = y3>0.094 | y2>0.079;
			end;
                elseif EDF.SIE.TH
                        for k=1:length(InChanSelect),K=InChanSelect(k);
                                OFCHK(bi(k)+1:bi(k+1),:)=(S(bi(k)+1:bi(k+1),:) < EDF.SIE.THRESHOLD(K,1)) + (S(bi(k)+1:bi(k+1),:)>= EDF.SIE.THRESHOLD(K,2));
                        end;
                        %S(OFCHK>0)=NaN;
                elseif nargout>2, %
                        fprintf(2,'Error SDFREAD (#EDF-RAW): No Overflow Threshold defined\n');
                end;
        %end;
        
	if Mode_REF
	    fprintf(2,'Warning SDFREAD: ReReferenzing "R" is not possible in combination with RAW "W"\n');
	end;
else
        if all(EDF.SPR(InChanSelect)==EDF.AS.MAXSPR)
                if ~OptiMEM % but OptiSPEED
                        [s, count]=fread(EDF.FILE.FID,[EDF.AS.spb,Records],datatyp);
                        count = (count/EDF.AS.spb);
                        S = zeros(maxspr*count,length(InChanSelect)); %%%
                        for k=1:length(InChanSelect), K=InChanSelect(k);
                                S(:,k)=reshape(s(bi(K)+1:bi(K+1),:),maxspr*count,1); %RS%
                        end;

                else %OptiMEM % but ~OptiSPEED
                        S = zeros(maxspr*Records,length(InChanSelect));%%%
                        idx0=0; count=0;
                        for l=1:ceil(Records/MAX_BLOCK_NUMBER),
                                tmp_norr=min(MAX_BLOCK_NUMBER,Records-count); % 16 = # of blocks read 
                                [s, C]=fread(EDF.FILE.FID,[EDF.AS.spb,tmp_norr],datatyp);
                                C = C/EDF.AS.spb;
                                for k=1:length(InChanSelect), K=InChanSelect(k);
                                        S(idx0+(1:maxspr*C),k) = reshape(s(bi(K)+1:bi(K+1),:),maxspr*C,1); %RS%
                                end;
                                idx0 = idx0 + maxspr*C;
                                count=count+C;
                        end;
                        S=S(1:idx0,:);
                end;
                %AFIR%
                if EDF.SIE.AFIR
                        EDF.AFIR.xin=S(:,EDF.AFIR.channel1); 
                end;
                
                % Overflow Check for EDF-all same sampling rate
                if nargout>2, %
                        OFCHK = zeros(size(S,1),length(InChanSelect));
                        if EDF.SIE.TH>1
                    		for k=1:length(InChanSelect),K=InChanSelect(k);
		            		[y1,EDF.Block.z1{k}]=filter([1 -1],1,tmp,EDF.Block.z1{k});
            				[y2,EDF.Block.z2{k}]=filter(ones(1,EDF.SPR(K)/EDF.Dur)/(EDF.SPR(K)/EDF.Dur),1,y1==0,EDF.Block.z2{k});
            				%y2=y2>0.079*(EDF.SPR(k)/EDF.Dur);
            				[y3,EDF.Block.z3{k}]=filter(ones(1,EDF.SPR(K)/EDF.Dur)/(EDF.SPR(K)/EDF.Dur),1,(tmp>=EDF.SIE.THRESHOLD(K,2)) | (tmp<=EDF.SIE.THRESHOLD(K,1)),EDF.Block.z3{k});
            				%y3=y3>0.094*(EDF.SPR(k)/EDF.Dur);
					OFCHK(bi(k)+1:bi(k+1),:) = y3>0.094 | y2>0.079;
				end;
            		elseif EDF.SIE.TH,
                                for k=1:length(InChanSelect),K=InChanSelect(k);
                                        OFCHK(:,k)=(S(:,k) <= EDF.SIE.THRESHOLD(K,1)) + (S(:,K)>= EDF.SIE.THRESHOLD(K,2));
                                end;
                                OFCHK = OFCHK*abs(EDF.Calib([InChanSelect+1],:)~=0); %RS%
                        else 
                                fprintf(2,'Error SDFREAD (#EDF-SPR=): No Overflow Threshold defined\n');
                        end;
                end;
                
                if 0; EDF.SIE.TH
                        for k=1:length(InChanSelect),K=InChanSelect(k);
                                tmp = S(:,k);
                                tmp = find(((tmp < EDF.SIE.THRESHOLD(K,1)) + (tmp >= EDF.SIE.THRESHOLD(K,2)))>0);
                %                S(tmp,k) = NaN;
                        end;
                end;
        else
                if ~OptiMEM % but OptiSPEED
                        [s, count]=fread(EDF.FILE.FID,[EDF.AS.spb,Records],datatyp);
                        count = (count/EDF.AS.spb);
                        S = zeros(maxspr*count,length(InChanSelect));
                        for k=1:length(InChanSelect), K=InChanSelect(k);
                                tmp=reshape(s(bi(K)+1:bi(K+1),:),EDF.SPR(K)*count,1);
                                if EDF.SPR(K) == maxspr
                                        S(:,k)=tmp;%reshape(tmp(:,ones(1,maxspr/EDF.SPR(K)))',maxspr*Records,1); %RS%
                                elseif EDF.SPR(K) > maxspr
                                        if rem(EDF.SPR(K)/maxspr,1)==0
                                                S(:,k)=rs(tmp,EDF.SPR(K),maxspr);
                                        else
                                                S(:,k)=rs(tmp,EDF.SIE.T); %(:,ones(1,maxspr/EDF.SPR(K)))',maxspr*Records,1); %RS%
                                        end;
                                elseif EDF.SPR(K) < maxspr
                                        S(:,k)=reshape(tmp(:,ones(1,maxspr/EDF.SPR(K)))',maxspr*count,1); %RS%
                                end;
                        end;    
                else %OptiMEM % but ~OptiSPEED
                        S = zeros(maxspr*Records,length(InChanSelect));
                        idx0=0; count=0;
                        for l=1:ceil(Records/MAX_BLOCK_NUMBER),
                        %while count < Records
                                tmp_norr = min(MAX_BLOCK_NUMBER,Records-count); % 16 = # of blocks read 
                                [s, C]=fread(EDF.FILE.FID,[EDF.AS.spb,tmp_norr],datatyp);
                                C = C/EDF.AS.spb;
                                for k=1:length(InChanSelect), K=InChanSelect(k); 
                                        tmp0=reshape(s(bi(K)+1:bi(K+1),:),EDF.SPR(K)*C,1);
                                        if EDF.SPR(K)==maxspr
                                                S(idx0+(1:maxspr*C),k)=tmp0;%reshape(tmp(:,ones(1,maxspr/EDF.SPR(K)))',maxspr*Records,1); %RS%
                                        elseif EDF.SPR(K)>maxspr
                                                if rem(EDF.SPR(K)/maxspr,1)==0
                                                        S(idx0+(1:maxspr*C),k)=rs(tmp0,EDF.SPR(K),maxspr);
                                                else
                                                        S(idx0+(1:maxspr*C),k)=rs(tmp0,EDF.SIE.T); %(:,ones(1,maxspr/EDF.SPR(K)))',maxspr*Records,1); %RS%
                                                end;
                                        elseif EDF.SPR(K)<maxspr
                                                S(idx0+(1:maxspr*C),k)=reshape(tmp0(:,ones(1,maxspr/EDF.SPR(K)))',maxspr*C,1); %RS%
                                        end;
                                end;
                                idx0 = idx0 + maxspr*C;
                                count= count+C;
                        end;
                        S=S(1:idx0,:);
                end;
                
                %AFIR%
                if EDF.SIE.AFIR
                        %EDF.AFIR.xin=S(:,find(EDF.AFIR.channel1==InChanSelect)); 
                        EDF.AFIR.xin=S(:,EDF.AFIR.channel1); 
                end;
                
                % Overflow Check for EDF with different sampling rates
                if nargout>2, %
                        if Mode_RS100
                                fprintf(2,'Warning SDFREAD: Overflow check for EDF files when reading channels with different sampling rate and resampling, is not implemented yet\n');
                        end;
                        OFCHK = S(:,InChanSelect);
                        if EDF.SIE.TH
                                for k=1:length(InChanSelect), K=InChanSelect(k);
                                        OFCHK(:,k)=(S(:,k) < EDF.SIE.THRESHOLD(K,1)) + (S(:,K)>= EDF.SIE.THRESHOLD(K,2));
                                end;
                        else 
                                fprintf(2,'Error SDFREAD (#EDF-SPR~): No Overflow Threshold defined\n');
                        end;
                end;
                if EDF.SIE.TH
                        for k=1:length(InChanSelect),K=InChanSelect(k);
                                tmp = S(:,k);
                                tmp = find(((tmp < EDF.SIE.THRESHOLD(K,1)) + (tmp >= EDF.SIE.THRESHOLD(K,2)))>0);
                %                S(tmp,k) = NaN;
                        end;
                end;
        end;
end;
end; % SDF

EDF.AS.numrec=count;
EDF.FILE.POS = EDF.AS.startrec + EDF.AS.numrec;
if EDF.AS.numrec~=Records, 
        fprintf(2,'Warning %s: %s only %i blocks instead of %i read\n',mfilename,EDF.FILE.Name,EDF.AS.numrec,Records);
end;

%%%%% Calibration of the signal 
if EDF.SIE.RR          % Autocalib 
        if Mode_RAW 
                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;
        else
               %S=[ones(size(S,1),1) S(:,InChanSelect)]*EDF.Calib([1 InChanSelect+1],:);
                S=[ones(size(S,1),1) S]*EDF.Calib([1 InChanSelect+1],:); % EDF.Calib must be sparse, otherwise overflow-check is incorrect. 
        end;
end;

%%%%% Removing ECG Templates
if EDF.SIE.TECG          
        if (EDF.AS.startrec+EDF.AS.numrec)~=(ftell(EDF.FILE.FID)-EDF.HeadLen)/EDF.AS.bpb;
                fprintf(2,'ERROR SDFREAD: Mode TECG requires update of EDF [S,EDF]=sdfread(EDF,...)\n');
                EDF=sdftell(EDF);
        end;
        pulse = zeros(EDF.AS.numrec*EDF.SPR(12),1);
        
        Index=[];
        while EDF.TECG.idx(EDF.TECG.idxidx) <= EDF.FILE.POS*EDF.SPR(12)
                Index=[Index EDF.TECG.idx(EDF.TECG.idxidx)-EDF.AS.startrec*EDF.SPR(12)];
                EDF.TECG.idxidx=EDF.TECG.idxidx+1;
        end;

        if ~isempty(Index)
                pulse(Index) = 1;
        end;
        
        %tmp=find(EDF.TECG.idx > EDF.AS.startrec*EDF.SPR(12) & EDF.TECG.idx <= EDF.FILE.POS*EDF.SPR(12));
        %if ~isempty(tmp)
        %        pulse(EDF.TECG.idx(tmp)-EDF.AS.startrec*EDF.AS.MAXSPR) = 1;
        %end;
        
        for i=1:size(S,2),
                [tmp,EDF.TECG.Z(:,i)] = filter(EDF.TECG.QRStemp(:,i),1,pulse,EDF.TECG.Z(:,i));
                S(:,i)=S(:,i)-tmp; % corrects the signal
        end;
end;

%%%%% Filtering
if EDF.SIE.FILT
        for k=1:size(S,2);
                [S(:,k),EDF.Filter.Z(:,k)]=filter(EDF.Filter.B,EDF.Filter.A,S(:,k),EDF.Filter.Z(:,k));
        end;
        % Resampling of Overflow Check 
        if nargout>2, %
                if EDF.SIE.TH,
                        for k=1:size(S,2);
                                [S(:,k),EDF.FilterOVG.Z(:,k)]=filter(EDF.Filter.B,EDF.Filter.A,S(:,k),EDF.FilterOVG.Z(:,k));
                        end;
                end;
        end;
end;

%AFIR%
% Implements Adaptive FIR filtering for ECG removal in EDF/SDF-tb.
% based on the Algorithm of Mikko Koivuluoma <[e-mail]>
if EDF.SIE.AFIR
        e=zeros(size(S,1),length(EDF.AFIR.channel2));
        
        xinre2 = [EDF.AFIR.x EDF.AFIR.xin'];
        
        ddin=[EDF.AFIR.d; S]; %(:,EDF.AFIR.channel2)];