0001 function [todss,pwr0,pwr1]=nt_dss0(c0,c1,keep1,keep2)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 nt_greetings;
0015
0016 if nargin<4||isempty(keep2); keep2=10.^-9; end
0017 if nargin<3; keep1=[]; end
0018 if nargin<2; error('needs at least two arguments'); end
0019
0020 if size(c0)~=size(c1); error('C0 and C1 should have same size'); end
0021 if size(c0,1)~=size(c0,2); error('C0 should be square'); end
0022
0023 if any(find(isnan(c0)))
0024 error('NaN in c0');
0025 end
0026 if any(find(isnan(c1)))
0027 error('NaN in c1');
0028 end
0029 if any(find(isinf(c0)))
0030 error('INF in c0');
0031 end
0032 if any(find(isinf(c1)))
0033 error('INF in c1');
0034 end
0035
0036 [topcs1,evs1]=nt_pcarot(c0,keep1,keep2);
0037 evs1=abs(evs1);
0038
0039
0040 if ~isempty(keep1); topcs1=topcs1(:,1:keep1); evs1=evs1(1:keep1); end
0041 if ~isempty(keep2); idx=find(evs1/max(evs1)>keep2); topcs1=topcs1(:,idx); evs1=evs1(idx); end
0042
0043
0044 N=diag(sqrt(1./(evs1)));
0045 c2=N'*topcs1'*c1*topcs1*N;
0046
0047
0048 [topcs2,evs2]=nt_pcarot(c2,keep1,keep2);
0049
0050
0051 todss=topcs1*N*topcs2;
0052 N2=diag(todss'*c0*todss);
0053 todss=todss*diag(1./sqrt(N2));
0054
0055
0056
0057 pwr0=sqrt(sum((c0'*todss).^2));
0058 pwr1=sqrt(sum((c1'*todss).^2));
0059
0060 if ~nargout
0061 figure(100); clf
0062 plot(pwr1./pwr0, '.-');
0063 xlabel('component');
0064 ylabel('score');
0065 clear todss pwr0 pwr1
0066 end
0067
0068