Answers to the problems in the Tony Davies column in Spectroscopy Europe 13/6 (2001)
The five lines of code
n=size(Tf,2);
TCc=zeros(1,n);
for
k=1:n
TCc(1,k)=(Tf(1,k)-32)*5/9;
end
may be replaced by the single line
TCc=(Tf-32)*5/9;
The operation of subtracting 32 and multiplying by 5/9 is automatically applied to each element of the whole vector Tf when it is written like this. It would apply to all the elements of the matrix were Tf a matrix.
We need the line n=size(X,1); but all the subsequent lines may be replaced by
Xd=X-ones(n,1)*mean(X);
The result of mean(X) is a row vector, with the same length as the rows of X, containing the means of the columns of X, which is the mean spectrum in this example (see the comments in the code for example 3 in the original article). Since ones(n,1) gives us as column vector of ones, of the same size as the columns of X, the product ones(n,1)*mean(X) is a matrix of the same size as X with the mean spectrum in each row. Subtracting this from X centres it.
Here is a rather more efficient version.
| % multiplicative scatter correction | |
| load('S12a'); | |
| X=spec; | |
| [r,c]=size(X); | % find the size of the X matrix |
| msp=mean(X); | % as in the previous example this is the mean spectrum |
| Gm=mean(msp); | % and this is the mean of all the values in the mean spectrum |
| cmsp=msp-Gm; | % center the mean spectrum |
| div=cmsp*cmsp'; | % sum of squares about the mean for the values in msp |
| B=(X*cmsp')/div; | % vector of r regression coefficients, one for each spectrum |
| A=mean(X')'-B*Gm; | % vector of r intercepts |
| Xmsc=(X-A*ones(1,c))./(B*ones(1,c)); | % do the MSC |
Most of this is very similar to the litmus example in the article. What you may not have seen before is the ./ operator on the final line. The two matrices before and after this operator are both of dimension r by c, and the ./ produces a result that is an element by element division of the first matrix by the second.