ica_test.m

来自「基于独立分量分析进行盲信号分离的原理简介和例程」· M 代码 · 共 143 行

% [e-mail]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Generate simulate data         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% data length
N = 200000;

% read data from wave file
s1 = wavread('w1.wav'); s1 = s1(1:N)';

s2 = wavread('w2.wav'); s2 = s2(1:N)';
s3 = wavread('w3.wav'); s3 = s3(1:N)';
s = [s1; s2; s3];

% mix the data
A = [   1       4       5
        3.4     2.4     2
        6       2       3.5];
x = A * s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we'll recover s1, s2 & s3 from x %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

len = size(x, 2);   		% number of data samples
IC_num = size(x, 1);		% number of ICs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Center the data with 0 mean       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for (cnt = 1: IC_num)
    x(cnt, :) = x(cnt, :) - mean(x(cnt, :)')';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Whiten the data            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Rx = x*x'/len;
[Ux Sx] = eig(Rx);  		% assert that Ux*Sx*Ux'-Rx = 0

Vx = inv(Sx).^0.5*Ux';		% Vx is the whiten matrix
z = Vx * x;				% z is whiten data
                        % assert that z*z'/len = I

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Recover ICs                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% What kind of G we will be used ?
% G(y) = (1/alpha) * ln(cosh(alpha*y))
% g(y) = tanh(alpha*y)
%       or
% G(y) = -exp(-y^2/2)
% g(y) = y * exp(-y^2/2)
% g'(y) = (1-y^2)*exp(-y^2/2)

w = eye(IC_num);
w_old = w;

eps = 1.0e-5;
flag = 1;
while (flag == 1)
    for (cnt = 1: IC_num)
        y = w(:, cnt)'*z;
        Gy = -exp(-y.^2/2);
        gy = y.*(-Gy);
        gy1 = (1-y.^2).*(-Gy);
        E1 = z*gy'/len;
        E2 = mean(gy1')';
    
        w(:, cnt) = E1-E2*w(:, cnt);
    end

    wt = w';
    [Uw Sw] = svd(wt*wt');
    wt = inv(Uw * sqrt(Sw)) * wt;
    w = wt';

    change = ((sum(sum((abs(w'*w_old) - eye(IC_num)).^2)))/IC_num)
    if ( change < eps)
        flag = 0;
    end
    
    w_old = w;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 		    Recover s              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s_n = w'*z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 		Play all the music         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Play mix sound 1
x(1,:)=x(1,:)/(max(s_n(1,:))-min(x(1,:)))*2;
disp('press a key to hear mix sound 1');
pause;
wavplay(x(1, :), 22050);

% Play mix sound 2
x(2,:)=x(2,:)/(max(x(2,:))-min(x(2,:)))*2;
disp('press a key to hear mix sound 2');
pause;
wavplay(x(2, :), 22050);

% Play mix sound 3
x(3,:)=x(3,:)/(max(x(3,:))-min(x(3,:)))*2;
disp('press a key to hear mix sound 3');
pause;
wavplay(x(3, :), 22050);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Play sound 1
s_n(1,:)=s_n(1,:)-mean(s_n(1,:));
s_n(1,:)=s_n(1,:)/(max(s_n(1,:))-min(s_n(1,:)))*2;
disp('press a key to hear sound 1');
pause;
wavplay(s_n(1, :), 22050);

% Play sound 2
s_n(2,:)=s_n(2,:)-mean(s_n(2,:));
s_n(2,:)=s_n(2,:)/(max(s_n(2,:))-min(s_n(2,:)))*2;
disp('press a key to hear sound 2');
pause;
wavplay(s_n(2, :), 22050);

% Play sound 3
s_n(3,:)=s_n(3,:)-mean(s_n(3,:));
s_n(3,:)=s_n(3,:)/(max(s_n(3,:))-min(s_n(3,:)))*2;
disp('press a key to hear sound 3');
pause;
wavplay(s_n(3, :), 22050);