0001 function [y,B,A]=nt_deglitch(x,glitchList,nfilt,nfit,advance);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 nt_greetings;
0019
0020 if nargin<2; error('!'); end
0021 if nargin<3; nfilt=[10,10]; end
0022 if nargin<4||isempty(nfit); nfit=50; end
0023 if nargin<5||isempty(advance); advance=3; end
0024
0025 if numel(nfilt)==1; nfilt=[nfilt,nfilt]; end
0026
0027 if ~6==exist('stmcb');
0028 error('Need function stmcb from signal toolbox');
0029 end
0030
0031
0032 y=x;
0033
0034 if size(x,2)==1 ~iscell(glitchList)
0035 glitchList={glitchList};
0036 end
0037
0038 for iChan=1:size(x,2)
0039
0040 events=glitchList{iChan};
0041
0042 for iEvent=1:numel(events)
0043
0044 iStart=events(iEvent)-advance;
0045 glitch=x(iStart+(1:nfit),iChan);
0046 glitch=glitch-mean(glitch(floor(nfit/2):end));
0047
0048
0049 try
0050 [B,A]=stmcb([glitch;zeros(size(glitch))],nfilt(1),nfilt(2));
0051 z=filter(B,A,[1;zeros(nfit-1,1)]);
0052 y(iStart+(1:nfit),iChan)=y(iStart+(1:nfit))-z;
0053 catch
0054 warning('stmcb failed (advance too large?)');
0055 end
0056 end
0057 end
0058
0059
0060 if ~nargout
0061
0062
0063 figure(1); clf; nt_banner('nt_destep');
0064 xx=[x(:),y(:)];
0065 lim=[min(xx(:)),max(xx(:))]; lim=[lim(1)-0.1*(lim(2)-lim(1)), lim(2)+0.1*(lim(2)-lim(1))];
0066 subplot 211; plot([x,y]); xlabel('samples'); ylim(lim); legend('raw','clean'); legend boxoff
0067 subplot 212; plot(y,'r'); xlabel('samples'); legend('clean'); legend boxoff
0068 figure(2); clf
0069 plot(filter(B,A,[1;zeros(199,1)]),'.-'); xlabel('samples'); title('estimated impulse response of antialiasing filter (last event)');
0070 clear y B A;
0071 end
0072
0073
0074
0075
0076 if 0
0077 N=8;
0078 Wn=0.25;
0079 nsamples=100;
0080 [B,A] = butter(N,Wn);
0081 y=filter(B,A,ones(nsamples,1));
0082 BETA0=[1, zeros(1,8), 1, zeros(1,8)];
0083 fun = @(beta,x)(filter([beta(1),beta(2),beta(3),beta(4),beta(5),beta(6),beta(7),beta(8),beta(9)],...
0084 [beta(10),beta(11),beta(12),beta(13),beta(14),beta(15),beta(16),beta(17),beta(18)],x));
0085 x=ones(nsamples,1);
0086 BETA = nlinfit(x,y,fun,BETA0);
0087 BETA=BETA/BETA(10);
0088 BB=BETA(1:9); AA=BETA(10:end);
0089 figure(1); clf;
0090 plot([y,filter(BB,AA,ones(size(y)))])
0091 end
0092
0093 if 0
0094 N=8;
0095 Wn=0.25;
0096 [B,A] = butter(N,Wn);
0097 x=[zeros(1100,1);ones(2000,1)];
0098 x=filter(B,A,x);
0099 x=x+0.01*randn(size(x));
0100 [y,stepList]=nt_destep(x);
0101 figure(1); clf;
0102 subplot 211;
0103 plot([x,y]); legend('raw','step removed')
0104 subplot 212;
0105 plot([y,nt_deglitch(y,stepList)]); legend('step removed', 'deglitched');
0106 end
0107 if 0
0108 N=8;
0109 Wn=0.25;
0110 [B,A] = butter(N,Wn);
0111 x=[zeros(1100,1);ones(1100,1);2*ones(1100,1);3*ones(1100,1);4*ones(1100,1);0*ones(1100,1)];
0112
0113 x=x+0.0*randn(size(x));
0114 nt_destep(x);
0115 end
0116 if 0
0117
0118 x=ft_read_data('/data/meg/litvak/phantom090715_BrainampDBS_20150709_18.ds'); x=permute(x,[2,1,3]); x=x(:,33:306,:);
0119 x=x(:,100,1);
0120 [y,stepList]=nt_destep(x,[],[],30);
0121 figure(1); subplot 211; plot(x);
0122 subplot 212; plot(nt_deglitch(y,stepList));
0123 end
0124
0125