% A simple demo program for EASI :
clear;
clc;
n=3; % number of sources
m=4; % number of sensors (not smaller than nsou !)
for t=1:5000 %给出源信号
%s(1,t)=sin(2*pi*90*t/10000);
s(1,t)=sin(2*pi*300*t/10000-6*cos(2*pi*60*t/5000));
s(3,t)=sin(2*pi*9*t/10000)*sin(2*pi*300*t/10000);
s(2,t)=sign(cos(2*pi*155*t/10000));
end
%s(5,:)=1-2*rand(1,T; %加入的噪声
[n,T]=size(s);
figure(1) %源信号图
for i=1:n
subplot(n,1,i);
plot(s(i,1:T));
end
A=1-2*rand(m,n); %混合矩阵
x=A*s; %根据混合矩阵A生成观测信号mixedsig
[m,T]=size(x);
figure(2) %观察信号图
for i=1:m
subplot(m,1,i);
plot(x(i,:));
end
%for i=1:n %观测信号零均值处理
% ave=mean(x(i,:));
% x(i,:)=x(i,:)-ave;
%end
u=0.005; % adaptation step
I=eye(m) ;
%B=0.2*eye(m,m); %初始化分离矩阵
B=0.5*eye(m);
for t=1:5000
% next line generates a vector of 'nsou' independent normalized QAM4 signals
y(:,t) =B*x(:,t);
%----------------- here is the algorithm ----------
%g = y(:,t) .* sqrt(diag(y(:,t)*y(:,t)'));
g=diag(y(:,t)*y(:,t)').*y(:,t) ; %%实现g=y(:,t).^3
% Other nonlinearities to be tried............
% g = y .* diag(0.1*y2).*diag(y2) ;
% g = y .* sqrt(diag(y2)) ;
% g = y .* log(diag(y2)) ;
% g = - y .* (diag(y2)<0.9) ;
% g = y ./ log(diag(y2)) ;
% g = - y ./ sqrt(diag(y2)) ;
% g = - y ./ diag(y2) ;
%----
%gy = g * y' ;
G=(I-y(:,t)*y(:,t)')/(1+u*(y(:,t))'*y(:,t)) +(y(:,t)*g' - g*(y(:,t))')/(1+u*abs((y(:,t))'*g)) ;
B=B+u*G*B ; % updating the separating matrix
%----------------- here ends the algorithm ----------
c=B*A; %系统矩阵
for p=1:m
max1(p)=abs(c(p,1));
for q=1:n
if max1(p)<=abs(c(p,q))
max1(p)=abs(c(p,q));
else max1(p)=max1(p);
end
end
end
s2=0;
for p=1:m
s1=0;
for q=1:n
s1=s1+abs(c(p,q))/max1(p);
end
s2=s2+abs(s1-1);
end
for q=1:n
max2(q)=abs(c(1,q));
for p=1:m
if max2(q)<=abs(c(p,q))
max2(q)=abs(c(p,q));
else max2(q)=max2(q);
end
end
end
s4=0;
for q=1:n
s3=0;
for p=1:m
s3=s3+abs(c(p,q))/max2(q);
end
s4=s4+abs(s3-1);
end
e(t)=s2/m+s4/n-(m-n)/n; %根据串音误差公式计算串音误差
end;
y=B*x;
figure(3) %估计源信号图,w0为前面估计的分离矩阵
for i=1:m %画出源信号的估计Y(t)的波形图
subplot(m,1,i);
plot(y(i,:));
end
figure(4) %画出串音误差图
plot(e)
评论3