附录3 本书自编MATLAB函数的源代码及注释

附录3 本书自编MATLAB函数的源代码及注释

本书介绍的一些算法和基本概念的应用需要在MATLAB下编辑命令得以实现。为便于读者使用和理解，将其汇编在本附录中。

1 精密度测试光谱的标准方差光谱的计算

本书第2章的（2-9）介绍了应用精密度测试光谱的标准方差评价光谱质量的公式与方法，其应用通过下述自编函数实现（该函数要求输入数据文件为csv格式）：

function [data,x,dirs,stdA]=SVSRS(filedir,filetype,k)

% 对多个光谱进行标准差分析（若各光谱的波长点不相等则该函数无法使用） % filedir—光谱文件所在目录，filetype—文件类型，两者必须以字符串形式输入。

% 例如SVSRS('c:/matlab/work/data/','*.csv'，0)即为对路径c:/matlab/work/data中的光谱进行标准方差计算

% k=0,对原始光谱进行标准方差分析；

% k=1，对光谱进行一阶导数处理后再进行标准方差分析； % k=2，对光谱进行SNV预处理后再进行标准方差分析； % k=3，对光谱进行MSC预处理后再进行标准方差分析 % x为光谱的横坐标（波数或波长），data为光谱吸光度矩阵（行为样本序号，列为光谱点序号）

% dirs为所计算路径下光谱文件名（字符串）组成的列向量 % stdA为指定路径下光谱的标准方差构成的行向量

d=dir([filedir filetype]); for i=1:length(d)

temp=load([filedir d(i).name]);

% 将filedir目录下第i个光谱文件中的数据赋给temp矩阵,temp的第1列为光谱横坐标（波数或波长点），第2列为该光谱的吸光度

data(i,:)=temp(:,2)';

% 将第i根光谱的吸光度赋给data矩阵，data的每一行为一根光谱

if i==1

dirs=d(i).name; else

dirs=strvcat(dirs,d(i).name); %将第i个光的文件赋给dirs end end

% 以上为读入光谱数据的命令

x=temp(:,1); % 将光谱横坐标赋给x

