对信号进行预白化处理，使信号的后处理更加容易，故障更加清晰matlab.zip

%%%第一段程序 clear all clc a=load('E:\date\a\06525193260001160313.TXT'); b=load('E:\date\a\06525193260001160314.TXT'); c=load('E:\date\a\06525193260001160315.TXT'); lev=3; xa=wden(a,'minimaxi','s','mln',lev,'db6'); xb=wden(b,'minimaxi','s','mln',lev,'db6'); xc=wden(c,'minimaxi','s','mln',lev,'db6'); K=3; fs=15360; t=1/fs; N=8192; f=[0:N/2-1]/(N*t); t1=[0:N-1]/fs; A=fft(xa); A1=abs(A)*2/N; B=fft(xb); B1=abs(B)*2/N; C=fft(xc); C1=abs(C)*2/N; figure(1) subplot(311); plot(t1,xa); xlabel('时间'); ylabel('振幅'); subplot(312); plot(t1,xb); xlabel('时间'); ylabel('振幅'); subplot(313); plot(t1,xc); xlabel('时间'); ylabel('振幅'); figure(2) subplot(311); plot(f,A1(1:N/2)); xlabel('频率'); ylabel('幅值'); subplot(312); plot(f,B1(1:N/2)); xlabel('频率'); ylabel('幅值'); subplot(313); plot(f,C1(1:N/2)); xlabel('频率'); ylabel('幅值'); %希尔伯特包络解调 figure(3) A2=abs(hilbert(xa)); A3=abs(fft(A2-mean(A2)))*2/N; subplot(211);plot(t1,A2);title('Hilbert包络波形'); subplot(212);plot([0:N-1]/(N*t),A3);title('Hilbert解调谱'); axis([0 1000 0 2]); figure(4) B2=abs(hilbert(xb)); B3=abs(fft(B2-mean(B2)))*2/N; subplot(211);plot(t1,B2);title('Hilbert包络波形'); subplot(212);plot([0:N-1]/(N*t),B3);title('Hilbert解调谱'); axis([0 1000 0 2]); figure(5) C2=abs(hilbert(xc)); C3=abs(fft(C2-mean(C2)))*2/N; subplot(211);plot(t1,C2);title('Hilbert包络波形'); subplot(212);plot([0:N-1]/(N*t),C3);title('Hilbert解调谱'); axis([0 1000 0 2]); A=randn(K); H=A*[A2;B2;C2]; for i=1:K H(i,:)=H(i,:)-1/N*sum(H(i,:));%分别对K个信号去均值，每一行是一个信号 end %figure(3); %plot(k,H); %title('观测信号') %实现对观察数据矩阵H的预白化 A=double(H); m=mean(A,2); A=A-m(:,ones(1,size(A,2)));%%上面三句把H转化为双精度型数据后在去均值，是预白化 covarianceMatrix=cov(A');%求协方差矩阵 covarianceMatrix=A*A'; [y,x]=eig(covarianceMatrix);%求特征值 x为由特征值构成的对角矩阵 y为特征向量构成的矩阵 whiteningMatrix=x^(-0.5)*y'; dewhiteningMatrix=y*sqrt(x); new_A=whiteningMatrix*A;%上面三句求白化矩阵 %figure(4) %plot(k,new_A); %title('whitened')%白化后的观测信号 %去相关 for i=1:3 for j=(i+1):3 temp=x(i,i); temp_2=y(:,i); if temp<x(j,j) x(i,i)=x(j,j); x(j,j)=temp; y(:,i)=y(:,j); y(:,j)=temp_2; end end end d=x; v=y; d.^(1/2); Y=inv(d.^(1/2))*v'*H; %inv求逆 %用FastICA算法实现对源信进行分离 epsilon=0.0001;%学习步长 W=rand(K); for p=1:K W(:,p)=W(:,p)/norm(W(:,p));%norm求最大的奇异值 exit=0; count=0; iter=1; while exit==0; count=count+1; temp=W(:,p); %记录上次迭代的值 W(:,p)=1/N*Y*((temp'*Y).^3)'-3*temp; ssum=zeros(K,1); for counter=1:p-1 ssum=ssum+(W(:,p)'*W(:,counter))*W(:,counter); end W(:,p)=W(:,p)-ssum; %正交化 W(:,p)=W(:,p)/norm(W(:,p)); if(abs((dot(W(:,p),temp)))<1+epsilon)&(abs((dot(W(:,p),temp)))>1-epsilon) %判断是否收敛 exit=1; end iter=iter+1; end end out=W'*Y; o1=out(1,:); o2=out(2,:); o3=out(3,:); figure(6); subplot(311); plot(t1,o1); subplot(312); plot(t1,o2); subplot(313); plot(t1,o3); No1=abs(fft(o1-mean(o1)))*2/N; No2=abs(fft(o2-mean(o2)))*2/N; No3=abs(fft(o3-mean(o3)))*2/N; figure(7) subplot(311);plot([0:N-1]/(N*t),No1); axis([0 1000 0 0.5]); subplot(312);plot([0:N-1]/(N*t),No2); axis([0 1000 0 0.5]); subplot(313);plot([0:N-1]/(N*t),No3); axis([0 1000 0 0.5]); %%第二段程序 clear all clc a=load('E:\date\a\06525193260001160313.TXT'); b=load('E:\date\a\06525193260001160314.TXT'); c=load('E:\date\a\06525193260001160315.TXT'); lev=3; xa=wden(a,'minimaxi','s','one',lev,'db6'); xb=wden(b,'minimaxi','s','one',lev,'db6'); xc=wden(c,'minimaxi','s','one',lev,'db6'); K=3; fs=15360; t=1/fs; N=8192; f=[0:N/2-1]/(N*t); t1=[0:N-1]/fs; A=fft(xa); A1=abs(A)*2/N; B=fft(xb); B1=abs(B)*2/N; C=fft(xc); C1=abs(C)*2/N; figure(4) subplot(311); plot(t1,a); xlabel('时间'); ylabel('振幅'); subplot(312); plot(t1,b); xlabel('时间'); ylabel('振幅'); subplot(313); plot(t1,c); xlabel('时间'); ylabel('振幅'); figure(1) subplot(311); plot(t1,xa); xlabel('时间'); ylabel('振幅'); subplot(312); plot(t1,xb); xlabel('时间'); ylabel('振幅'); subplot(313); plot(t1,xc); xlabel('时间'); ylabel('振幅'); figure(2) subplot(311); plot(f,A1(1:N/2)); xlabel('频率'); ylabel('幅值'); subplot(312); plot(f,B1(1:N/2)); xlabel('频率'); ylabel('幅值'); subplot(313); plot(f,C1(1:N/2)); xlabel('频率'); ylabel('幅值'); %希尔伯特包络解调 %figure(3) A2=abs(hilbert(xa)); A3=abs(fft(A2-mean(A2)))*2/N; %subplot(211);plot(t1,A2);title('Hilbert包络波形'); %subplot(212);plot([0:N-1]/(N*t),A3);title('Hilbert解调谱'); %axis([0 1000 0 2]); %figure(4) B2=abs(hilbert(xb)); B3=abs(fft(B2-mean(B2)))*2/N; %subplot(211);plot(t1,B2);title('Hilbert包络波形'); %subplot(212);plot([0:N-1]/(N*t),B3);title('Hilbert解调谱'); %axis([0 1000 0 2]); %figure(5) C2=abs(hilbert(xc)); C3=abs(fft(C2-mean(C2)))*2/N; %subplot(211);plot(t1,C2);title('Hilbert包络波形'); %subplot(212);plot([0:N-1]/(N*t),C3);title('Hilbert解调谱'); %axis([0 1000 0 2]); figure(3); subplot(311); plot([0:N-1]/(N*t),A3);title('Hilbert解调谱'); axis([0 1000 0 3]); subplot(312); plot([0:N-1]/(N*t),B3);title('Hilbert解调谱'); axis([0 1000 0 3]); subplot(313); plot([0:N-1]/(N*t),C3);title('Hilbert解调谱'); axis([0 1000 0 3]); % 将其组成矩阵 MixedS=[A2;B2;C2]; % 信号个数即为变量数 % 因此S_all是一个变量个数＊采样个数的矩阵 %Sweight=rand(size(S,1)); % 取一随机矩阵，作为信号混合的权矩阵 %MixedS=Sweight*S; % 得到三个信号的混合信号矩阵 % 将混合矩阵重新排列并输出 %subplot(4,3,4),plot(MixedS(1,:)),title('混合信号1'),axis([0,100,-4,4]); %subplot(4,3,5),plot(MixedS(2,:)),title('混合信号2'),axis([0,100,-4,4]); %subplot(4,3,6),plot(MixedS(3,:)),title('混合信号3'),axis([0,100,-4,4]); MixedS_bak=MixedS; % 将混合后的数据备份，以便在恢复时直接调用 %%%%%%%%%%%%%%%%%%%%%%%%%% 标准化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MixedS_mean=zeros(3,1); for i=1:3 MixedS_mean(i)=mean(MixedS(i,:)); end % 计算MixedS的均值 for i=1:3 for j=1:size(MixedS,2) MixedS(i,j)=MixedS(i,j)-MixedS_mean(i); end end %%%%%%%%%%%%%%%%%%%%%%%%%%% 白化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MixedS_cov=cov(MixedS'); % cov为求协方差的函数 [E,D]=eig(MixedS_cov); % 对信号矩阵的协方差函数进行特征值分解 Q=inv(sqrt(D))*(E)'; % Q为白化矩阵 MixedS_white=Q*MixedS; % MixedS_white为白化后的信号矩阵 IsI=cov(MixedS_white'); % IsI应为单位阵 %%%%%%%%%%%%%%%%%%%%%%%%　FASTICA算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% X=MixedS_white; % 以下算法将对X进行操作 [VariableNum,SampleNum]=size(X); numofIC=VariableNum; % 在此应用中，独立元个数等于变量个数 B=zeros(numofIC,VariableNum); % 初始化列向量w的寄存矩阵,B=[b1 b2 ... bd] for r=1:numofIC i=1;maxIterationsNum=100; % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数） IterationsNum=0; b=rand(numofIC,1)-.5; % 随机设置b初值 b=b./norm(b); % 对b标准化 norm(b):向量元素平方和开根号 while i<=maxIterationsNum+1 if i == maxIterationsNum % 循环结束处理 fprintf('\n第%d分量在%d次迭代内并不收敛。', r,maxIterationsNum); break; end bOld=b; a2=1; u=1; t=X'*b; g=t.*exp(-a2*t.^2/2); dg=(1-a2*t.