浅谈盲源分离算法在故障诊断中的应用（附matlab代码）

clcclearclose all%% 仿真信号load A2（A2为文字内容中的混合矩阵A）N=1024*2; % 信号长度fs=1024;   % 采样频率f0=50;     %  基频t=0:1/fs:(N-1)/fs;s1=sin(2*pi*f0*t)+sin(2*pi*2*f0*t);s2=sin(2*pi*f0*t)+sin(2*pi*0.33*f0*t);s3=sin(2*pi*0.45*f0*t);S=[s1;s2;s3]; %源信号%仿真源信号时域、频域图像figure(1)subplot(321) plot(S(1,:))subplot(323) plot(S(2,:))subplot(325) plot(S(3,:))subplot(322)hua_fft(s1,fs,1,0,500)subplot(324)hua_fft(s2,fs,1,0,500)subplot(326)hua_fft(s3,fs,1,0,500)%% 混合信号及图像X=A*S+0.5*randn(2,N);%混合信号时域、频域图像figure(3)subplot(221) plot(X(1,:))subplot(223) plot(X(2,:))subplot(222) hua_fft(X(1,:),fs,1,0,500)subplot(224) hua_fft(X(2,:),fs,1,0,500)%% 盲源分离分解信号X_fft=X*dftmtx(N);%% 源信号信号在傅里叶基下的恢复(最短路径法)eS_fft1=SourceRecovery_ShortestPathMethon(X_fft,A);eS_1=real(eS_fft1/dftmtx(N));%% 恢复源信号时域图像figure(5)subplot(321) plot(eS_1(1,:))subplot(323) plot(eS_1(2,:))subplot(325) plot(eS_1(3,:))subplot(322) hua_fft(eS_1(1,:),fs,1,0,500)subplot(324) hua_fft(eS_1(2,:),fs,1,0,500)subplot(326) hua_fft(eS_1(3,:),fs,1,0,500)S_corr1=corr(S',eS_1');%% 适当降噪处理eS_fft2=abs(eS_fft1);c=max(eS_fft2,[],2);for j=1:3  eS_fft1(j,find(eS_fft2(j,:)<(c(j,1))/3))=0;endeS_1=real(eS_fft1/dftmtx(N));%% 降噪后源信号时域图像figure(6)subplot(321) plot(eS_1(1,:))subplot(323) plot(eS_1(2,:))subplot(325) plot(eS_1(3,:))subplot(322) hua_fft(eS_1(1,:),fs,1,0,500)subplot(324) hua_fft(eS_1(2,:),fs,1,0,500)subplot(326) hua_fft(eS_1(3,:),fs,1,0,500)S_corr2=corr(S',eS_1')；%%%%用到的函数function eS = SourceRecovery_ShortestPathMethon(X,A)% 把矩阵A映射到C_N^(m)个子空间中来恢复源，使得重构误差最小、% X --  观测信号（时频点复数）% A --  混合矩阵% eS -- 估计出的源信号n = size(A,2);m = size(X,1);N = size(X,2);Combine = nchoosek(1:n, m); % 从1~n中，每次取m，所有组合为Combine混合矩阵A构成m*m型矩阵的几种可能组合eS = zeros(n, size(X,2));  residual = zeros(N, size(Combine,1)); % 误差系数for j = 1:size(Combine,1) residual(:,j)=sum(abs(pinv(A(:,Combine(j,:)))* X).^2,1);end[~,u] = min(residual,[],2);for i = 1:size(Combine,1)    idx = find(u==i); % 找出每种组合中的最小误差对应的    eS(Combine(i,:),idx) = pinv(A(:,Combine(i,:)))* X(:,idx); % X(:,idx)表示误差最小的时频点endfunction hua_fft(y,fs,style,varargin)% 作用：画信号的幅值谱和功率谱% 输入：% y   --  输入信号（一个向量，行或列向量）% fs  --  采样频率% style - 当style=1,画幅值谱；当style=2,画功率谱;当style=其他的数字，那么画幅值谱和功率谱%       - 当style=1时，还可以多输入2个可选参数，如：hua_fft(y,fs,1,0,500)表示查看[0,500]段的频率图像（第一个为设置的频率起点，第二个为终点）% 输出：%         直接绘制出对应的图像%%y=y-mean(y);                  % 去除信号中的直流分量，防止对频谱分析造成干扰nfft= 2^nextpow2(length(y));  % 自动计算最佳FFT步长nfft）y_ft=fft(y,nfft);% 计算功率谱（相当于y轴）y_p=y_ft.*conj(y_ft)/nfft;    % conj()函数是求y函数的共轭复数，实数的共轭复数是他本身。如：conj(1+2i)=1-2i        % 幅值谱和功率谱的横坐标（频率轴）y_f=fs*(0:nfft/2-1)/nfft;     % DFT变换后对应的频率的序列         if style==1         % 画幅值谱    if nargin==3        plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));        ylabel('幅值');xlabel('频率');title('信号幅值谱');else           % 画出指定频率区间的幅值谱        f1=varargin{1};        fn=varargin{2};        ni=round(f1 * nfft/fs+1);        na=round(fn * nfft/fs+1);        plot(y_f(ni:na),abs(y_ft(ni:na)*2/length(y)));    end   elseif style==2    % 画功率谱    plot(y_f,y_p(1:nfft/2));    ylabel('功率谱密度');xlabel('频率');title('信号功率谱');   else               % 画幅值谱和功率谱    subplot(211);    plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));    ylabel('幅值');xlabel('频率');title('信号幅值谱');    subplot(212);    plot(y_f,y_p(1:nfft/2));    ylabel('功率谱密度');xlabel('频率');title('信号功率谱');endend