sdfopen.m

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

                                        end;
                                else
                                        EDF.AFIR.channel1 = channel1;
                                end;
                                
                                if all(EDF.SIE.InChanSelect~=channel1)
                                        EDF.SIE.InChanSelect=[EDF.SIE.InChanSelect channel1];        
                                end;
                        end;
                end;


                if isempty(findstr(upper(arg4),'NOTCH50')) 
                        EDF.Filter.A=1;
                        EDF.Filter.B=1;
                else
                        EDF.SIE.FILT=1;
                        EDF.Filter.A=1;
                        EDF.Filter.B=1;
                        %if all(EDF.SampleRate(EDF.SIE.ChanSelect)==100)
                        if EDF.AS.MAXSPR/EDF.Dur==100
                                EDF.Filter.B=[1 1]/2;
                                %elseif all(EDF.SampleRate(EDF.SIE.ChanSelect)==200)
                        elseif EDF.AS.MAXSPR/EDF.Dur==200
                                EDF.Filter.B=[1 1 1 1]/4;
                                %elseif all(EDF.SampleRate(EDF.SIE.ChanSelect)==256)
                        elseif EDF.AS.MAXSPR/EDF.Dur==400
                                EDF.Filter.B=ones(1,8)/8;
                                %elseif all(EDF.SampleRate(EDF.SIE.ChanSelect)==256)
                        elseif EDF.AS.MAXSPR/EDF.Dur==256
                                EDF.Filter.B=poly([exp([-1 1 -2 2]*2*pi*i*50/256)]); %max(EDF.SampleRate(EDF.SIE.ChanSelect)))]);
                                EDF.Filter.B=EDF.Filter.B/sum(EDF.Filter.B);
                        else
                                fprintf(EDF.FILE.stderr,'Warning SDFOPEN: 50Hz Notch does not fit\n');
                        end;
                end;



                if ~isempty(findstr(upper(arg4),'NOTCH60')) 
                        fprintf(EDF.FILE.stderr,'Warning SDFOPEN: option NOTCH60 not implemented yet.\n');    
                end;


                if ~isempty(findstr(upper(arg4),'HPF')),  % high pass filtering
                        if EDF.SIE.FILT==0; EDF.Filter.B=1; end;
                        EDF.SIE.FILT=1;
                        EDF.Filter.A=1;
		end;
                if ~isempty(findstr(upper(arg4),'HPF0.0Hz')),  % high pass filtering
                        EDF.Filter.B=conv([1 -1], EDF.Filter.B);
                elseif ~isempty(findstr(upper(arg4),'TAU')),  % high pass filtering / compensate time constant
                        tmp=findstr(upper(arg4),'TAU');
                        TAU=strtok(upper(arg4(tmp:length(arg4))),'S');
                        tau=str2num(TAU);
                        if isempty(tau)
                                fprintf(EDF.FILE.stderr,'Warning SDFOPEN: invalid tau-value.\n');
                        else
                                EDF.Filter.B=conv([1 (EDF.Dur/EDF.AS.MAXSPR/tau-1)], EDF.Filter.B);
                        end;
			
                %%%% example 'HPF_1.0Hz_Hamming',  % high pass filtering
                elseif ~isempty(findstr(upper(arg4),'HPF')),  % high pass filtering
			    tmp=findstr(upper(arg4),'HPF');
			    FilterArg0=arg4(tmp+4:length(arg4));
			    %[tmp,FilterArg0]=strtok(arg4,'_');
			    [FilterArg1,FilterArg2]=strtok(FilterArg0,'_');
			    [FilterArg2,FilterArg3]=strtok(FilterArg2,'_');
			    tmp=findstr(FilterArg1,'Hz');
			    F0=str2num(FilterArg1(1:tmp-1));				    
			    B=feval(FilterArg2,F0*EDF.AS.MAXSPR/EDF.Dur);
			    B=B/sum(B);
			    B(ceil(length(B)/2))=(B(ceil(length(B)/2)))-1;
			    
                        EDF.Filter.B=conv(-B, EDF.Filter.B);
                end;


                if ~isempty(findstr(upper(arg4),'UNITS_BLOCK'))
			EDF.SIE.TimeUnits_Seconds=0; 
                end;

        end; % end nargin >3
        
        if EDF.SIE.FILT==1;
		EDF.Filter.Z=[];
                for k=1:length(EDF.SIE.ChanSelect),
                        [tmp,EDF.Filter.Z(:,k)]=filter(EDF.Filter.B,EDF.Filter.A,zeros(length(EDF.Filter.B+1),1));
                end;
    		EDF.FilterOVG.Z=EDF.Filter.Z;
        end;
        
        if EDF.SIE.REGC
                FN=[lower(EDF.FILE.Name) 'cov.mat'];
                if exist(FN)~=2
                        fprintf(EDF.FILE.stderr,'Warning %s: Covariance-file %s not found.\n',EDF.AS.Method,FN);
                        EDF.SIE.REGC=0;   
                else
                        if exist('OCTAVE_VERSION')==5
                                load(file_in_loadpath(FN));
                        else
                                load(FN);
                        end;
                        if exist('XC') == 1
                                %EDF.SIE.COV = tmp.XC;
                                %[N,MU,COV,Corr]=decovm(XC);
                                N=size(XC,2);
                                COV=(XC(2:N,2:N)/XC(1,1)-XC(2:N,1)*XC(1,2:N)/XC(1,1)^2);
                                
                                %clear tmp;
                                %cov = diag(EDF.Cal)*COV*diag(EDF.Cal);
                                mcov = M'*diag(EDF.Cal)*COV*diag(EDF.Cal);
                                %mcov(~())=0;
                                EDF.SIE.REG = eye(EDF.NS) - M*((mcov*M)\(mcov));
                                EDF.SIE.REG(channel1,channel1) = 1; % do not remove the regressed channels
                                %mcov, EDF.SIE.REG, 
                        else
                                fprintf(EDF.FILE.stderr,'Error %s: Regression Coefficients for ECG minimization not available.\n',EDF.AS.Method);
                        end;
                end;
                
        end;
        
        if EDF.SIE.TECG == 1; 
                % define channels that should be corrected
                if isfield(QRS,'Version')
                        OutChanSelect=[1:11 13:EDF.NS];
                        if EDF.SIE.REGC % correct templates
                                QRS.Templates=QRS.Templates*EDF.SIE.REG;
                                fprintf(EDF.FILE.stderr,'Warning SDFOPEN: Mode TECG+RECG not tested\n');
                        end;
                        if QRS.Version~=2
                                fprintf(EDF.FILE.stderr,'Warning SDFOPEN Mode TECG: undefined QRS-version\n');
                        end;
                else
                        %OutChanSelect=find(EDF.ChanTyp=='E' | EDF.ChanTyp=='O');
                        OutChanSelect=[1:9 ];
                        if any(EDF.SIE.ChanSelect>10)
                                fprintf(EDF.FILE.stderr,'Warning SDFOPEN: Mode TECG: Only #1-#9 are corrected\n');
                        end;
                        if EDF.SIE.REGC, % correct the templates
                                QRS.Templates=QRS.Templates*EDF.SIE.REG([1:9 12],[1:9 12]);
                                fprintf(EDF.FILE.stderr,'Warning SDFOPEN: Mode TECG+RECG not tested\n');
                        end;
                        
                end;
                fs = EDF.SPR(12)/EDF.Dur; 
                QRS.Templates=detrend(QRS.Templates,0); %remove mean
                EDF.TECG.idx = [(QRS.Index-fs/2-1) (EDF.NRec+1)*EDF.SPR(12)]; %include terminating element
                EDF.TECG.idxidx = 1; %pointer to next index

                % initialize if any spike is detected before first window    
	        pulse = zeros(length(QRS.Templates),1);
    		Index=[];
	        while EDF.TECG.idx(EDF.TECG.idxidx) < 1,
    		        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+length(QRS.Templates)) = 1;  
	        end;
        
                for i=1:length(EDF.SIE.InChanSelect),
                        k=find(OutChanSelect==EDF.SIE.InChanSelect(i));
                        if isempty(k)
                                EDF.TECG.QRStemp(:,i) = zeros(fs,1);
                        else
                                EDF.TECG.QRStemp(:,i) = QRS.Templates(0.5*fs:1.5*fs-1,k).*hanning(fs);
                        end;
                        [tmp,EDF.TECG.Z(:,i)] = filter(EDF.TECG.QRStemp(:,i),1,pulse);
                end;
        end; % if EDF.SIE.TECG==1
        
        %syms Fp1 Fp2 M1 M2 O2 O1 A1 A2 C3 C4
        for k=EDF.SIE.ChanSelect,
                %fprintf(1,'#%i: ',k);
                tmp=find(EDF.SIE.ReRefMx(:,k))';
                
                if EDF.SIE.ReRefMx(tmp(1),k)==1,
                        x=EDF.Label(tmp(1),:);
                else
                        x=sprintf('%3.1f*%s',EDF.SIE.ReRefMx(tmp(1),k),deblank(EDF.Label(tmp(1),:)));
                end;
                for l=2:length(tmp), L=tmp(l);
                        if (EDF.SPR(tmp(l-1),:)~=EDF.SPR(tmp(l),:))  
                                fprintf(EDF.FILE.stderr,'Warning %s: SampleRate Mismatch in "%s", channel #%i and #%i\n',upper(EDF.AS.Method),EDF.FILE.Name,tmp(l-1),tmp(l));
                        end;
                        if ~strcmp(EDF.PhysDim(tmp(l-1),:),EDF.PhysDim(tmp(l),:))  
                                fprintf(EDF.FILE.stderr,'Warning %s: Dimension Mismatch in "%s", channel #%i and #%i\n',upper(EDF.AS.Method),EDF.FILE.Name,tmp(l-1),tmp(l));
                        end;
                        if ~strcmp(EDF.Transducer(tmp(l-1),:),EDF.Transducer(tmp(l),:))  
                                fprintf(EDF.FILE.stderr,'Warning %s: Transducer Mismatch in "%s", channel #%i and #%i\n',upper(EDF.AS.Method),EDF.FILE.Name,tmp(l-1),tmp(l));
                        end;
                        if ~strcmp(EDF.PreFilt(tmp(l-1),:),EDF.PreFilt(tmp(l),:))  
                                fprintf(EDF.FILE.stderr,'Warning %s: PreFiltering Mismatch in "%s", channel #%i and #%i\n',upper(EDF.AS.Method),EDF.FILE.Name,tmp(l-1),tmp(l));
                        end;
                        x=[x sprintf('+(%3.1f)*(%s)',EDF.SIE.ReRefMx(tmp(l),k),deblank(EDF.Label(tmp(l),:)))];
                end;
                EDF.Label(k,1:length(x))=x; %char(sym(x))
        end;
        %Label,
        EDF.SIE.ReRefMx = EDF.SIE.ReRefMx(:,EDF.SIE.ChanSelect);
	if exist('OCTAVE_VERSION')>=5,
            EDF.Calib = (EDF.Calib*EDF.SIE.REG*EDF.SIE.ReRefMx); % important to be sparse, otherwise overflow-check does not work correctly.
	    fprintf(EDF.FILE.stderr,'Warning SDFOPEN: Overflow Check does not work without SPARSE\n');
	else
            EDF.Calib = sparse(EDF.Calib*EDF.SIE.REG*EDF.SIE.ReRefMx); % important to be sparse, otherwise overflow-check does not work correctly.
	end;
	        
        if 1, % ??? not sure, whether it has any advantage
        EDF.PhysDim = EDF.PhysDim(EDF.SIE.ChanSelect,:);
        EDF.PreFilt = EDF.PreFilt(EDF.SIE.ChanSelect,:);
        EDF.Transducer = EDF.Transducer(EDF.SIE.ChanSelect,:);
        end;
        
        if EDF.SIE.RS,
		tmp = EDF.AS.MAXSPR/EDF.Dur;
                if arg5==0
                        %arg5 = max(EDF.SPR(EDF.SIE.ChanSelect))/EDF.Dur;
                        
                elseif ~rem(tmp,arg5); % The target sampling rate must divide the source sampling rate   
			EDF.SIE.RS = 1;
			tmp=tmp/arg5;
			EDF.SIE.T = ones(tmp,1)/tmp;
                elseif arg5==100; % Currently, only the target sampling rate of 100Hz are supported. 
                        EDF.SIE.RS=1;
                        tmp=EDF.AS.MAXSPR/EDF.Dur;
                        if exist('OCTAVE_VERSION')
                                load resample_matrix4octave.mat T256100 T200100;
                        else
                                load('resample_matrix');
                        end;
                        if 1,
                                if tmp==400,
                                        EDF.SIE.T=ones(4,1)/4;
                                elseif tmp==256,
                                        EDF.SIE.T=T256100;
                                elseif tmp==200,
                                        EDF.SIE.T=T200100; 
                                elseif tmp==100,
                                        EDF.SIE.T=1;
                                else
                                        fprintf('Warning %s-READ: sampling rates should be equal\n',upper(EDF.AS.Method));     
                                end;	
                        else
                                tmp=EDF.SPR(EDF.SIE.ChanSelect)/EDF.Dur;
                                if all((tmp==256) | (tmp<100)) 
                                        EDF.SIE.RS = 1;
                                        %tmp=load(RSMN,'T256100');	
                                        EDF.SIE.T = T256100;	
                                elseif all((tmp==400) | (tmp<100)) 
                                        EDF.SIE.RS = 1;
                                        EDF.SIE.T = ones(4,1)/4;
                                elseif all((tmp==200) | (tmp<100)) 
                                        EDF.SIE.RS = 1;
                                        %tmp=load(RSMN,'T200100');	
                                        EDF.SIE.T = T200100;	
                                elseif all(tmp==100) 
                                        %EDF.SIE.T=load('resample_matrix','T100100');	
                                        EDF.SIE.RS=0;
                                else
                                        EDF.SIE.RS=0;
                                        fprintf('Warning %s-READ: sampling rates should be equal\n',upper(EDF.AS.Method));     
                                end;
                        end;
                else
                        fprintf(EDF.FILE.stderr,'Error %s-READ: invalid target sampling rate of %i Hz\n',upper(EDF.AS.Method),arg5);
                        EDF.SIE.RS=0;
			EDF.ErrNo=[EDF.ErrNo,];
			%EDF=sdfclose(EDF);
			%return;
                end;
        end;
        
        FN=[lower(EDF.FILE.Name), 'th.mat'];
        if exist(FN)~=2,
	        if EDF.SIE.TH, % && ~exist('OCTAVE_VERSION'),
                        fprintf(EDF.FILE.stderr,'Warning %s: THRESHOLD-file %s not found.\n',EDF.AS.Method,FN);
                        EDF.SIE.TH=0;   
                end;
        else
                if exist('OCTAVE_VERSION')==5
                        tmp=load(file_in_loadpath(FN));
                else
                        tmp=load(FN);
                end;
    	        if isfield(tmp,'TRESHOLD') 
        	        EDF.SIE.THRESHOLD = tmp.TRESHOLD;
                %else
            	        %fprintf(EDF.FILE.stderr,'Error %s: TRESHOLD''s not found.\n',EDF.AS.Method);
                end;
        end;
    	if EDF.SIE.TH>1, % Failing electrode detector 
	        fprintf(2,'Warning SDFOPEN: FED not implemented yet\n');
                for k=1:length(InChanSelect),K=InChanSelect(k);
	        %for k=1:EDF.NS,
   % 		        [y1,EDF.Block.z1{k}] = filter([1 -1], 1, zeros(EDF.SPR(K)/EDF.Dur,1));
    %        		[y2,EDF.Block.z2{k}] = filter(ones(1,EDF.SPR(K)/EDF.Dur)/(EDF.SPR(K)/EDF.Dur),1,zeros(EDF.SPR(K)/EDF.Dur,1));
                
    %       		[y3,EDF.Block.z3{k}] = filter(ones(1,EDF.SPR(K)/EDF.Dur)/(EDF.SPR(K)/EDF.Dur),1,zeros(EDF.SPR(K)/EDF.Dur,1));
    		end;
        end;

% Initialization of Bufferblock for random access (without EDF-blocklimits) of data 
	if ~EDF.SIE.RAW & EDF.SIE.TimeUnits_Seconds
                EDF.Block.number=[0 0 0 0]; %Actual Blocknumber, start and end time of loaded block, diff(EDF.Block.number(1:2))==0 denotes no block is loaded;
                                            % EDF.Blcok.number(3:4) indicate start and end of the returned data, [units]=samples.
		EDF.Block.data=[];
		EDF.Block.dataOFCHK=[];
	end;
        
%end; % end of SDFOPEN-READ



%%%%%%% ============= WRITE ===========%%%%%%%%%%%%        

elseif (arg2=='w') | (arg2=='w+')
%        fprintf(EDF.FILE.stderr,'error EDFOPEN: write mode not possible.\n'); 
        H1=[]; H2=[];
%        return;