MATLAB实现信号的盲源分离~fastICA方法~源信号与混合信号个数均为2~

clc;clear all;close all; % 读入混合前的原始图片并显示 [S1,Fs1]=wavread ('music.wav'); S1 = S1'; [S2,Fs2]=wavread ('man.wav'); S2 = S2'; subplot(3,2,1),plot(S1),title('输入声音1'), subplot(3,2,2),plot(S2),title('输入声音2'), % 将其组成矩阵 S=[S1;S2]; %weight=rand(size(S,1)); % 取一随机矩阵，作为信号混合的矩阵 % MixedS=weight*S; % 模拟麦克风获取信号 % 读取麦克风获取的信号 M1 = wavread ('MixedS1.wav')'; M2 = wavread ('MixedS2.wav')'; MixedS = [M1;M2]; % 将混合声音显示 subplot(3,2,3),plot(MixedS(1,:)),title('混合声音1'), subplot(3,2,4),plot(MixedS(2,:)),title('混合声音2'), % wavwrite(MixedS(1,:),Fs1,'MixedS1.wav'); % wavwrite(MixedS(2,:),Fs1,'MixedS2.wav'); MixedS_bak=MixedS; % 将混合后的数据备份，以便在恢复时直接调用 %%%%%%%%%%%%%%%%%%%%%%%%%% 中心化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MixedS_mean=zeros(2,1); MixedS_mean=mean(MixedS,2); MixedS = MixedS-repmat(MixedS_mean,1,size(MixedS,2)); %%%%%%%%%%%%%%%%%%%%%%%%%%% 白化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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; % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数） b=rand(numofIC,1)-.5; % 随机设置b初值 b=b/norm(b); % 对b标准化 norm(b):向量元素平方和开根号 while i<=maxIterationsNum+1 bOld=b; t=X'*b; g=t.*exp(-t.^2/2); dg=(1-t.^2).*exp(-t.^2/2); b=X*g/SampleNum-mean(dg)*b; b=b-B*B'*b; % 对b正交化 b=b/norm(b); if abs(abs(b'*bOld)-1)<1e-9 % 如果收敛，则 B(:,r)=b; % 保存所得向量b break; end i=i+1; end end if i == maxIterationsNum+1 % 循环结束处理 fprintf('\n第%d分量在%d次迭代内并不收敛。', r,maxIterationsNum); break; end %%%%%%%%%%%%%%%%%%%%%%%%%% ICA计算的数据复原并构图 %%%%%%%%%%%%%%%%%%%%%%%%% ICAedS=B'*Q*MixedS_bak; % 计算ICA后的矩阵 % 将混合矩阵重新排列并输出 subplot(3,2,5),plot(ICAedS(1,:)),title('ICA解混声音1'), subplot(3,2,6),plot(ICAedS(2,:)),title('ICA解混声音2'), wavwrite(ICAedS(1,:),Fs1,'SS1.wav'); wavwrite(ICAedS(2,:),Fs1,'SS2.wav');