0001 function [y,filt]=nt_destep(x,thresh,filt,guard,depth);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027 nt_greetings;
0028
0029 if nargin<1; error('!'); end
0030 if nargin<2||isempty(thresh); thresh=0.1; end
0031 if nargin<3; filt=[]; end
0032 if nargin<4||isempty(guard); guard=1000; end
0033 if nargin<5 || isempty(depth); depth=6; end
0034
0035 if isempty(filt)
0036 if ~6==exist('stmcb');
0037 error('Need function stmcb from signal toolbox (or nlinfit from statistics toolbox)');
0038 end
0039 end
0040 if ischar(filt)
0041 load(filt);
0042 if ~1==exist('A') || ~1==exist('B');
0043 error(['file >',filt,'< does not contain B and A']);
0044 end
0045 filt={B,A};
0046 end
0047
0048 y=x;
0049
0050 nfit=50;
0051 B=1; A=1;
0052
0053
0054 for iChan=1:size(x,2)
0055
0056 [splitI,score_vector,score]=nt_split(x(:,iChan),depth,thresh,guard);
0057
0058 if ~isempty(splitI)
0059 splitI=splitI-4;
0060 splitI=[0,splitI,size(x,1)];
0061 for iSplit=2:numel(splitI)-1
0062 m1=mean(y(splitI(iSplit-1)+1:splitI(iSplit),iChan));
0063 m2=mean(y(splitI(iSplit)+1:splitI(iSplit+1),iChan));
0064 step=(m2-m1);
0065 step_response=y(splitI(iSplit):splitI(iSplit+1),iChan);
0066 if size(step_response,1)>nfit; step_response=step_response(1:nfit); end
0067 step_response=[step_response;ones(size(step_response,1),1)*mean(step_response(round(size(step_response,1)/2):end))];
0068
0069 if isempty(filt)
0070
0071 The filter is estimated anew for each step and each
0072 channel.
0073
0074 xx=ones(size(step_response,1),1)*step; yy=[step_response];
0075 if 1
0076
0077 [B,A]=stmcb(diff([step_response]/step),10,10);
0078 else
0079 BETA0=[1,zeros(1,5),1,zeros(1,5)];
0080 fun = @(beta,x)(filter([beta(1),beta(2),beta(3),beta(4),beta(5),beta(6)],...
0081 [beta(7),beta(8),beta(9),beta(10),beta(11),beta(12)],x));
0082 BETA = nlinfit(xx,yy,fun,BETA0);
0083 BETA=BETA/BETA(7);
0084 B=BETA(1:6); A=BETA(7:end);
0085 end
0086 else
0087 B=filt{1}; A=filt{2};
0088 end
0089
0090 y(splitI(iSplit)+1:end,iChan)=...
0091 y(splitI(iSplit)+1:end,iChan)-...
0092 filter(B,A,step*ones(size(y,1)-splitI(iSplit),1));
0093 end
0094 end
0095 end
0096
0097
0098 if ~nargout
0099
0100 disp(['steps at: ', num2str(splitI(2:end-1))]);
0101 figure(1); clf; nt_banner('nt_destep');
0102 xx=[x(:),y(:)];
0103 lim=[min(xx(:)),max(xx(:))]; lim=[lim(1)-0.1*(lim(2)-lim(1)), lim(2)+0.1*(lim(2)-lim(1))];
0104 subplot 211; plot([x,y]); xlabel('samples'); ylim(lim); legend('raw','clean'); legend boxoff
0105 subplot 212; plot(y,'r'); xlabel('samples'); legend('clean'); legend boxoff
0106 figure(2); clf
0107 plot(filter(B,A,ones(200,1)),'.-'); xlabel('samples'); title('estimated step response of antialiasing filter (last step)');
0108 clear y, filt;
0109 end
0110
0111
0112
0113
0114 if 0
0115 N=8;
0116 Wn=0.25;
0117 nsamples=100;
0118 [B,A] = butter(N,Wn);
0119 y=filter(B,A,ones(nsamples,1));
0120 BETA0=[1, zeros(1,8), 1, zeros(1,8)];
0121 fun = @(beta,x)(filter([beta(1),beta(2),beta(3),beta(4),beta(5),beta(6),beta(7),beta(8),beta(9)],...
0122 [beta(10),beta(11),beta(12),beta(13),beta(14),beta(15),beta(16),beta(17),beta(18)],x));
0123 x=ones(nsamples,1);
0124 BETA = nlinfit(x,y,fun,BETA0);
0125 BETA=BETA/BETA(10);
0126 BB=BETA(1:9); AA=BETA(10:end);
0127 figure(1); clf;
0128 plot([y,filter(BB,AA,ones(size(y)))])
0129 end
0130
0131 if 0
0132 N=8;
0133 Wn=0.25;
0134 [B,A] = butter(N,Wn);
0135 x=[zeros(1100,1);ones(2000,1)];
0136 x=filter(B,A,x);
0137 x=x+0.01*randn(size(x));
0138 nt_destep(x);
0139 end
0140 if 0
0141 N=8;
0142 Wn=0.25;
0143 [B,A] = butter(N,Wn);
0144 x=[zeros(1100,1);ones(1100,1);2*ones(1100,1);3*ones(1100,1);4*ones(1100,1);0*ones(1100,1)];
0145 x=filter(B,A,x);
0146 x=x+0.0*randn(size(x));
0147 nt_destep(x);
0148 end
0149 if 0
0150 x=ft_read_data('/data/meg/phantom090715_BrainampDBS_20150709_18.ds'); x=permute(x,[2,1,3]); x=x(:,33:306,:);
0151 x=x(:,100,1);
0152 nt_destep(x,[],[],30);
0153 end
0154
0155