求解这段指令的意思，越详细越好，谢谢啦

function [Kp,T2]=KPCA(ax,ay)

[Nx]=size(ax);

mean_X =mean(ax);

axb=ax;

std_X=std(ax);

ax=ax-mean_X(ones(Nx,1),:);

std_X(find(std_X==0))=1;%数据预处理

ax=ax./std_X(ones(Nx,1),:);

c=10000;

% gama=0.05;

% ni=1;

% F1=ax(1,:);

% F=F1';

% for ib=2:Nx

% for i=1:ni

% for j=1:ni

%

% batar1(ib).block(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);

% end

% batar2(i,1)=exp(-norm(ax(i,:)-ax(ib,:))^2/c);

% batar3(1,i)=exp(-norm(ax(ib,:)-ax(i,:))^2/c);

% end

% s1=exp(-norm(ax(ib,:)-ax(ib,:))^2/c);

% batar= batar3(1,i)*inv(batar1(ib).block)* batar2(i,1);

% s=(s1- batar)/s1;

% if s>sqrt(gama)

% ni=ni+1;

% F=[F ax(ib,:)'];

% end

%

%

% end

% AX=F'%训练集基的提取结束

[N]=size(ax,1);

for i=1:N

for j=1:N

K(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);%求核矩阵

end

end

n1=ones(N,N);

N1=1/N*n1;

Kp=K-N1*K-K*N1+N1*K*N1;

[u,s,v]=svd(Kp);

lamda=s*(1/(N-1));

P=v;

lamda=diag(lamda);

B=length(find(lamda>1e-10)); %求非零的特征值个数

%求主元个数

npc=1;

while sum(lamda(1:npc))/sum(lamda(1:B))<0.9

npc=npc+1;

end

%求特征空间有效维数

DimFS=1;

while sum(lamda(1:DimFS))/sum(lamda(1:B))<=0.99

DimFS=DimFS+1;

end

lamda=diag(lamda);

for i=1:B

% P(:,i)=P(:,i)/norm(P(:,i)*s(i,i));

% P(:,i)=P(:,i)/(norm(P(:,i))*sqrt(N*lamda(i,i)));

P(:,i)=P(:,i)/(sqrt((N-1)*lamda(i,i)));

end

[Ny]=size(ay,1);

mean_X =mean(axb);

std_X = std(axb);

[num_sample] = Ny;

ay = ay-mean_X(ones(num_sample,1),:);

ay = ay./std_X(ones(num_sample,1),:);

% mean_y = mean(ay);

% std_y=std(ay);

% ay = ay-mean_y(ones(Ny,1),:);

% std_y(find(std_y==0))=1;%数据处理

% ay = ay./std_y(ones(Ny,1),:);

for i=1:Ny

for j=1:N

Ky(i,j)=exp(-norm(ay(i,:)-ax(j,:))^2/c);

end

end

t1=ones(1,N);

t11=1/N*t1;

for i=1:Ny

kp1(i,:)= Ky(i,:)-t11*K- Ky(i,:)*N1+t11*K*N1;

end

for i=1:Ny

for k=1:B

t(i,k)=P(:,k)'*kp1(i,:)';

end

end

% 求T2,SPE

% covtyb=inv(t'*t);

for i=1:Ny

T2(i)=t(i,1:npc)*inv(lamda(1:npc,1:npc))*t(i,1:npc)'; %也可以

% SPE(i)=t(i,1:npc)*t(i,1:npc)';

% T2(1,i)=t(i,1:npc)*(covtyb(1:npc,1:npc))*t(i,1:npc)';

% SPE(i)=t(i,(npc+1):B)*t(i,(npc+1):B)';

SPE(i)=1-(2/N)*ones(N,1)'*Ky(i,:)'+(1/(N^2))*ones(N,1)'*K*ones(N,1)-t(i,1:npc)*t(i,1:npc)';

end

%T2,SPE控制线

t2cl=npc*(N-1)*(N+1)*icdf('f',0.99,npc,N-npc)/(N*(N-npc));

for i=1:3

theta(i)=trace((lamda(npc+1:DimFS,npc+1:DimFS))^i);

end

h0=1-2*theta(1)*theta(3)/(3*theta(2)^2);

ca=icdf('norm',0.99,0,1);

s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0)*4;

for i=1:Ny

t2cl1(i)=t2cl;

end

for i=1:Ny

s_cl1(i)=s_cl;

end

figure

plot(1:Ny,T2(1:Ny),'k');

hold on;

plot(1:Ny,t2cl1(1:Ny),'r');

title('T2');

hold off;

figure

plot(1:Ny,SPE(1:Ny),'k')

hold on;

plot(1:Ny, s_cl1(1:Ny),'r');

title('SPE');

hold off;