📄 sdfopen.m
字号:
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=str2double(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=str2double(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 k=1:length(EDF.InChanSelect),
k=find(OutChanSelect==EDF.InChanSelect(k));
if isempty(k)
EDF.TECG.QRStemp(:,k) = zeros(fs,1);
else
EDF.TECG.QRStemp(:,k) = QRS.Templates(0.5*fs:1.5*fs-1,k).*hanning(fs);
end;
[tmp,EDF.TECG.Z(:,k)] = filter(EDF.TECG.QRStemp(:,k),1,pulse);
end;
end; % if EDF.SIE.TECG==1
%syms Fp1 Fp2 M1 M2 O2 O1 A1 A2 C3 C4
if 0, % ??? not sure, whether it has any advantage
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);
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 SDFOPEN: THRESHOLD-file %s not found.\n',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;
%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 (%s): FED not implemented yet\n',EDF.FileName);
for k=1:length(EDF.InChanSelect),K=EDF.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 any(arg2=='w') % | (arg2=='w+')
% fprintf(EDF.FILE.stderr,'error EDFOPEN: write mode not possible.\n');
H1=[]; H2=[];
% return;
EDF.SIE.RAW = 0;
if ~isstruct(arg1) % if arg1 is the filename
EDF.FileName=arg1;
if nargin<3
tmp=input('SDFOPEN: list of samplerates for each channel? ');
EDF.EDF.SampleRate = tmp(:);
else
EDF.EDF.SampleRate = arg3;
end;
EDF.NS=length(EDF.EDF.SampleRate);
if nargin<4
tmp=input('SDFOPEN: Duration of one block in seconds: ');
EDF.Dur = tmp;
EDF.SPR=EDF.Dur*EDF.SampleRate;
else
if ~isempty(findstr(upper(arg4),'RAW'))
EDF.SIE.RAW = 1;
else
EDF.Dur = arg4;
EDF.SPR=EDF.Dur*EDF.SampleRate;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -