FASTICA盲源信号分离代码Matlab

FASTICA盲源信号分离代码Matlab

clc;clear all;close all;
% 读入混合前的原始图片并显示
[S1Fs1]=wavread （‘music.wav‘）; 
S1 = S1‘;
[S2Fs2]=wavread （‘man.wav‘）;
S2 = S2‘;
subplot（321）plot（S1）title（‘输入声音1‘）
subplot（322）plot（S2）title（‘输入声音2‘）
% 将其组成矩阵
S=[S1;S2];     
% weight=rand（size（S1））;               % 取一随机矩阵，作为信号混合的矩阵
% MixedS=weight*S;                      % 模拟麦克风获取信号
% 读取麦克风获取的信号
M1 = wavread （‘MixedS1.wav‘）‘;
M2 = wavread （‘MixedS2.wav‘）‘;
MixedS = [M1;M2];
% 将混合声音显示
subplot（323）plot（MixedS（1:））title（‘混合声音1‘）
subplot（324）plot（MixedS（2:））title（‘混合声音2‘）
% wavwrite（MixedS（1:）Fs1‘MixedS1.wav‘）;
% wavwrite（MixedS（2:）Fs1‘MixedS2.wav‘）;
MixedS_bak=MixedS;                         % 将混合后的数据备份，以便在恢复时直接调用
%%%%%%%%%%%%%%%%%%%%%%%%%%  中心化  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_mean=zeros（21）;
MixedS_mean=mean（MixedS2）;
MixedS = MixedS-repmat（MixedS_mean1size（MixedS2））;
%%%%%%%%%%%%%%%%%%%%%%%%%%%  白化  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MixedS_cov=cov（MixedS‘）;                    % cov为求协方差的函数
[ED]=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进行操作
[VariableNumSampleNum]=size（X）;
numofIC=VariableNum;                       % 在此应用中，独立元个数等于变量个数
B=zeros（numofICVariableNum）;              % 初始化列向量w的寄存矩阵B=[b1  b2  ...   bd]
for r=1:numofIC
    i=1;maxIterationsNum=100;               % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数）
    b=rand（numofIC1）-.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次迭代内并不收敛。‘ rmaxIterationsNum）;
      break;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%  ICA计算的数据复原并构图  %%%%%%%%%%%%%%%%%%%%%%%%%
ICAedS=B‘*Q*MixedS_bak;                     % 计算ICA后的矩阵

% 将混合矩阵重新排列并输出
subplot（325）plot（ICAedS（1:））title（‘ICA解混声音1‘）
subplot（326）plot（ICAedS（2:））title（‘ICA解混声音2‘）
wavwrite（ICAedS（1:）Fs1‘SS1.wav‘）;
wavwrite（ICAedS（2:）Fs1‘SS2.wav‘）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3027  2011-05-29 00:28  BBS\ICADemo.m

     文件     882044  2006-05-28 15:31  BBS\man.wav

     文件     882044  2011-05-29 00:20  BBS\MixedS1.wav

     文件     882044  2011-05-29 00:20  BBS\MixedS2.wav

     文件     882142  2006-05-28 15:32  BBS\music.wav

     文件     882044  2011-05-29 00:32  BBS\SS1.wav

     文件     882044  2011-05-29 00:32  BBS\SS2.wav

     目录          0  2011-05-29 00:32  BBS

----------- ---------  ---------- -----  ----

              5295389                    8