newica.asv

一种新的ICA算法

%block-extended-Infomax+fastICA
function [A,W,Y]=newICA(whitendsig,whiteningmatrix,dewhiteningmatrix,Numb,...
    elta,g,a1,a2,numIC,epsilon,numiteration,maxnumiteration);
%A，估计的混合矩阵
%W，解混矩阵
%Y,分解出的独立分量
%whitendsig,白化后的矩阵
%whiteningmatrix,白化矩阵
%dewhiteningmatrix,解白化矩阵
%Numb,分块的大小
%elta,学习率
%g,使用的非线性函数：pow3,tanh,gauss
%a1,tanh的参数
%a2,gauss的参数
%numIC,估计的独立分量数，缺省值为观测矩阵的维数
%epsilon,收敛标准，缺省值0.0001
%numiteration，Infomax算法的循环数，缺省值6
%maxnumiteration，fastICA算法的最大循环数，缺省值1000

if nargin<5,error('参数不够！');end
[Dim,numofsample] = size(whitendsig);
if nargin<12,maxnumiteration=10000;end
if nargin<11,numiteration=6;end
if nargin<10,epsilon=0.0001;end
if nargin<9,numIC=Dim;end
if nargin<8,a2=1;end
if nargin<7,a1=1;end
if nargin<6,g='pow3';end

%先进行Infomax，得到一个W初始值
W1=eye(Dim);%推荐使用单位矩阵
W1old=zeros(size(W1));
I=eye(Dim);
K=eye(Dim);
numofblock=fix(numofsample/Numb) %把整个样本分为numofblock个块进行处理
global U;

for N = 1:numofblock  
    for round1 = 1:numiteration ;
        W1old = W1;
        U(:,((N-1)*Numb+1):1:((N-1)*Numb+Numb)) = W1*whitendsig(:,((N-1)*Numb+1):1:((N-1)*Numb+Numb));
        for i=1:Dim
            ui=U(i,((N-1)*Numb+1):1:((N-1)*Numb+Numb));
            k=sign(mean((sech(ui)).^2)*mean(ui.^2)-mean(ui.*tanh(ui)));
            K(i,i) = k;
        end
        dW1 = elta*(I-K*tanh(U(:,((N-1)*Numb+1):1:((N-1)*Numb+Numb)))*U(:,((N-1)*Numb+1):1:((N-1)*Numb+Numb))'-U(:,((N-1)*Numb+1):1:((N-1)*Numb+Numb))*U(:,((N-1)*Numb+1):1:((N-1)*Numb+Numb))')*W1;
        W1=W1+dW1
    end
end

%再进行fastICA

%检查所使用的非线性函数
if strcmp(lower(g),'pow3'),
    usedNonlinearity=1;
elseif strcmp(lower(g),'tanh'),
    usedNonlinearity=2;
elseif strcmp(lower(g),'gauss'),
    usedNonlinearity=3;
else 
    error('错误的非线性函数“g"= %s',g);
end
fprintf('使用非线性函数 %s',g);

%将Infomax得到的W1作为W的初始值
A=zeros(Dim);
B=orth((W1*dewhiteningmatrix)');
Bold=zeros(size(B));

%固定点算法循环
for round = 1:maxnumiteration + 1;
    if round == maxnumiteration + 1;
        fprintf('在 %d 步后未收敛',maxnumiteration);
    break;
    end
    B = B*real(sqrtm(inv(B'* B)));
    %判断是否收敛
    e = min(abs(diag(B'*Bold)));
    if(1-e < epsilon)
        fprintf('在 %d 步后收敛',round);     
    %计算混合矩阵
    A = dewhiteningmatrix * B;
    W = B'*whiteningmatrix;
    Y= B'*whitendsig;
    break;
    end
    Bold = B;
    %计算独立分量U1=W*whitendsig,U的列向量为独立分量
    U1 = whitendsig'*B;
    if usedNonlinearity==1,
       B = (whitendsig * (U1.^3))/numofsample - 3 * B;
    elseif usedNonlinearity==2,
       hypTan = tanh(a1*U1);
       B = whitendsig * hypTan/numofsample - ones(size(B,1),1) * sum(ones(size(U1)).*(1-hypTan.^2)).* B/numofsample*a1;
    else usedNonlinearity==3,
        Usquared=U1.^2;
        gauss = U1.*exp(-a2*Usquared/2);
        dGauss = (1-a2*Usquared).*exp(-a2*Usquared/2);
        B = whitendsig*gauss/numofsample - ones(size(B,1),1) * sum(ones(size(U1)).*dGauss).* B/numofsample;
    end
    B=B/norm(B);
    %计算解混矩阵    
end