快速ICA程序代码

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-------------------------------