function [CenterX,MahalaD,MahalaDt,Xt]=PCAMahala(X,Xa,e)
% X (m x n)為預處理之后的校正集光譜矩陣(row is samples,column is spectrum)
% Xa(p x n)為預處理之后校正集樣品之外的所有樣品的光譜矩陣(row is samples,column is spectrum)
% Xt(l x n)為預處理之后經過馬氏距離篩選的預測集光譜矩陣(row is samples,column is spectrum)
% MahalaD(m x 1)為校正集樣品的馬氏距離
% MahalaDt(p x 1)為校正集樣品之外的所有樣品的光譜矩陣
% l為經過馬氏距離篩選出來的樣品數
% e(設定設置閾值范圍)馬氏距離標準差的權重系數
[m,n]=size(X);
[p,n]=size(Xa);
Xt=zeros(p,n);
Xm=mean(X);
CenterX=X-ones(m,1)*Xm;
M=CenterX'*CenterX./(m-1);
Minv=inv(M);
% (計算校正集樣品的馬氏距離)
for i=1:m
    MahalaD(i)=sqrt((X(i,:)-Xm)*Minv*(X(i,:)-Xm)');
end
MahalaDm=mean(MahalaD);
Dstd=std(MahalaD);
l=1;
for i=1:p
    MahalaDt(i)=sqrt((Xa(i,:)-Xm)*Minv*(Xa(i,:)-Xm)'); %(計算預測集的馬氏距離)
    if MahalaDt(i)<(MahalaDm+e*Dstd)                   %(篩選預測集樣品)
        Xt(l,:)=Xa(i,:);
        l=l+1;
    end
end
for i=p:-1:l
    Xt(i,:)=[];
end