PCA.zip_FAULT DETECTION _pca fault detection_pca故障检测_proper7aj_故

close all;clear all; clc; X1 = load('.\data\66138891TE-data\d00.dat'); X2 = load('.\data\66138891TE-data\d14_te.dat'); X_train = X1'; X_test = X2; [n,d] = size(X_train); E_train=mean(X_train); %列均值 STD=1./std(X_train); %列标准差 X_train=X_train-repmat(E_train,n,1); %扩展为m行 X_train=X_train*diag(STD); %方差标准化 [n1,d1]=size(X_test); X_test=X_test-repmat(E_train,n1,1); %测试数据做同样处理 X_test=X_test*diag(STD); Rx = (1/(n-1))*(X_train')*(X_train); Rt = (1/(n1-1))*(X_test)'*(X_test); [U,S,V] = svd(Rx); [U1,S1,V1] = svd(Rt); len = length(S); allcon = sum(sum(S)); con = 0; for i = 1:len con = con + sum(S(i,:)); rate = con/allcon; if rate >=0.85 break end end Ek = i; P = U(1:Ek,:); %测试 % for i = 1:d1 % avgx = mean(X_test(:,i)); % for j = 1:n1 % X_test(j,i) = X_test(j,i) - avgx; % end % end % % for i = 1:d1 % stdx = std(X_test(:,i)); % X_test(:,i) = X_test(:,i)/stdx; % end % % % Tex = (1/(n1-1))*(X_test')*(X_test); % [U1,S1,V1] = svd(Tex); EX_test = P'*P*X_test'; E = X_test' - EX_test; len2 = length(E); VV = sum(S); VV = VV(1:Ek); VV = diag(VV); Tt = X_test *V(:,1:Ek); for i = 1:len2 SPE(i) = E(:,i)'*E(:,i); end cta1 = 0;cta2 = 0;cta3 = 0; for j=Ek+1:len cta1 = cta1+sum(S(j,:)); cta2 = cta2+(sum(S(j,:)))^2; cta3 = cta3+(sum(S(j,:)))^3; end h0 = 1 - (2*cta1*cta3)/(3*(cta2)^2); limit2 = 5; %取置信度为95%； Qalf = cta1*power((((limit2*h0*sqrt(2*cta2))/cta1)+1+((cta2*h0*(h0-1))/(cta1^2))),(1/h0)); figure; plot(1:length(SPE),SPE,'b'); hold on; plot(Qalf*ones(1,len2),'r'); xlabel('Samples'); ylabel('SPE'); SN = inv(VV); for i=1:n1 T_2(i) = Tt(i,:)*SN*Tt(i,:)'; end Falf = 1.5 ; %取a=0.05 T_2a = Falf*(Ek*(n1^2-1))/(n1*(n1-Ek)); T_2ai = T_2a*ones(n1,1); figure; plot(1:n1,T_2,'b'); hold on; plot(1:n1,T_2ai,'r'); xlabel('Samples'); ylabel('T^2'); Right = 0; for i = 1:161 if SPE(i)<=Qalf Right = Right +1; end end for i = 162:len2 if SPE(i)>=Qalf Right = Right +1; end end rate = Right/len2