www.gusucode.com > 基于PCA的故障信号检测源码程序 > 基于PCA的故障信号检测源码程序/pca/db5_PCA1.m

    
close;clc;clear all;

%--------------------------------------------------------------------------
% 载入数据

d1=xlsread('zc.xls');  
d2=xlsread('gz.xls');  
%--------------------------------------------------------------------------
% 数据的标准化

stantData=d1(:,:); %引入标准数据
c=d2;              %引入待检测数据
c1=ones(960,1);
Xmean=mean(stantData);%标准数据的均值
Xstd=std(stantData);  %标准数据的方差          
X=(c-c1*Xmean)./(c1*Xstd); %用待检测数据减去标准数据的均值并除以方差后，得到标准化的待检测数据  
 %X=zscore(c); %标准化                         %在此之所以用标准数据，是为了知道检测数据的偏离程度，因而减去标准数据的均值和方差
                           %为什么在此用标准数据的均值和方差
[P, T, LATENT, TSQUARED] = princomp(X);   %  LATENT是X的特征值；T=T得分向量；P=P负荷矩阵；（即T=XP）
%--------------------------------------------------------------------------
%确定主元素个数k

percent = input('确定方差贡献率限(0～1）：');       %  alpha0 是预设的贡献率,通常为85%
k=0;
for i=1:size(LATENT,1);                     %应用size函数返回LATENT的行维数（语法r=size(A,n)返回n指定的A的方向的维数：n=1为行方向的维数；n=2为列方向的维数）
    alpha(i)=sum(LATENT(1:i))/sum(LATENT); %求解贡献率
    if alpha(i)>=percent ;                  %判断是否大于给定贡献率 
        k=i;                               %把主元个数提供给K
        break;       
    end;
end;
%--------------------------------------------------------------------------
%存储新的变量

Xp=zeros(size(X));    

for j=1:k;   
    Xp=Xp+T(:,j)*P(:,j)';  %Xp是新的变量
end
%--------------------------------------------------------------------------

beta = input('确定统计阈值置信水平(0～1）：');%0.9 --0.99

theta=zeros(4,1);
for i=1:3 ;       
    for j=k+1:size(X,2) ;    %j=k+1:52  
        theta(i)=theta(i)+LATENT(j)^(i);    
    end
end
h0=1-2*theta(1)*theta(3)/[3*(theta(2)^2)];%h0=1-2*theta(1)*theta(3)/3/(theta(2)^2);
SPEbeta=theta(1)*(norminv(beta)*(2*theta(2)*h0^2)^0.5/theta(1)+1+theta(2)*h0*(h0-1)/(theta(1)^2))^(1/h0);%求取SPE的阈值
T2knbeta=finv(beta,k,(size(X,1)-k))*k*(size(X,1)^2-1)/(size(X,1)*(size(X,1)-k));   %利用F分布来求T2统计的阈值
%--------------------------------------------------------------------------
%求取T2的公式

T2=zeros(960,1);
 for i=1:960;
   T2(i)=X(i,:)*P(:,1:k)*inv(diag(LATENT(1:k)))*P(:,1:k)'*X(i,:)'; % 书上的公式
 end
%--------------------------------------------------------------------------
%  求取Q的公式

 Q=zeros(960,1);
for i=1:960;
  Q(i)=X(i,1:k)*(eye(k)-P(1:k,1:k)*P(1:k,1:k)')'*(eye(k)-P(1:k,1:k)*P(1:k,1:k)')*X(i,1:k)';   % 书上的公式
 end
%--------------------------------------------------------------------------

figure(1);
plot(T2(1:960));
hold on;
TS=T2knbeta*ones(960,1);  %产生全1矩阵{size(X,1) x 1}维，目的是画一条直线
plot(TS,'r--');
xlabel('采样次数');
ylabel('T2统计图');
hold off;
print -djpeg 144T2;  %保存t2统计图像
k 
SPEbeta
T2knbeta
%--------------------------------------------------------------------------

figure(2);
plot(Q(1:960));   
hold on;
SPES=SPEbeta*ones(960,1);%产生全1矩阵{size(X,1) x 1}维，目的是画一条直线
plot(SPES,'r--');
xlabel('采样次数');
ylabel('SPE统计图');
hold off;
print -djpeg 144SPE;  %保存SPE统计图像