matlab hsv提取s_Matlab预处理高光谱和近外光谱（食品医药数据）

1. 微分处理光谱

一阶微分公式：

x_(i,1st)=(x_i-x_(i+g))/g

二阶微分公式：

x_(i,2st)=(x_i+x_(i+g)-2x_i)/g^2

式中，x_i是谱图数据中i波数下的透过率，g为窗口宽度

一阶导数MATLAB代码：

X1st=diff(X,1);%X 为输入光谱矩阵，X1st 为输出一阶导数光谱矩阵

二阶导数MATLAB代码：

X2st=diff(X,2); %X 为输入光谱矩阵，X2st 为输出二阶导数光谱矩阵

2. 多元散射校正处理谱图

（1）计算平均光谱：

x ̅_(i,j)=(∑_(i=1)^n▒X_(i,j) )/n

（2）一元线性回归：

X_i=m_i X ̅+b_i

（3）多元散射校正：

X_i(MSC) =((X_i-b_i))/m_i

公式中X表示n×p维定标光谱数据矩阵，n为样品数，p为波点数。X ̅表示所有样品的原始近红外光谱在各个波长点处求平均值所得到的平均光谱矢量，X_i是1×p维矩阵，表示单个样品光谱矢量，m_i和b_i分别表示各样品近红外光谱X_i与平均光谱X进行一元线性回归后得到的相对偏移系数和平移变量。I代表第i个样本，j代表第j个波数。

MATLAB代码如下：

me=mean(X);

[m,~]=size(X);

for i=1:m,

p=polyfit(me,X(i,:),1);

Xmsc(i,:)=(X(i,:)- p(2)*ones(1,n))./(p(1)*ones(1,n));

End

代码中输入光谱矩阵为 X，输出光谱矩阵为 Xmsc。

3. 标准正态变量变换（SNV）

X_(i,SNV)=(x_(i,k)-x_i)/√((∑_(k=i)^m▒〖(x_(i,k)-x_i)〗^2 )/((m-1)))

式中x_i是第i样品光谱的平均值，k=1,2,…,m。m为波长点数；i=1,2,..,n。n为校正样品数；X_(i,SNV)是变换后的光谱。

SNV的MATLAB代码如下：

[~,n]=size(X);

rmean=mean(X,2);

dr=X-repmat(rmean,1,n);

Xsnv=dr./repmat(sqrt(sum(dr.^2,2)/(n-1)),1,n)

代码中 X 为光谱输入矩阵，Xsnv 为 SNV 处理后光谱输出矩阵。

4. 小波变换处理谱图

傅里叶变换是小波变换的基础，傅里叶变换是一种分析信号的方法。傅里叶变换的基本公式为：

F(ω)=∫_(-∞)^(+∞)▒〖f(t)e^(-iωt) dt〗

公式中 f (t )是输入信号,在这里可看作光谱信号，t是时域信号，可以看作是波数，F (ω)是在信号频谱。

傅里叶变换原理表明：任何连续测量的时序或信号，都可以表示为不同频率正弦信号的无限叠加。傅里叶变换将原来难以处理的时域信号转换成了易于分析的频域信号，在频域表示的信号图称为信号的频谱。在频谱中可以看到时域谱图中看不到的特征，处理这些特征后利用傅里叶逆变换可以将这些频域信号转换成时域信号。

小波变换是傅里叶变换的延伸。傅里叶变换有很多缺点：不能求解变系数线性微分方程，只适用于平移不变的线性系统，不能反映信号在局部时间范围内的信息。小波变换继承了傅里叶变换的长处，同时解决了傅里叶变换的一些缺点。利用小波变换能够“放大”信号局部时间、频率范围内的信息，消除信号噪声。连续小波变换基本公式为：

Υ(s,τ)=∫▒〖f(t)φ_(s,τ)^* (t)dt〗

小波逆变换为：

f(t)=∬▒〖Υ(s,τ) φ_(s,τ) (t)dsdτ〗

公式中*表示共轭，φ_(s,τ)表示小波基函数。小波基函数是由一个基本小波φ(t)放缩、平移形成。

φ_(s,τ) (t)=1/√s φ((t-τ)/s)

本研究使用 MATLAB 的 wavelet 工具箱里的连续小波变换函数，MATLAB 代码如下：

Xcwt(i,:)=CWT(X(i,:),scale,basis);

代码中 X(i,:)为第i个样本原始光谱矩阵，scale是小波分解尺度，取值范围1到正无穷，basis是小波基函数，Xcwt(i,:)是第i个样本光谱小波变换后的小波系数。

5. 提取图谱主成分

主成分分析的实质是对原坐标系进行平移和转换，使得新坐标原点与样本点集合的重心重合，新坐标的第一主轴与数据变异的最大方向对应，新坐标的第二主轴与第一主轴标准正交，并且对应于数据变异的第二大方向。其基本步骤如下：

（1）对数据进行标准化处理，得到矩阵X。也可以不处理，但是标准化处理可以使数据计算方便。

（2）计算矩阵X的协方差矩阵V。这时V又是X的相关系数矩阵。

（3）求V的前n个特征值Υ_1>Υ_2>…>Υ_n，以及特征值对应的特征向量α_1>α_2>…>α_m要求它们相互正交。

（4）求第h主成分F_h，有：

F_h=Xα_h=∑_(j=1)^p▒〖α_hj x_j 〗

式中：α_hj是主轴α_h的第j个分量。所以主成分F_h是原变量x_1,x_2, ,…, x_p的线性组合，组合系数恰好为α_hj。从这个角度，又可以说F_h是一个新的综合变量。

MATLAB代码如下：

[~,XSelected_SCORE,~]=princomp(X);

f1=XSelected_SCORE(:,1);

f2=XSelected_SCORE(:,2);

.

.

fn=XSelected_SCORE(:,n);

Xprincompn=[f1,f2,…,fn];

X 为输入光谱矩阵，XSelected_SCORE 是提取的全部主成分得分，f1,f2,…,fn 是第 1，2，…，n 个主成分。Xprincompn 是提取前 n 个主成分。

6. 相关系数法度量相似度

相关系数（通常用 R 表示）表明变量之间线性相关程度。在这里计算两个光谱谱图之间的相关系数来评价相似度，相关系数小于等于 1，越接近 1 越相似。公式如下：

R=(∑_(i=1)^n▒(x_ji-x ̅_j )(y_ki-y ̅_k ) )/√(∑_(i=1)^n▒(x_ji-x ̅_j )^2 ∑_(i=1)^n▒(y_ki-y ̅_k )^2 )

R为相关系数，x_ji，y_ki是j，k两个谱图第i点的透射率，x ̅_j和y ̅_k是j，k两个谱图在各个波数的均值，n根据谱图选取。MATLAB 实现代码如下：

r= corrcoef(X1, X2);

R=r(1,2);

代码中 X1,X2 代表两个光谱谱图数据，R 为要求的相关系数。corrcoef 是 MATLAB 计算相关系数矩阵的内置函数。