bs-ng.m

盲信号自然梯度法盲分离程序

lose all;  clear;   clc;

%% sources
load d18d3;   load h526h7;
s11 = d18d3;  s22 = h526h7;
p = 1;
q = 1;
s1 = resample(s11,p,q);        s2 = resample(s22,p,q);
N = 4000;                     sources = [s1(1:N);s2(1:N)];
[row_s,line_s] = size(sources);

%% delete the mean 
ss = sources;
for i = 1:row_s
    s = ss(i,:);
    s = s-mean(s).*ones(1,line_s);
    s = s./std(s);
    ss(i,:) = s;
end
sources = ss;

%% mixing matrix
d_lamda = 0.5;
theta1 = 4;          theta2 = 8;
tao1 = 2*pi*d_lamda*sin(theta1); tao2 = 2*pi*d_lamda*sin(theta2);
A = [1 1;exp(-j*tao1) exp(-j*tao2)];

%% mixtures
x1 = A*sources;
x = x1-mean(x1,2)*ones(1,N);

%% whiten signals
R = zeros(row_s,row_s);
for i = 1:row_s
    for j = 1:row_s
       R(i,j) =(x(i,:)*x(j,:)')/N;
    end
end      
[F,D]=eig(R);
for k=1:size(D,1)
    D(k,k)=1/sqrt(D(k,k));
end
V=D*F';                                       
v=V*x;

%% separate
I = eye(row_s);   P = I;   W = I;  lamda = 0.01;
for sep_t = 1:line_s
    y = W*v(:,sep_t);
    ys(:,sep_t) = y;
    if kurt(ys(1,:))<0; g = y-tanh(y); else g = y+tanh(y); end;
    W = W+lamda*(I-g*y')*W;
end

%% plot
figure(1);
for i = 1:row_s
    subplot(row_s,1,i);
    plot(sources(i,:));
end; 
subplot(row_s,1,1),title('Sources signal');
figure(2);
for i = 1:row_s
    subplot(row_s,1,i);
    plot(real(x(i,:)));
end; 
subplot(row_s,1,1),title('Mixtures signal');
figure(3);
for i = 1:row_s
    subplot(row_s,1,i);
    plot(real(ys(i,:)));
end; 
subplot(row_s,1,1),title('recovered signal');