ngica.asv

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

function y=ngica(x)

%基于自然梯度的ICA
%预白化（零均值，白化）
%激活函数（非线性函数）
%f(y)=y.^3
%学习速率u(t)
%u(t)=0.001*exp(-0.1*t)
%dw=u(t)*(I-fy(:,t)*y(:,t)')*W;
%其中I=diag(diag(fy(:,t)*y(:,t)'))
%性能指标函数为F
%F=1/N*sum(1/(n-1)*sum(abs(p)/max(abs(p))-1)
%P=W*A
%
%2005.1.6
%author:xiaobing zhang

%预处理
[n,T] = size(x);
m=mean(x')'*ones(1,T);
x = x - m;% Remove mean
[F,D]=eig((x*x')/T);
v=F*(D^(-0.5))*F';
x=v*x;
%取初值
W=rand(3);
I=eye(3);
y=zeros(3,length(x));
fy=zeros(3,length(x));
%主要运算
for i=1:100
    for t=1:T
        eta=0.001*exp(-0.1*t);
        y(:,t)=W*x(:,t);
        fy(:,t)=y(:,t).^3;
        I=diag(diag(fy(:,t)*y(:,t)'));
        dw=eta*(I-fy(:,t)*y(:,t)')*W;
        W=W+dw;
        W=W/norm(W);   
    end
    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  
end
y=W*x;