0001 function [y,idx,w]=nt_tsr(x,ref,shifts,wx,wref,keep,thresh)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 nt_greetings;
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047 if nargin<2; error('too few arguments'); end
0048 if nargin<3 || isempty(shifts); shifts=0; end
0049 if nargin<4; wx=[]; end
0050 if nargin<5; wref=[]; end
0051 if nargin<6 || isempty(keep); keep=[]; end
0052 if nargin<7 || isempty(thresh); thresh=10.^-20; end
0053
0054
0055 if size(x,1)~=size(ref,1); error('X and REF should have same nrows'); end
0056 if size(x,3)~=size(ref,3); error('X and REF should have same npages'); end
0057 if ~isempty(wx) && size(x,3)~=size(wx,3); error('X and WX should have same npages'); end
0058 if ~isempty(wx) && size(x,1)~=size(wx,1); error('X and WX should have same nrows'); end
0059 if ~isempty(wref) && size(ref,1)~=size(wref,1); error('REF and WREF should have same nrows'); end
0060 if ~isempty(wref) && size(ref,3)~=size(wref,3); error('REF and WREF should have same npages'); end
0061 if max(shifts)-min(0,min(shifts)) >= size(x,1); error('X has too few samples to support SHIFTS'); end
0062 if ~isempty(wx) && size(wx,2)~=1; error('wx should have ncols=1'); end
0063 if ~isempty(wref) && size(wref,2)~=1; error('wref should have ncols=1'); end
0064 if ~isempty(wx) && sum(wx(:))==0; error('weights on x are all zero!'); end
0065 if ~isempty(wref) && sum(wref(:))==0; error('weights on ref are all zero!'); end
0066
0067 if numel(shifts)>1000; error(['numel(shifts)=',num2str(numel(shifts)), ' (if OK comment out this line)']); end
0068
0069
0070
0071
0072
0073 offset1=max(0,-min(shifts));
0074 idx=1+offset1:size(x,1);
0075 x=x(idx,:,:);
0076 if ~isempty(wx); wx=wx(idx,:,:); end
0077 shifts=shifts+offset1;
0078
0079
0080 offset2=max(0,max(shifts));
0081 idx=1: size(ref,1)-offset2;
0082 x=x(idx,:,:);
0083 if ~isempty(wx); wx=wx(idx,:,:); end
0084
0085 [mx,nx,ox]=size(x);
0086 [mref,nref,oref]=size(ref);
0087
0088
0089 w=zeros([mx,1,oref]);
0090 if isempty(wx) && isempty(wref)
0091 w(1:mx,:,:)=1;
0092 elseif isempty(wref);
0093 w(:,:,:)=wx(:,:,:);
0094 elseif isempty(wx)
0095 for k=1:ox
0096 wr=wref(:,:,k);
0097 wr=ts_multishift(wr,shifts);
0098 wr=min(wr,[],2);
0099 w(:,:,k)=wr;
0100 end;
0101 else
0102 for k=1:ox
0103 wr=wref(:,:,k);
0104 wr=nt_multishift(wr,shifts);
0105 wr=min(wr,[],2);
0106 wr=min(wr,wx(1:size(wr,1),:,k));
0107 w(:,:,k)=wr;
0108 end
0109 end
0110 wx=w;
0111 wref=zeros(mref,1,oref);
0112 wref(idx,:,:)=w;
0113
0114
0115 [x,mn1]=nt_demean(x,wx);
0116 ref=nt_demean(ref,wref);
0117
0118
0119 ref=nt_normcol(ref,wref);
0120 ref=nt_pca(ref,0,[],10^-6);
0121 ref=nt_normcol(ref,wref);
0122
0123
0124 [cref,twcref]=nt_cov(ref,shifts,wref);
0125 [cxref,twcxref]=nt_xcov(x,ref,shifts,wx);
0126
0127
0128 r=nt_regcov(cxref/twcxref,cref/twcref,keep,thresh);
0129
0130
0131
0132
0133 y=zeros(mx,nx,ox);
0134 for k=1:ox
0135 z=nt_multishift(ref(:,:,k),shifts)*r;
0136
0137 y(:,:,k)=x(1:size(z,1),:,k)-z;
0138 end
0139 [y,mn2]=nt_demean(y,wx);
0140
0141
0142 idx=1+offset1:size(y,1)+offset1;
0143 mn=mn1+mn2;
0144 w=wref;
0145
0146