FastICA盲源分离算法语音信号处理

盲源分离FastICA算法仿真

这个例程展示了采用FastICA算法进行两个语音信号的解混合。

%% FastICA算法
% 输入信号为两段语音
% 采用
close all,clear all;clc;
[S1,fs1] = audioread('ICA\sound1.wav'); % 读取原始语音信号
[S2,fs2] = audioread('ICA\sound2.wav');
subplot(3,2,1),plot(S1),title('输入信号1');
subplot(3,2,2),plot(S2),title('输入信号2');

% 将其组成矩阵
s1 = S1';
s2 = S2';
S=[s1;s2];  
Sweight = rand(size(S,1));  
   
MixedS=Sweight*S;     % 将混合矩阵重新排列并输出

subplot(3,2,3),plot(MixedS(1,:)),title('混合信号1');
subplot(3,2,4),plot(MixedS(2,:)),title('混合信号2');

MixedS_bak=MixedS;                  
%%%%%%%%%%%%%%%%%%%%%%%%%%  标准化  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_mean=zeros(2,1);
for i=1:2
    MixedS_mean(i)=mean(MixedS(i,:));
end                                        % 计算MixedS的均值

for i=1:2
    for j=1:size(MixedS,1)
        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=1000;               % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数）
    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.^2).*exp(-a2*t.^2/2);
        b=((1-u)*t'*g*b+u*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

for k=1:numofIC
    W(:,k)=B(:,k)/5^k;          % 得到解混矩阵W
end
%%%%%%%%%%%%%%%%%%%%%%%%%%  ICA计算的数据复原并构图  %%%%%%%%%%%%%%%%%%%%%%%%%
ICAedS=W'*Q*MixedS_bak;                     % 计算ICA后的矩阵

% 将混合矩阵重新排列并输出
subplot(3,2,5),plot(ICAedS(1,:)),title('ICA解混信号1');
subplot(3,2,6),plot(ICAedS(2,:)),title('ICA解混信号2');
%%%%%%%%%%%%%%%%%%%%%%%%% 源信号、混合信号以及解混合之后的信号的读写以及播放%%%%%%%%%%%%%%%%%%%
% 对源信号进行播放
sound(S1);%对源信号1进行播放
sound(S2);%对源信号2进行播放

% 对混合后的音频信号进行保存并播放
audiowrite('ICA\mixed_wav1.wav',MixedS(1,:),8000);%保存wav数据
audiowrite('ICA\mixed_wav2.wav',MixedS(2,:),8000);%保存wav数据
sound(MixedS(1,:));%对混合信号1其进行播放
sound(MixedS(2,:));%对混合信号2其进行播放

%对解混合后的信号进行保存并播放
audiowrite('ICA\fastica_wav1.wav',ICAedS(1,:),8000);%保存wav数据
audiowrite('ICA\fastica_wav2.wav',ICAedS(1,:),8000);%保存wav数据
sound(ICAedS(1,:));%对解混信号1其进行播放
sound(ICAedS(2,:));%对解混信号2其进行播放

其中的音频文件的素材和仿真程序产生的混合及解混合的音频都可以在以下链接下载。
链接：https://pan.baidu.com/s/18rWQgNKO7_79WvwmAx_KrA
提取码：emeq