close;?
clc;?
clear?all;?
d00_ori=load?（‘D:\TE_data\d00_te.dat‘）;?%960*52
%d00_ori=load?D:\TE_data\d00_te.dat;?%960*52
d05_ori=load（‘D:\TE_data\d21_te.dat‘）;%960*52


dataNormal=d00_ori（:[123456789101112131415161718192021224243444546474849505152]）;
dataF=d05_ori（:[123456789101112131415161718192021224243444546474849505152]）;


%%%---------------------ICA----------------------------------開始
X=dataNormal;?
[Xmusigma]=zscore（X）;?%注意：正確的注解：cov（X‘）的對角線元素全為1，即diag（cov（X‘）?）=diag（?（X‘?*?X）/（n-1）?）=?I，
[icasigAW]?=?fastica?（X‘?‘approach‘?‘symm‘?‘g‘?‘tanh‘）;%mixedsig的行是信號種類，列是樣本數量，

ICs=9;

icasigDom=icasig（1:ICs:）;%%[9??960]取前面ICs個獨立成分，則得到的I2的數值達到8000類似于葛志強文章的結果
for?k=1:size（icasigDom2）
????I2Dom（k）=icasigDom（:k）‘*icasigDom（:k）;??%公式（21）
end
I2min=min（I2Dom）;?%進行kernel?density?estimation
I2max=max（I2Dom）;
xi?=?linspace（I2minI2max201）;?%201是插值點總數
fKde?=?ksdensity（I2Domxi‘function‘‘cdf‘）;
kdePoint=find（fKde>0.9887&fKde<0.9903）;

figure（1）;
subplot（311）;
plot（I2Dom（?kdePoint（1）?）?*?ones（size（dataNormal1）1）‘r--‘）
hold?on

dataF2=dataF;
XF2=dataF2;
for?j=1:size（dataF22）
????XF2_std（:j）=（XF2（:j）-mu（j））/sigma（j）;?%%歸一化方法來自文獻《combined?multivariate?statistical?process?control》，即利用正常工況數據的均值、方差
end
s_newd=W（1:ICs:）*XF2_std‘;%%取前面ICs個獨立成分，則得到的I2的數值達到8000類似于葛志強文章的結果

for?k=1:size（s_newd2）
????I2（k）=s_newd（:k）‘*s_newd（:k）;??%公式（21）
end
plot（I2）%%故障數據歸一化，得到的圖形很像文獻的Fig5的（c）
ylabel（‘I^2‘）;
hold?off
weak?=?find（（I2?-?I2Dom（?kdePoint（1）））<0）;
weak?=?find（weak>160）;
weak?=?size（weak2）;?%%%列的大小即為個數
fprintf（2‘I^2漏報比例百分之?%g?\n‘weak./8）;%%%%%%%%%%%%%%%%后800個對應于漏報
%%%---------------------ICA-------------------------------