Home > NoiseTools > nt_deglitch.m

# nt_deglitch

## PURPOSE

[y,B,A]=nt_deglitch(x,glitchList,filt,nfit) - remove temporally localized glitch

## SYNOPSIS

function [y,B,A]=nt_deglitch(x,glitchList,nfilt,nfit,advance);

## DESCRIPTION

    [y,B,A]=nt_deglitch(x,glitchList,filt,nfit) - remove temporally localized glitch

      y: clean data
      B,A: filter coefficients for last estimated event

      x: data to clean (time * channels)
      glitchList: list of glitch positions (cell array if x multichannel)
      nfilt: number of filter coefficients to estimate [default: 10, 10]
      nfit: number of samples over which to fit filter [default: 50]
      advance: number of samples to include before event [default: 3]

     If x is multichannel glitchI should be a cell array.

     See also nt_detrend.

     NoiseTools.

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_banner h=nt_banner(text,varargin) - annotate with text at head of figure
- nt_deglitch [y,B,A]=nt_deglitch(x,glitchList,filt,nfit) - remove temporally localized glitch
- 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

This function is called by:

- nt_deglitch [y,B,A]=nt_deglitch(x,glitchList,filt,nfit) - remove temporally localized glitch
- nt_destep [y,stepList]=nt_destep(x,thresh,guard,depth,minstep) - remove step glitch from MEG data

## SOURCE CODE

    0001 function [y,B,A]=nt_deglitch(x,glitchList,nfilt,nfit,advance);
    0002 %[y,B,A]=nt_deglitch(x,glitchList,filt,nfit) - remove temporally localized glitch
    0003 %
    0004 %  y: clean data
    0005 %  B,A: filter coefficients for last estimated event
    0006 %
    0007 %  x: data to clean (time * channels)
    0008 %  glitchList: list of glitch positions (cell array if x multichannel)
    0009 %  nfilt: number of filter coefficients to estimate [default: 10, 10]
    0010 %  nfit: number of samples over which to fit filter [default: 50]
    0011 %  advance: number of samples to include before event [default: 3]
    0012 %
    0013 % If x is multichannel glitchI should be a cell array.
    0014 %
    0015 % See also nt_detrend.
    0016 %
    0017 % NoiseTools.
    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         %The filter is estimated anew for each event.
    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     % don't return values, just plot
    0062     %disp(['steps at: ', num2str(stepI(2:end-1))]);
    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 % test code
    0076 if 0
    0077     N=8;
    0078     Wn=0.25; % CTF corner frequency is 1/4 of sampling rate
    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; % CTF corner frequency is 1/4 of sampling rate
    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; % CTF corner frequency is 1/4 of sampling rate
    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     %x=filter(B,A,x);
    0113     x=x+0.0*randn(size(x));
    0114     nt_destep(x);
    0115 end
    0116 if 0
    0117     % this example needs fixing
    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