bksa.m

这是盲信号的代码 都已经通过编译了 做这方面的同仁可以参考一下 我觉得蛮惯用的

function [y,PI]=bksa()
%
%基于峰度的盲源开关算法
%系统仿真学报 2005.1
%copyrigth：2005.4.13
%author：lucky zhang
%*********init******************
fs=8000;
 f=14.324;
 load chirp;
 y=y(1:1000);
 k=1:1000;
s1=sin(2*pi*f*k/fs);
s2=y';
s=[s1;s2];
A=[1.5 1;1 1.4];
%A=randn(2);
x=A*s;
x=baihua(x);
x1=x(1,:);
x2=x(2,:);
w=eye(2);
afa=0.001;
%********迭代程序**********************
for i=1:50
y=w*x;
y1=y(1,:);
k1=peak(y1);
J1=sign(k1);
y2=y(2,:);
k2=peak(y2);
J2=sign(k2);
J=[J1,J2];
J=diag(J);
sum=0;
for k=1:1000
    D(:,:,k)=y(:,k)*y(:,k)'+J*tanh(y(:,k)*y(:,k)');
    af(:,:,k)=diag(diag(D(:,:,k)));
    dw=afa*(af(:,:,k)-D(:,:,k))*w;
    w=w+dw;
end
%%%%%%%%%%%性能指标%%%%%%%%%%%%%%%%%%%%%%
P=w*A;
n=length(P);
c=0;
for j=1:2
    a=0;
    b=0;
    for k=1:2
    a=a+abs(P(j,k))/max(abs(P(j,:)));
    b=b+abs(P(k,j))/max(abs(P(:,j)));
    end
    a=a-1;
    b=b-1;
   c=c+a+b;
end
PI(i)=c/(n*(n-1));
xn(i)=xnzb(P);
end
%*********快速ICA算法*******************


%********画图程序**********************
subplot(521);
plot(s1);
ylabel('s1');
subplot(522);
plot(s2);
ylabel('s2');
subplot(523);
plot(x1);
ylabel('x1');
subplot(524);
plot(x2);
ylabel('x2');
subplot(525);
plot(y1);
ylabel('y1');
text(400,1, '本算法')
subplot(526);
plot(y2);
ylabel('y2');