0001 function y=nt_inpaint(x,w,w1)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
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
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
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;
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
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,:));
0080 iGood=find(wDict(iPattern,:));
0081 if isempty(iBad); continue; end
0082 iInPattern=find(iBack==iPattern);
0083 iOutPattern=find(iBack~=iPattern);
0084
0085
0086 cc = c;
0087 [V,D]=eig(cc(iGood,iGood)); V=real(V); D=real(D);
0088 topcs=V(:,find(D/max(D) > PCA_THRESH));
0089 b=cc(iBad,iGood)*topcs/(topcs'*cc(iGood,iGood)*topcs);
0090 xx=(x(iInPattern,iGood));
0091 y(iInPattern,iBad)=xx*(topcs*b')-mn(iInPattern,iBad);
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
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
0139
0140
0141 [ww,~,IC,nOccurrences]=nt_unique(w,'rows');
0142 [nOccurrences,iSort]=sort(nOccurrences, 'descend');
0143 [~,iReverse]=sort(iSort);
0144 ww=ww(iSort,:);
0145 iBack=iReverse(IC);
0146
0147 function [w,x]=coalesceWeights(w,x)
0148
0149
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
0158 idx=find(~w1&w2&~w3);
0159 w2(idx)=0; w(2:end-1,:)=w2;
0160 idx=find(w1&~w2&w3);
0161 w2(idx)=1; w(2:end-1,:)=w2;
0162 x2(idx)=(x1(idx))+x3(idx)/2; x(2:end-1,:)=x2;
0163
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;
0172 else
0173 iFix=iIsolated+1;
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