Home > NoiseTools > nt_inpaint.m

nt_inpaint

PURPOSE ^

y=nt_inpaint(x,w,cw1) - interpolate missing samples of matrix

SYNOPSIS ^

function y=nt_inpaint(x,w,w1)

DESCRIPTION ^

y=nt_inpaint(x,w,cw1) - interpolate missing samples of matrix

  y: interpolated data

  x: data to interpolate (time*channels)
  w: weight matrix to control interpolation (same size as x)
  w1: weight matrix to calculate covariance (if missing use min(w,2))

 The missing parts of each channel x(:,i) (flagged by w(:,i)) are replaced
 from intact channels, using projection matrix calculated from w1.

  NoiseTools
 See nt_regw, nt_detrend_robust, nt_star2.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function y=nt_inpaint(x,w,w1)
0002 %y=nt_inpaint(x,w,cw1) - interpolate missing samples of matrix
0003 %
0004 %  y: interpolated data
0005 %
0006 %  x: data to interpolate (time*channels)
0007 %  w: weight matrix to control interpolation (same size as x)
0008 %  w1: weight matrix to calculate covariance (if missing use min(w,2))
0009 %
0010 % The missing parts of each channel x(:,i) (flagged by w(:,i)) are replaced
0011 % from intact channels, using projection matrix calculated from w1.
0012 %
0013 %  NoiseTools
0014 % See nt_regw, nt_detrend_robust, nt_star2.
0015 %
0016 
0017 nt_greetings;
0018 PCA_THRESH=0.000000001;
0019 
0020 if nargin<2; error('!'); end
0021 if nargin<3; w1=[]; end
0022 
0023 if nargout==0;
0024     % don't return result, just plot
0025     y=nt_inpaint(x,w,c);
0026     figure 1; clf;
0027     subplot 211; 
0028     plot(nt_unfold(x)); title('raw');
0029     subplot 212; 
0030     plot(nt_unfold(y)); title('processed');
0031     clear y; return
0032 end
0033 
0034 %% check/fix sizes
0035 if any(size(w) ~= size(x)); error('w must be of same size as x'); end
0036 if any(w ~= 0 & w ~= 1); error('weights must be 0 or 1'); end
0037 m=size(x,1);
0038 x=nt_unfold(x);
0039 w=nt_unfold(w);
0040 iSkipChannels=find(all(w==0));
0041 if ~isempty (iSkipChannels); 
0042     disp(['ignore channels marked as bad: ', num2str(iSkipChannels)]);
0043     x0=x;
0044     x(:,iSkipChannels)=[];
0045     w(:,iSkipChannels)=[];
0046 end
0047 if ~isempty(w1)
0048     if size(w1,2)~=1; 
0049         disp(size(w1,2)); 
0050         error('!'); 
0051     end
0052 else
0053     disp('calculating covariance based on weight matrix');
0054     w1=min(w,[],2);
0055     propGood=mean(w1);
0056     disp(['proportion good samples: ',num2str(propGood)]);
0057     if ~propGood;
0058         error('no good samples')
0059     end
0060 end
0061 y=x;
0062 x=nt_demean(x,w1);
0063 mn=x-y; % save mean
0064 s=sum(bsxfun(@times, x, w1));
0065 c=nt_cov(bsxfun(@times, x, w1))/sum(w1);
0066 disp([sum(w1), numel(find(isnan(c)))]);
0067 
0068 
0069 %% find different patterns of weights, merge to reduce computation
0070 [wDict,nOccurrences,iBack]=patternDict(w); 
0071 disp(['distinct patterns, mean weight: ', num2str(size(wDict,1)), ', ', num2str(mean(w(:)))]);
0072 [w,x]=coalesceWeights(w,x);
0073 [w,x]=coalesceWeights(w,x);
0074 [w,x]=coalesceWeights(w,x);
0075 [wDict,nOccurrences,iBack]=patternDict(w); 
0076 disp([' ... reduced to: ', num2str(size(wDict,1)), ', ', num2str(mean(w(:)))]);
0077 
0078 for iPattern=1:size(wDict,1);
0079     iBad=find(~wDict(iPattern,:)); % channel indices
0080     iGood=find(wDict(iPattern,:)); 
0081     if isempty(iBad); continue; end
0082     iInPattern=find(iBack==iPattern); % sample indices
0083     iOutPattern=find(iBack~=iPattern);
0084     %ss=s-sum(bsxfun(@times, x(iInPattern,:), w1(iInPattern,:)));
0085     %mm=ss/numel(iOutPattern);
0086     cc = c;% - mm'*mm;%*numel(iOutPattern);
0087     [V,D]=eig(cc(iGood,iGood)); V=real(V); D=real(D);
0088     topcs=V(:,find(D/max(D) > PCA_THRESH)); % discard weak dims
0089     b=cc(iBad,iGood)*topcs/(topcs'*cc(iGood,iGood)*topcs); % projection matrix
0090     xx=(x(iInPattern,iGood));
0091     y(iInPattern,iBad)=xx*(topcs*b')-mn(iInPattern,iBad); % projection
0092 end
0093 
0094 if ~isempty(iSkipChannels)
0095     yy=x0;
0096     yy(:,setdiff(1:size(x0,2),iSkipChannels))=y;
0097     y=yy;
0098 end
0099 
0100 %%% TEST %%%
0101 if 0
0102     x0=sin(2*pi*(1:10000)'*(1:10)/10000);
0103     x=x0*randn(10,20)+10;
0104     w=ones(size(x));  
0105     w(1:4000,1)=0; w(5001:6000,3)=0; w(6000:6500,4)=0; w(7000:7500,5)=0; w(8000:8500,6)=0; w(9000:9500,7)=0;
0106     x(6000:6500,4)=0;x(7000:7500,5)=10;
0107     y=nt_inpaint(x,w);
0108     figure(2); clf
0109     subplot 211; plot([x(:,1),y(:,1)]); legend('original','processed'); legend boxoff
0110     subplot 212; plot([x(:,1)-y(:,1)]); legend('difference'); legend boxoff
0111 end
0112 if 0
0113     x0=sin(2*pi*(1:10000)'*(1:10)/10000);
0114     x=x0*randn(10,20)+10;
0115     w=ones(size(x));  
0116     w(1:4000,1)=0; 
0117     w(:,8)=0; 
0118     x(1:2000,1)=1;
0119     y=nt_inpaint(x,w);
0120     figure(2); clf
0121     subplot 211; plot([x(:,1),y(:,1)]); legend('original','processed'); legend boxoff
0122     subplot 212; plot([x(:,1)-y(:,1)]); legend('difference'); legend boxoff
0123 end
0124 if 0
0125     [p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/mc/MC_aespa_speech_47.mat');
0126     x=x'; x=x(:,1:128); x=x(1:10000,:);
0127     x=nt_demean(x);
0128     [x,w]=nt_detrend_robust(x,10);
0129     y=nt_inpaint(x,w);
0130     figure(1); clf
0131     subplot 311; plot(x); subplot 312; plot(y); subplot 313; plot(x-y);
0132     figure(2); clf
0133     ch=35;subplot 311; plot([x(:,ch),y(:,ch)]); subplot 312; plot([x(:,ch)-y(:,ch)]); subplot 313; plot(w(:,ch), '.-');
0134 end
0135 
0136 
0137 function [ww,nOccurrences,iBack]=patternDict(w);
0138 % ww: dictionary of patterns
0139 % nOccurrences: number of times each pattern occurred
0140 % iBack: index to reconstruct input from dictionary
0141 [ww,~,IC,nOccurrences]=nt_unique(w,'rows');
0142 [nOccurrences,iSort]=sort(nOccurrences, 'descend'); % sort by decreasing number
0143 [~,iReverse]=sort(iSort); %
0144 ww=ww(iSort,:); % same order for patterns, w = ww(iReverse1(IC),:)
0145 iBack=iReverse(IC); % w = ww(iBack,:)
0146 
0147 function [w,x]=coalesceWeights(w,x)
0148 % reduce the number of weight patterns based on pruning & interpolation, to
0149 % save time
0150 nchans=size(w,2);
0151 w1=w(1:end-2,:); 
0152 w2=w(2:end-1,:); 
0153 w3=w(3:end,:); 
0154 x1=x(1:end-2,:); 
0155 x2=x(2:end-1,:); 
0156 x3=x(3:end,:); 
0157 % weights that correspond to a single sample
0158 idx=find(~w1&w2&~w3); % isolated good, set to bad
0159 w2(idx)=0; w(2:end-1,:)=w2;
0160 idx=find(w1&~w2&w3); % isolated bad, interpolate, set to good
0161 w2(idx)=1; w(2:end-1,:)=w2; 
0162 x2(idx)=(x1(idx))+x3(idx)/2; x(2:end-1,:)=x2;
0163 % weight patterns that occur rarely
0164 [wDict,nOccurrences,iBack]=patternDict(w); 
0165 idx=find(nOccurrences<=3);
0166 isolated=iBack(idx);
0167 for k=1:numel(isolated)
0168     iIsolated=isolated(k);
0169     if iIsolated-1>=1 && iIsolated+1<=size(w,1)            
0170         if sum(abs(w(iIsolated-1,:)-w(iIsolated,:)))<sum(abs(w(iIsolated,:)-w(iIsolated+1,:)));
0171             iFix=iIsolated-1;% preceding is closest
0172         else
0173             iFix=iIsolated+1;% following is closest
0174         end
0175         idx=find(w(iIsolated,:)&~w(iFix,:));
0176         w(iIsolated,idx)=0;
0177         idx=find(~w(iIsolated,:)&w(iFix,:));
0178         w(iIsolated,idx)=1;
0179         x(iIsolated,idx)=(x(iIsolated-1,idx)+x(iIsolated+1,idx))/2;
0180     end
0181 end
0182

Generated on Sun 11-Dec-2016 18:36:17 by m2html © 2005