function KFDA()
clear
clc
width=120; %图像宽高
height=91;
nClass=15; %类别总数15类
per_ImgN=11;%每类共11个图像
trainN=8; %每类8个训练图像
allTrainImg=[];%所有训练图像
nType=[];%记录各训练图像所属类别
k=1;
for i=1:nClass
for j=1:trainN
imgdata=imread(strcat('E:\学习\资料\课件\模式识别\Fisher\Yale\',num2str(i),'\',num2str(j),'.BMP'));%读入图片
row=imgdata(1:width*height); % row是行矢量 1×N,其中N=120*91,提取顺序是先列后行,即从上到下,从左到右
column=row'; %column为列向量N*1
column=double(column);
allTrainImg=[allTrainImg,column]; %一列数据表示一个训练图像数据
nType(k)=i;
k=k+1;
end
end
%计算KPCA的核矩阵
ploya=2;ployd=2; %多项式核函数
kpcaMatrix = ((allTrainImg' * allTrainImg) + ploya).^ployd;
% 中心化核矩阵
sk = size(kpcaMatrix,1); % kpcaMatrix是方阵
row = sum(kpcaMatrix)/sk; % 列均值(行向量)
all = sum(kpcaMatrix(:))/(sk*sk); % 总均值
kpcaMatrix = kpcaMatrix - repmat(row,[sk 1]) - repmat(row',[1 sk]) + repmat(all,[sk sk]);
%求特征向量与特征值
[eigV eigD]=eig(kpcaMatrix);
d=diag(eigD); %提取矩阵的对角元素(即特征值组成的列向量)
[d1 index]=sort(d,'descend'); %降序排序,d1为排序后的向量,index为d1中每一项对应于d中项的索引
cols=length(find(d));%寻找非零特征值的个数
for i=1:cols
Vecsort(:,i) = eigV(:,index(i)); % vsort 是一个M*cols阶矩阵,保存的是按降序排列的特征向量,每一列构成一个特征向量
end %完成降序排列
d1=d1(1:cols);
Valsort=diag(d1);
%归一化特征向量
Vecsort = Vecsort*inv(sqrtm(Valsort));
%训练图像投影
trainF = Vecsort'*kpcaMatrix;
%在训练样本的KPCA特征矩阵上求得散布矩阵Sb、Sw、St
trainMean=mean(trainF,2);%求所有训练图像的均值(列向量)
classLabel = unique(nType);
priorProb=1/nClass; %先验概率
classMean=[];%求每类的均值
for i=1:nClass
index = find(nType==classLabel(i));
classMean(:,i)=mean(trainF(:,index),2);
end
Sb=0;
%求类间散布矩阵Sb
for i=1:nClass
Sb=Sb+(trainMean-classMean(:,i))*(trainMean-classMean(:,i))';
end
Sb=priorProb.*Sb;
Sw=0;%求类内散布矩阵Sw
for i=1:nClass
for j=trainN*(i-1)+1:trainN*i
g=trainF(:,j)-classMean(:,i);
Sw=Sw+g*g';
end
end
Sw=priorProb.*Sw;
St=Sb+Sw;
[V,S]=eig(Sw);
s=diag(S); %提取矩阵的对角元素(即特征值组成的列向量)
[s,index]=sort(s,'descend');%降序排序
Rank_Sw=rank(Sw);
%取Sw的零空间
k=1;
for i=(Rank_Sw+1):size(Sw,1)
P1(:,k)=V(:,index(i));
k=k+1;
end
Sb1=P1'*Sb*P1;
Rank_SbI=rank(Sb1);
[V1,S1]=eig(Sb1);
s1=diag(S1);
[s1,index1]=sort(s1,'descend');
for i=1:Rank_SbI
Z1(:,i)=V1(:,index1(i));
end
W1=P1*Z1;
%取P1的正交补空间
for i=1:Rank_Sw
P2(:,i)=V(:,index(i));
end
Sb2=P2'*Sb*P2;
St2=P2'*St*P2;
[V2,S2]=eig(Sb2,St2);
s2=diag(S2);
[s2,index2]=sort(s2,'descend');
k=1;
for j=1:size(s2,1)
if s2(j)>1e-6
Z2(:,k)=V2(:,index2(j));
k=k+1;
end
end
W2=P2*Z2;
W=[W1,W2];
ProjectSet=W;
trainF = ProjectSet'*trainF;%训练图像投影
allTestImg=[];%所有测试图像
for i=1:nClass
for j=trainN+1:per_ImgN
imgdata=imread(strcat('E:\学习\资料\课件\模式识别\Fisher\Yale\',num2str(i),'\',num2str(j),'.BMP'));%读入图片
row=imgdata(1:width*height); % row是行矢量 1×N,其中N=120*91,提取顺序是先列后行,即从上到下,从左到右
column=row'; %column为列向量N*1
column=double(column);
allTestImg=[allTestImg,column]; %一列数据表示一个训练图像数据
end
end
dim=size(allTrainImg,2);%allTrainImg的列数
kpcaTestMatrix = ((allTrainImg' * allTestImg ) + ploya).^ployd;
testF =ProjectSet'* Vecsort'*(kpcaTestMatrix - repmat(sum(kpcaTestMatrix,1),[dim 1])/dim);%测试图像投影
right=0;%做测试,记录正确识别数
for i=1:nClass
for j=1:(per_ImgN-trainN)
for k=1:trainN*nClass
juli(k)=norm(testF(:,(i-1)*(per_ImgN-trainN)+j)-trainF(:,k));%计算图片与训练样本之间的距离,如果距离最小,认为测试图片与训练图片属于同一类
end
[r t]=min(juli); %测试图片与训练图片t为一类,即属于[(t-1)/trainN]+1类
if fix((t-1)/trainN)==(i-1)
right=right+1;
end
end
end
right=right/(nClass*(per_ImgN-trainN))