nt_eyeblink

[y,z,mask]=nt_eyeblink(x,eyechans,nremove) - project out eyeblinks

  y: clean data
  z: eyeblink components
  mask: mask used to weight eyeblink intervals

  x: data to clean (time X channels)
  eyechans: ocular channel numbers, or eog signals
  nremove: number of components to remove

 NoiseTools

0001 function [y,z,mask]=nt_eyeblink(x,eyechans,nremove,sr)
0002 %[y,z,mask]=nt_eyeblink(x,eyechans,nremove) - project out eyeblinks
0003 %
0004 %  y: clean data
0005 %  z: eyeblink components
0006 %  mask: mask used to weight eyeblink intervals
0007 %
0008 %  x: data to clean (time X channels)
0009 %  eyechans: ocular channel numbers, or eog signals
0010 %  nremove: number of components to remove
0011 %
0012 % NoiseTools
0013 nt_greetings;
0014 
0015 if nargin<2; error('!'); end
0016 if nargin<3||isempty(nremove); nremove=1; end
0017 if nargin<4; sr=[]; end
0018 
0019 if nargout==0
0020     % plot results
0021     [y,z,mask]=nt_eyeblink(x,eyechans,nremove,sr);
0022     figure(201); clf;
0023     subplot 211; plot(x); title('raw');
0024     subplot 212; plot(y); title('clean');
0025     figure(202); clf;
0026     subplot 211; plot(mask); title('mask');
0027     subplot 212; plot(z(:,1:3)); title('eyeblink components 1:3');
0028     disp(size(z)); 
0029     disp(size(kurtosis(z)));
0030     figure(203); clf;
0031     plot(kurtosis(z));
0032     drawnow
0033     clear y z mask
0034 end
0035 
0036 if size(eyechans,1) ~= size(x,1)
0037     % interpret as channel numbers
0038     eyechans=x(:,eyechans);
0039     if isempty(sr); error ('!'); end
0040 end
0041 
0042 % demean, detrend, highpass eyechans to emphasize blinks
0043 eyechans=nt_demean(eyechans);
0044 ORDER=1;
0045 eyechans=nt_detrend(eyechans,ORDER);
0046 if sr
0047     HPF=1;
0048     [b,a]=butter(2,HPF/(sr/2),'high');
0049     eyechans=filtfilt(b,a,eyechans);
0050 end
0051 eyechans=nt_normcol(eyechans);
0052 
0053 % PCA to merge across eye channels
0054 topcs=nt_pca0(nt_normcol(eyechans));
0055 eyechans=eyechans*topcs(:,1);
0056 
0057 % mask to emphasize eyeblink intervals
0058 mask=mean(eyechans.^2,2);
0059 quantile=.8;
0060 tmp=sort(mask);
0061 mask=min(mask,tmp(round(size(mask,1)*quantile))); % avoid extreme weight
0062 
0063 % DSS to emphasize eyeblink components
0064 topcs=nt_pca0(x);
0065 xx=x*topcs(:,1:min(10,size(x,2)));
0066 c0=nt_cov(xx);
0067 c1=nt_cov(nt_demean(bsxfun(@times,xx,mask)));
0068 [todss,pwr0,pwr1]=nt_dss0(c0,c1); 
0069 %figure(99); clf; plot(pwr1./pwr0,'.-'); title('nt_eyeblink'); ylabel('score'); xlabel('component');
0070 z=nt_mmat(xx,todss(:,1:nremove));
0071 
0072 % figure(10); clf; subplot 211; plot(eyechans);
0073 % subplot 212; plot(z);
0074 % pause
0075 
0076 y=nt_tsr(x,z);
0077 
0078 % figure(11); clf; subplot 211; plot(x);
0079 % subplot 212; plot(nt_tsr(x,y));
0080 % pause
0081 
0082 
0083 
0084 
0085 
0086 
0087 
0088 
0089