[n,m]=size(data); %将样本个数赋给n，光谱点个数赋给m if k<1 % 若k<1计算原始光谱的标准方差并作图 stdA=std(data); plot(x,stdA');

elseif k==1 % 若k=1对光谱进行一阶导数预处理后求其标准方差并作图 b1=gradient(data); %对光谱求一阶导数 for i=1:n

A(i,:)=smooth(b1(i,:)',11,'sgolay')'; %对导数光谱进行11点光滑 end

stdA=std(A); plot(x,stdA');

elseif k==2 %若k=2对光谱进行SNV预处理后求其标准方差并作图 A=(data-repmat(mean(data,2),1,m))./repmat(std(data,0,2),1,m);

% 对光谱数据进行SNV预处理并赋给A

stdA=std(A); plot(x,stdA')

else k==3 % 若k=3对光谱进行MSC预处理后求其标准方差并作图 % 以下对光谱数据矩阵data进行msc预处理 msp = mean(data); for i=1:n

B =regress(data(i,:)',[msp;ones(1,m)]',0.01); Amsc(i,:)=(data(i,:)-repmat(B(2,1),1,m))./B(1,1); end

% 完成对data的msc预处理 stdA=std(Amsc); plot(x,stdA'); end

% end of 标准方差光谱计算

例1：在MATLAB下分别输入如下命令：

[data,x,dirs,stdA]=SVSRS('E:\\基础化学计量学及其应用\\repeat/','*.csv',0); [data,x,dirs,stdA]=SVSRS('E:\\基础化学计量学及其应用\\repeat/','*.csv',1); [data,x,dirs,stdA]=SVSRS('E:\\基础化学计量学及其应用\\repeat/','*.csv',2); [data,x,dirs,stdA]=SVSRS('E:\\基础化学计量学及其应用\\repeat/','*.csv',3);

打开dirs可看到该路径下同一广藿香挥发油样品连续5次测得的5根近红外光谱文件名：

dirs=

ghx2001.CSV ghx2002.CSV ghx2003.CSV ghx2004.CSV ghx2005.CSV

分别得到下面4个标准方差光谱图：

广藿香挥发油精密度测试原始光谱的标准方差 广藿香挥发油精密度测试一阶导数光谱的标准方差

广藿香挥发油精密度测试SNV光谱的标准方差 广藿香挥发油精密度测试MSC光谱的标准方差

2 简化的KNN模式识别方法

8.4.1介绍了KNN模式识别方法及其简化的类重心模式识别方法并在例8-6与例8-7中演示了在MATLAB下处理该例的步骤与命令。基于主成分空间下马氏距离的类重心判别方法（即例8-7中的方法）完整和通用的M文件源代码如下：

function [cgroup_model,percent,vgroup_model,Acr,D]=SKNN(usage,feature,group,cum)

% 以离散变量为模式特征，以主成分空间下样本到各个类重心的马氏距离为判别依据的简化KNN模式识别方法

% usage 是样本的用途矩阵，1为训练集，2为预测集或检验集

% feature 为样本的特征变量矩阵,该矩阵行数为样本个数，列数为特征个数

% group为提供样本类别信息的列向量，对应的元素可以是数字也可以是字符串 % 例：对于3类问题，其类别分别用1，2，3表示，则group=[1 1 2 3]'表示第1，2个样本为第1类，第3个样本为第2类，第4个样本为第3类

% cum为设定的前m个主成分累积贡献率的阈值，一般大于0.9 % cgroup 是预测的检验集样本类别

% percent 是检验集预测正确率（即模型判别检验集正确的样本数与检验集样本个数之比）

% vgroup 检验集样本类别信息

% Acr 整体预测正确率（即模型判别建模集+检验集样本正确个数与总样本数之比）

% D存放每个样本到各类中心的马氏距离，行为样本序号，列为类别序号

if nargin < 3

error('本函数至少需要 3 个输入参数'); end

[n,d] = size(feature); %n为总样本数，d为特征个数 if size(group,1) ~= n

error('GROUP 的长度必需等于 feature 的长度'); end

[pc,score,latent]=princomp(feature);

% 对特征进行主成分分析，消除特征间的相关性 yt=cumsum(latent)./sum(latent);

%求每个特征值下对应的累计贡献率并赋给yt % 以下确定主成个数 M=[]; for i=1:d

if yt(i)>=cum m=i;

M=[M;m]; end

pcm=min(M); end pcm

lmd=latent(1:pcm)'

% 完成主成分个数的确定

smpidx=find(usage==2); %将检验集样本的序号赋给smpidx trnidx=find(usage==1); %将训练集样本的序号赋给trnidx if length(trnidx)<3

error ('训练集至少要含有3个样本'); end

score=score(:,1:pcm); %将选定的pcm个主成分得分赋给score

trnscore=score(trnidx,:); %将训练集的pcm个主成分得分赋给trnscore

tgroup=group(trnidx,:) %将训练集的类别信息赋给tgroup并屏幕显示结果 vgroup=group(smpidx,:) %将检验集的类别信息赋给vgroup并屏幕显示结果

[gindext,groupt] = grp2idx(tgroup);%将训练集样本的类别信息赋给列向量gindext

ngroupt = length(groupt); %将训练集中的总类别数赋给ngroupt gtsize = hist(gindext,1:ngroupt);

%统计训练集每个类中的样本数目并赋给行向量gsize，gsize中的每一列表示对应类的练训集样本数目

ntrn = length(trnidx); %将训练集样本数赋给ntrn

[gindexv,groupv] = grp2idx(vgroup); %将检验集样本的类别信息赋给列向量gindexv ngroupv = length(groupv); %将检验集中的总类别数赋给ngroupv gvsize = hist(gindexv,1:ngroupv);

nv = length(smpidx); %将检验集样本数赋给nv

% 以下在主成分空间计算训练集中每个类的类重心

gmeans = repmat(NaN, ngroupt, pcm); %开一个行数为类别数、列数为主成分个数的矩阵gmeans

for k = 1:ngroupt

gmeans(k,:) = mean(trnscore(find(gindext == k),:),1)

%将第k类训练集样本的主成分平均得分赋给gmeans的第k行 end

% 结束每个类的重心坐标（该类样本的主成分得分平均值）的计算并将结果存放在矩阵gmeans中

% 以下开始计算每个样本到各个类重心的马氏距离

D = repmat(NaN,n,ngroupt); %开一个行数为样本数、列数为类别数的矩阵D for k1 = 1:ngroupt

temp=gmeans(k1,:);

x1=(score - repmat(temp,n,1)).^2; x2=x1./repmat(lmd,n,1);

D(:,k1) = sqrt(sum(x2,2)); end

% 每个样本到各个类重心的马氏距离计算结束，存储在矩阵D中

[Dmin,G]=min(D,[],2); %将矩阵D的每一行的最小值赋给Dmin，D中该行最小值所在的列号赋给G

gaindex=repmat(NaN,n,1); %开一个行数为样本数n的列向量gaindex

gaindex(trnidx,:)=gindext; %将训练集类别信息按训练集的样本序号放入gaindex的对应行

gaindex(smpidx,:)=gindexv; %将检验集类别信息按检验集的样本序号放入gaindex的对应行

Acrn=sum(G==gaindex); %统计判别正确的样本总数（训练集+检验集） cgroup_model=G(trnidx);%将预测的训练集类别信息放在cgroup_model Acr=Acrn/n %计算整体正确率

vgroup_model=G(smpidx);%将预测的检验集类别信息放在vgroup_model crt=sum(vgroup==vgroup_model);%统计判断正确的检验集样本个数 if nv>0

percent=crt/nv %计算检验集预测正确率

本文来源：https://www.bwwdw.com/article/3962.html