matlabpcaspectrummixing

Matlab: Principle component analysis on signal (spectral unmixing)


I have a spectrum (wavelength(x) versus absorption(y)) which is a mix of two signals (xa,ya) and (xb,yb). I am trying to use PCA (code I found online) to unmix the signals in (x,y):

%step 1, input data
numdata=length(data(:,1));
x=data(:,1);
y=data(:,1);

%step 2, finding a mean and subtracting
xmean=mean(x);
ymean=mean(y);

xnew=x-xmean*ones(numdata,1);
ynew=y-ymean*ones(numdata,1);

subplot(3,1,1);
plot(x,y, 'o');
title('Original Data');

%step 3, covariance matrix
covariancematrix=cov(xnew,ynew);

%step 4, Finding Eigenvectors
[V,D] = eig(covariancematrix);
D=diag(D);
maxeigval=V(:,find(D==max(D)));


%step 5, Deriving the new data set
%finding the projection onto the eigenvectors

finaldata=maxeigval'*[xnew,ynew]';
subplot(3,1,2);
stem(finaldata, 'DisplayName', 'finaldata', 'YDataSource', 'finaldata');
title('PCA 1D output ')
%we do a classification now
subplot(3,1,3);
title('Final Classification')
hold on
for i=1:size(finaldata,2)
    if  finaldata(i)>=0
        plot(x(i),y(i),'o')
        plot(x(i),y(i),'r*')

    else
        plot(x(i),y(i),'o')
        plot(x(i),y(i),'g*')
    end

end

How best to apply the PCA output to unmix (y) into components ya and yb? I don't have experience with PCA and cannot find any good tutorials online for this application. Is it best to generate the eigenvectors for training spectrum and then compare with test spectrum? Thanks

enter image description here


Solution

  • Section 3.3 of this thesis is informative: https://brage.bibsys.no/xmlui//bitstream/handle/11250/2371385/12296_FULLTEXT.pdf?sequence=1&isAllowed=y

    "PCA itself is not a classification method but this was done based on the knowledge of which absorption spectra belonged to which material. PCA can, however, be used as a tool in classification. In order to do so, one needs training data. PCA is performed on the training data, and and some test data is projected on to the basis of the training data. This is called PCA decomposition."

    So I think I can use above code as a starting point.