bksa.m

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

function PI=bksa()
%
%牛龙的基于峰度的盲源开关算法
%编程时间：2005.2.20
%作者：张小兵
%*********初始化程序******************
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];
x=A*s;
[x,wz]=baihua(x);
x1=x(1,:);
x2=x(2,:);
w=eye(2);
afa=0.001;
%********迭代程序**********************
for i=1:15
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)));
    %af(:,:,k)=eye(2);
    dw=afa*(af(:,:,k)-D(:,:,k))*w;
    w=w+dw;
end
%%%%%%%%%%%性能指标%%%%%%%%%%%%%%%%%%%%%%
P=w*wz*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));
end
%*********快速ICA算法*******************
z=fastica(x);
z1=z(1,:);
z2=z(2,:);

%********画图程序**********************
subplot(421);
plot(s1);
ylabel('s1');
subplot(423);
plot(s2);
ylabel('s2');
subplot(425);
plot(x1);
ylabel('x1');
subplot(427);
plot(x2);
ylabel('x2');
subplot(422);
plot(y1);
ylabel('y1');
text(400,1, '本算法')
subplot(424);
plot(y2);
ylabel('y2');
subplot(426);
plot(z1);
ylabel('z1');
text(400,1, 'fastica');
subplot(428)
plot(z2);
ylabel('z2');