Home > NoiseTools > nt_destep.m

nt_destep

PURPOSE ^

[y,B,A]=nt_destep(x,thresh,filt,guard,depth) - remove step glitch from MEG data

SYNOPSIS ^

function [y,filt]=nt_destep(x,thresh,filt,guard,depth);

DESCRIPTION ^

[y,B,A]=nt_destep(x,thresh,filt,guard,depth) - remove step glitch from MEG data

  y: clean data
  filt: estimated filter coefficients (filt={B,A})

  x: data to clean (time * channels)
  thresh: threshold for variance reduction [default: 0.1]
  filt: antialiasing filter, if [] estimate [default:[]]
  guard: minimum duration of stable interval in samples [default: 1000]
  depth: recursion depth for nt_split [default:6], determines number of steps

 Searches for large steps, removes them if variance ratio less than 
 "thresh", and both intervals shorter than "guard".

 Parameter "filt" can be:
  []: estimate from the data
  {B,A}: filter coefficients
  file name: mat file containg B and A
 
 Estimation from data might sometimes give erroneous results. Consider running
 nt_destep on a clean step, save the filter estimates, and use them.
 
 See also nt_detrend_robust.

 NoiseTools.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [y,filt]=nt_destep(x,thresh,filt,guard,depth);
0002 %[y,B,A]=nt_destep(x,thresh,filt,guard,depth) - remove step glitch from MEG data
0003 %
0004 %  y: clean data
0005 %  filt: estimated filter coefficients (filt={B,A})
0006 %
0007 %  x: data to clean (time * channels)
0008 %  thresh: threshold for variance reduction [default: 0.1]
0009 %  filt: antialiasing filter, if [] estimate [default:[]]
0010 %  guard: minimum duration of stable interval in samples [default: 1000]
0011 %  depth: recursion depth for nt_split [default:6], determines number of steps
0012 %
0013 % Searches for large steps, removes them if variance ratio less than
0014 % "thresh", and both intervals shorter than "guard".
0015 %
0016 % Parameter "filt" can be:
0017 %  []: estimate from the data
0018 %  {B,A}: filter coefficients
0019 %  file name: mat file containg B and A
0020 %
0021 % Estimation from data might sometimes give erroneous results. Consider running
0022 % nt_destep on a clean step, save the filter estimates, and use them.
0023 %
0024 % See also nt_detrend_robust.
0025 %
0026 % NoiseTools.
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(nglitch); 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); % should contain B and A
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; % number of samples for IR fit
0051 B=1; A=1; % avoid error if no step
0052 %if nfit*2>guard; error('!'); end
0053 
0054 for iChan=1:size(x,2)
0055     % find step indices
0056     [splitI,score_vector,score]=nt_split(x(:,iChan),depth,thresh,guard);
0057     
0058     if ~isempty(splitI)
0059         splitI=splitI-5;
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)*m2]; % helps ensure stable filter
0068 %             step_response=y(splitI(iSplit)+(1:nfit),iChan);
0069 %             step_response=[step_response;ones(nfit,1)*m2]; % helps ensure stable filter
0070 
0071             %figure(1); clf; plot(step_response); pause
0072             
0073             if isempty(filt)  % estimate filter by fitting step response
0074                 %{
0075                 The filter is estimated anew for each step and each
0076                 channel.
0077                 %}
0078                 xx=ones(size(step_response,1),1)*step; yy=[step_response];
0079                 if 1
0080                     [B,A]=stmcb(diff([step_response]/step),10,10);
0081                 else
0082                     BETA0=[1,zeros(1,5),1,zeros(1,5)];
0083                     fun = @(beta,x)(filter([beta(1),beta(2),beta(3),beta(4),beta(5),beta(6)],...
0084                         [beta(7),beta(8),beta(9),beta(10),beta(11),beta(12)],x));
0085                     BETA = nlinfit(xx,yy,fun,BETA0);
0086                     BETA=BETA/BETA(7);
0087                     B=BETA(1:6); A=BETA(7:end);
0088                 end
0089             else
0090                 B=filt{1}; A=filt{2};
0091             end
0092                 
0093             y(splitI(iSplit)+1:end,iChan)=...
0094                 y(splitI(iSplit)+1:end,iChan)-...
0095                 filter(B,A,step*ones(size(y,1)-splitI(iSplit),1));
0096         end
0097     end
0098 end
0099 
0100 
0101 if ~nargout
0102     % don't return values, just plot
0103     disp(['steps at: ', num2str(splitI(2:end-1))]);
0104     figure(1); clf; nt_banner('nt_destep');
0105     xx=[x(:),y(:)];
0106     lim=[min(xx(:)),max(xx(:))]; lim=[lim(1)-0.1*(lim(2)-lim(1)), lim(2)+0.1*(lim(2)-lim(1))];
0107     subplot 211; plot([x,y]); xlabel('samples'); ylim(lim); legend('raw','clean'); legend boxoff
0108     subplot 212;     plot(y,'r'); xlabel('samples'); legend('clean'); legend boxoff
0109     figure(2); clf
0110     plot(filter(B,A,ones(200,1)),'.-'); xlabel('samples'); title('estimated step response of antialiasing filter');
0111     clear y, filt;
0112 end
0113 
0114 
0115 
0116 % test code
0117 if 0
0118     N=8;
0119     Wn=0.25; % CTF corner frequency is 1/4 of sampling rate
0120     nsamples=100;
0121     [B,A] = butter(N,Wn);
0122     y=filter(B,A,ones(nsamples,1));
0123     BETA0=[1, zeros(1,8), 1, zeros(1,8)];
0124     fun = @(beta,x)(filter([beta(1),beta(2),beta(3),beta(4),beta(5),beta(6),beta(7),beta(8),beta(9)],...
0125         [beta(10),beta(11),beta(12),beta(13),beta(14),beta(15),beta(16),beta(17),beta(18)],x));
0126     x=ones(nsamples,1);
0127     BETA = nlinfit(x,y,fun,BETA0);
0128     BETA=BETA/BETA(10);
0129     BB=BETA(1:9); AA=BETA(10:end);
0130     figure(1); clf;
0131     plot([y,filter(BB,AA,ones(size(y)))])
0132 end
0133 
0134 if 0
0135     N=8;
0136     Wn=0.25; % CTF corner frequency is 1/4 of sampling rate
0137     [B,A] = butter(N,Wn);
0138     x=[zeros(1100,1);ones(2000,1)];
0139     x=filter(B,A,x);
0140     x=x+0.01*randn(size(x));
0141     nt_destep(x);
0142 end
0143 if 0
0144     N=8;
0145     Wn=0.25; % CTF corner frequency is 1/4 of sampling rate
0146     [B,A] = butter(N,Wn);
0147     x=[zeros(1100,1);ones(1100,1);2*ones(1100,1);3*ones(1100,1);4*ones(1100,1);0*ones(1100,1)];
0148     x=filter(B,A,x);
0149     x=x+0.0*randn(size(x));
0150     nt_destep(x);
0151 end
0152 if 0
0153     x=ft_read_data('/data/meg/phantom090715_BrainampDBS_20150709_18.ds'); x=permute(x,[2,1,3]); x=x(:,33:306,:);
0154     x=x(:,100,1);
0155     nt_destep(x,[],[],30);
0156 end
0157 
0158

Generated on Mon 30-Jan-2017 18:59:11 by m2html © 2005