Home > NoiseTools > nt_destep.m

# nt_destep

## PURPOSE

[y,stepList]=nt_destep(x,thresh,guard,depth,minstep) - remove step glitch from MEG data

## SYNOPSIS

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

## DESCRIPTION

    [y,stepList]=nt_destep(x,thresh,guard,depth,minstep) - remove step glitch from MEG data

      y: step-removed data
      stepList: indices of steps

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

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

     See also nt_detrend, nt_deglitch, nt_split.

     NoiseTools.

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_banner h=nt_banner(text,varargin) - annotate with text at head of figure
- nt_demean [y,mn]=nt_demean(x,w) - remove weighted mean over cols
- nt_destep [y,stepList]=nt_destep(x,thresh,guard,depth,minstep) - remove step glitch from MEG data
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_split [idx,score_vector,score]=nt_split(x,depth,thresh,guard,minstep) - split time series into intervals

This function is called by:

- nt_destep [y,stepList]=nt_destep(x,thresh,guard,depth,minstep) - remove step glitch from MEG data

## SOURCE CODE

    0001 function [y,stepList]=nt_destep(x,thresh,guard,depth,minstep);
    0002 %[y,stepList]=nt_destep(x,thresh,guard,depth,minstep) - remove step glitch from MEG data
    0003 %
    0004 %  y: step-removed data
    0005 %  stepList: indices of steps
    0006 %
    0007 %  x: data to clean (time * channels)
    0008 %  thresh: threshold for variance reduction [default: 0.1]
    0009 %  guard: minimum duration of stable interval in samples [default: 1000]
    0010 %  depth: recursion depth for nt_split [default:6], determines number of steps
    0011 %  minstep: if step size smaller skip step
    0012 %
    0013 % Searches for large steps, removes them if variance ratio less than
    0014 % "thresh", and both intervals shorter than "guard".
    0015 %
    0016 % See also nt_detrend, nt_deglitch, nt_split.
    0017 %
    0018 % NoiseTools.
    0019 nt_greetings;
    0020
    0021 if nargin<1; error('!'); end
    0022 if nargin<2||isempty(thresh); thresh=0.1; end
    0023 if nargin<3||isempty(guard); guard=100; end
    0024 if nargin<4 || isempty(depth); depth=6; end
    0025 if nargin<5 ; minstep=[]; end
    0026
    0027 if isempty(minstep); minstep=(max(x(:))-min(x(:)))*0.0000001; end
    0028 y=x;
    0029
    0030 %disp(numel(thresh))
    0031
    0032 for iChan=1:size(x,2)
    0033     % find step indices
    0034     [stepList,score_vector,score]=nt_split(x(:,iChan),depth,thresh,guard,minstep);
    0035
    0036     if ~isempty(stepList)
    0037         stepList=[0,stepList,size(x,1)];
    0038         for iSplit=2:numel(stepList)-1
    0039             y1=y(stepList(iSplit-1)+1:stepList(iSplit),iChan); % plateau before
    0040             y2=y(stepList(iSplit)+1:stepList(iSplit+1),iChan); % plateau after
    0041             step=(mean(y2)-mean(y1));
    0042             y(stepList(iSplit)+1:end,iChan)=y(stepList(iSplit)+1:end,iChan)-step;
    0043         end
    0044     end
    0045 end
    0046 stepList=stepList(2:end-1);
    0047
    0048
    0049 if ~nargout
    0050     % don't return values, just plot
    0051     disp(['steps at: ', num2str(stepList(2:end-1))]);
    0052     figure(1); clf; nt_banner('nt_destep');
    0053     xx=[x(:),y(:)];
    0054     lim=[min(xx(:)),max(xx(:))]; lim=[lim(1)-0.1*(lim(2)-lim(1)), lim(2)+0.1*(lim(2)-lim(1))];
    0055     subplot 211; plot([x,y]); xlabel('samples'); ylim(lim); legend('raw','clean'); legend boxoff
    0056     subplot 212;     plot(y,'r'); xlabel('samples'); legend('clean'); legend boxoff
    0057     clear y, stepList;
    0058 end
    0059
    0060
    0061
    0062 % test code
    0063 % if 0
    0064 %     N=8;
    0065 %     Wn=0.25; % CTF corner frequency is 1/4 of sampling rate
    0066 %     nsamples=100;
    0067 %     [B,A] = butter(N,Wn);
    0068 %     y=filter(B,A,ones(nsamples,1));
    0069 %     BETA0=[1, zeros(1,8), 1, zeros(1,8)];
    0070 %     fun = @(beta,x)(filter([beta(1),beta(2),beta(3),beta(4),beta(5),beta(6),beta(7),beta(8),beta(9)],...
    0071 %         [beta(10),beta(11),beta(12),beta(13),beta(14),beta(15),beta(16),beta(17),beta(18)],x));
    0072 %     x=ones(nsamples,1);
    0073 %     BETA = nlinfit(x,y,fun,BETA0);
    0074 %     BETA=BETA/BETA(10);
    0075 %     BB=BETA(1:9); AA=BETA(10:end);
    0076 %     figure(1); clf;
    0077 %     plot([y,filter(BB,AA,ones(size(y)))])
    0078 % end
    0079 if 0
    0080     N=8;
    0081     Wn=0.25; % CTF corner frequency is 1/4 of sampling rate
    0082     [B,A] = butter(N,Wn);
    0083     x=[zeros(1100,1);ones(2000,1)];
    0084     x=filter(B,A,x);
    0085     x=x+0.01*randn(size(x));
    0086     [y,stepList]=nt_destep(x);
    0087     figure(1); clf;
    0088     subplot 211;
    0089     plot([x,y]); legend('raw','destep');legend boxoff
    0090     subplot 212;
    0091     plot([y,nt_deglitch(y,stepList)]);
    0092     legend('destep','deglitch');legend boxoff
    0093 end
    0094 if 0
    0095     N=8;
    0096     Wn=0.25; % CTF corner frequency is 1/4 of sampling rate
    0097     [B,A] = butter(N,Wn);
    0098     x=[zeros(1000,1);ones(1000,1);2*ones(1000,1);3*ones(1000,1);4*ones(1000,1);0*ones(1000,1)];
    0099     x=filter(B,A,x);
    0100     x=x+0.0*randn(size(x));
    0101     [y,stepList]=nt_destep(x);
    0102     figure(1); clf;
    0103     subplot 211;
    0104     plot([x,y]); legend('raw','destep'); legend boxoff
    0105     subplot 212;
    0106     plot([y,nt_deglitch(y,stepList)]);
    0107     legend('destep','deglitch');legend boxoff
    0108 end
    0109 if 0
    0110     x=ft_read_data('/data/meg/litvak/phantom090715_BrainampDBS_20150709_18.ds'); x=permute(x,[2,1,3]); x=x(:,33:306,:);
    0111     x=nt_demean(x(1:2300,100,1));
    0112     x=x*10^12; % -->pT
    0113     [y,stepList]=nt_destep(x);
    0114     figure(1); clf;
    0115     subplot 211;
    0116     plot([x,y]); legend('raw','destep'); legend boxoff
    0117     subplot 212;
    0118     plot([y,nt_deglitch(y,stepList)]);
    0119     legend('destep','deglitch');legend boxoff
    0120 end
    0121
    0122