波分析-Fast Ica(独立成分分析)算法

I1 = wavread('speechmixture1.wav');
I2 = wavread('speechmixture2.wav');
I3 = wavread('speechmixture3.wav');
I4 = wavread('speechmixture4.wav');
subplot(4,3,1),plot(I1),title('Input1');
subplot(4,3,2),plot(I2),title('Input2');
subplot(4,3,3),plot(I3),title('Input3');
subplot(4,3,3),plot(I4),title('Input4');
 
% 将其组成矩阵
II1 = I1';
II2 = I2';
II3 = I3';
II4 = I4';
S=[II1;II2;II3;II4];    
Sweight=rand(size(S,1));      
MixedS=Sweight*S;     % 将混合矩阵重新排列并输出
subplot(4,3,4),plot(MixedS(1,:)),title('Input1');
subplot(4,3,5),plot(MixedS(2,:)),title('Input2');
subplot(4,3,6),plot(MixedS(3,:)),title('Input3');
subplot(4,3,6),plot(MixedS(4,:)),title('Input4');
MixedS_bak=MixedS;                  

MixedS_mean=zeros(3,1);
for i=1:4
    MixedS_mean(i)=mean(MixedS(i,:));
end                                     
for i=1:4
    for j=1:size(MixedS,2)
        MixedS(i,j)=MixedS(i,j)-MixedS_mean(i);
    end
end

MixedS_cov=cov(MixedS');                 
[E,D]=eig(MixedS_cov);            
Q=inv(sqrt(D))*(E)';                
MixedS_white=Q*MixedS;                   
IsI=cov(MixedS_white');                 
 

X=MixedS_white;                           
[VariableNum,SampleNum]=size(X);
numofIC=VariableNum;                   
B=zeros(numofIC,VariableNum);          
for r=1:numofIC
    i=1;maxIterationsNum=100;        
    IterationsNum=0;
    b=rand(numofIC,1)-.5;               
    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/norm(b); 
        if abs(abs(b'*bOld)-1)<1e-9      
             B(:,r)=b;                   
             break;
         end
        i=i+1;        
    end
%    B(:,r)=b;  
end

ICAedS=B'*Q*MixedS_bak;

subplot(4,3,7),plot(ICAedS(1,:)),title('ICA解混信号1');
subplot(4,3,8),plot(ICAedS(2,:)),title('ICA解混信号2');
subplot(4,3,9),plot(ICAedS(3,:)),title('ICA解混信号3');
subplot(4,3,9),plot(ICAedS(4,:)),title('ICA解混信号4');
%保存wav
% wavwrite(ICAedS(1,:),8000,8,'ica_speechmixture1.wav');
% wavwrite(ICAedS(2,:),8000,8,'ica_speechmixture2.wav');
% wavwrite(ICAedS(3,:),8000,8,'ica_speechmixture3.wav');
% wavwrite(ICAedS(3,:),8000,8,'ica_speechmixture4.wav');