www.pudn.com > image_mva_0.rar > simplisma.m

function [purspec,purint,purity_spec]=simplisma(data,varlist,offset,n,data2); 
% function [purspec,purint,purity_spec]=simplisma(data,varlist,offset,n,data2); 
% 
% It is a short non interactive version of SIMPLISMA taken from Windig's 
% article Chemometrics and Intelligent Laboratory Systems, 36, 1997, 3-16.  
% 
%   INPUT: 
%   data contains the data matrix (spectra in rows) 
%   data 2 can be ignored or empty. 
%       For second derivative applications data contains the conventional  
%       data and data2 contains the inverted 2nd data. 
%       to create data2 use function: 
%       data2=invder(data); 
%   Varlist contains the variable identifiers  
%   Offset is a correction factor for low intensity variables (1- no 
%   offset, 15 - large offset) 
%   n is a number of components 
%  
%   OUTPUT: 
%       purespec contains the pure spectra 
%       purint contains the intensities ('concentrations') of the pure 
%       spectra in the mixtures 
%       purity_spec - spectra containing purity spectra 
% 
% The programm will plot the purity and standard deviation spectra, where 
% the pure variables selected will be marked by a '*'. After each plot, 
% any key needs to be pressed to continue. 
 
%INITIALIZE; 
if nargin==5; 
   temp=data;data=data2;data2=temp;clear temp 
end 
[nspec,nvar]=size(data); purvarindex=[]; 
if nargin==4; 
   data2=[]; 
end; 
 
%CACULATE STATISTICS 
 
stddata=std(data)*sqrt(nspec-1)/sqrt(nspec); 
meandata=mean(data); 
meandataoffset=meandata+((offset/100)*max(meandata)); 
lengthdata=sqrt((stddata.*stddata+meandataoffset.*meandataoffset)*... 
   sqrt(nspec)); 
lengthmatrix=lengthdata(ones(1,nspec),:); 
datalengthscaled=data./lengthmatrix; 
puredata=stddata./meandataoffset; 
 
%DETERMINE PURE VARIABLES 
purity_spec=0*[1:nvar]; 
max_index=0; 
for i=1:n+1; 
   purvar=datalengthscaled(:,purvarindex); 
   for j=1:nvar; 
      addcolumn=datalengthscaled(:,j); 
      purvartest=[purvar addcolumn]; 
      matrix=purvartest'*purvartest; 
      weight(j)=det(matrix); 
   end; 
   purityspec=weight.*puredata; 
	purity_spec=[purity_spec; purityspec]; 
   maxindex=find(purityspec==max(purityspec)); 
   maxindex=maxindex(1); 
   figure(2) 
   subplot(3,2,1); plot(varlist,purityspec,'g',varlist(maxindex),... 
      purityspec(maxindex),'g*'); 
   max_index=[max_index, maxindex]; 
   axis([sort([varlist(1) varlist(length(varlist))]) 0 1.1*max(purityspec)]); 
   if varlist(1)>varlist(2); 
      set(gca,'Xdir','reverse'); 
   end; 
   title(['purity spectrum # ', num2str(i)]); 
   stdspec=weight.*stddata; 
   subplot(3,2,2);plot(varlist,  stdspec,'g',varlist(maxindex),... 
      stdspec(maxindex),'g*'); 
   axis([sort([varlist(1) varlist(length(varlist))]) 0 1.1*max(stdspec)]) ; 
   if varlist(1)>varlist(2); 
      set(gca,'Xdir','reverse'); 
   end; 
title(['standard deviation spectrum # ', num2str(i)]); 
    
   pause 
    
   purvarindex=[purvarindex maxindex]; 
end 
close(2) 
purvarindex(n+1)=[]; 
 
%RESOLVE SPECTRA 
 
purematrix=(data(:,purvarindex)); 
if isempty(data2) 
   purspec=purematrix\data; 
else; 
   purspec=purematrix\data2; 
end; 
 
%RESOLVE INTENSITIES 
 
if isempty(data2); 
   purint=data/purspec; 
else; 
   purint=data2/purspec; 
end; 
 
%SCALE 
 
if isempty(data2); 
   tsi=sum(data')'; 
else; 
   tsi=sum(data2')'; 
end; 
a=purint\tsi; 
purint=purint*diag(a); 
purspec=inv(diag(a))*purspec; 
H2.Position=[264 188 339 423]; 
figure(H2) 
subplot(2,1,1),plot(varlist,purspec), set(gca,'Xdir','reverse') 
title ('pure spectra') 
subplot(2,1,2), plot(purint), title ('pure intensity') 
H3.Position=[616 190 339 423]; 
figure(H3) 
for i=1:n+1; 
   subplot(n+1,1,i), plot(abs(varlist),purity_spec(i+1,:)) 
   hold on, plot(abs(varlist(max_index(i+1))),purity_spec(max_index(i+1)),'g*'); 
   set(gca,'Xdir','reverse') 
   hold off 
end