fenli1.m

自己写的一个fastica程序

clc;clear all;close all;

I1=wavread ('E:\fz\left.wav')';         
I2=wavread ('E:\fz\right.wav')';


subplot(2,2,1),plot(I1),title('混合声音1'),
subplot(2,2,2),plot(I2),title('混合声音2'),


MIXS=[I1;I2];                         
                                    
MIXS_bak=MIXS;                       
% 标准化  
MIXS_mean=zeros(2,1);
for i=1:2
    MIXS_mean(i)=mean(MIXS(i,:));
end                                        % 计算MIXS的均值

for i=1:2
    for j=1:size(MIXS,2)
        MIXS(i,j)=MIXS(i,j)-MIXS_mean(i);
    end
end

%  白化 

MIXS_cov=cov(MIXS');                    % cov为求协方差的函数
[E,D]=eig(MIXS_cov);                      % 对矩阵的协方差函数进行特征值分解
Q=inv(sqrt(D))*(E)';                        % Q为白化矩阵
MIXS_white=Q*MIXS;                      % MixedS_white为白化后的矩阵
IsI=cov(MIXS_white');                     % IsI应为单位阵            

%FASTICA算法  
X=MIXS_white;                           
[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标准化
    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
%    B(:,r)=b;                                % 保存所得向量b
end

ICAS=B'*Q*MIXS_bak;                     % 计算ICA后的矩阵

subplot(2,2,3),plot(ICAS(1,:)),title('ICA解混声音1'),
subplot(2,2,4),plot(ICAS(2,:)),title('ICA解混声音2'),
wavwrite(ICAS(1,:),'E:\fz\ica1.wav');
wavwrite(ICAS(2,:),'E:\fz\ica2.wav');