exinfomax.m

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

function [y,PI]=ExInfomax(x)

%
%扩展信息最大化
%编程时间：2006.1.12
%作者：张小兵
%*********初始化程序******************
%超高斯亚高斯
%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];
%加载两语音
%load wtwo
%s=s(:,5000:25000-1);
%s1=s(1,:);
%s2=s(2,:);
%混合
%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;
B=200;
[N,P]=size(x);
%********迭代程序**********************
for i=1:200
sum=0;

noblocks=fix(P/B);
for t=1:B:noblocks*B,  
    y=w*x(:,t:t+B-1);
    y1=y(1,:);
    k1=peak(y1);
    J1=sign(k1);
    y2=y(2,:);
    k2=peak(y2);
    J2=sign(k2);
    J=[J1,J2];
    J=diag(J);
    D=y*y'+J*tanh(y)*y';
    %af=diag(diag(D));
    af=eye(2);
    dw=afa*(af-D)*w;
    w=w+dw;
end;

%%%%%%%%%%%性能指标%%%%%%%%%%%%%%%%%%%%%%
G=w*wz*A;
n=length(P);
c=0;
for j=1:2
    a=0;
    b=0;
    for k=1:2
    a=a+abs(G(j,k))/max(abs(G(j,:)));
    b=b+abs(G(k,j))/max(abs(G(:,j)));
    end
    a=a-1;
    b=b-1;
   c=c+a+b;
end
PI(i)=c/(n*(n-1));
end
toc
y=w*x;
y1=y(1,:);
y2=y(2,:);
%*********快速ICA算法*******************
z=fastica(x);
z1=z(1,:);
z2=z(2,:);

%********画图程序**********************
figure(1);
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');