%block-extended-Infomax+fastICA
function [A,W,Y]=newICA(whitendsig,whiteningmatrix,dewhiteningmatrix,Numb,...
    elta,g,a1,a2,numIC,epsilon,numiteration,maxnumiteration);
%A，估計(jì)的混合矩陣
%W，解混矩陣
%Y,分解出的獨(dú)立分量
%whitendsig,白化后的矩陣
%whiteningmatrix,白化矩陣
%dewhiteningmatrix,解白化矩陣
%Numb,分塊的大小
%elta,學(xué)習(xí)率
%g,使用的非線性函數(shù)：pow3,tanh,gauss
%a1,tanh的參數(shù)
%a2,gauss的參數(shù)
%numIC,估計(jì)的獨(dú)立分量數(shù)，缺省值為觀測(cè)矩陣的維數(shù)
%epsilon,收斂標(biāo)準(zhǔn)，缺省值0.0001
%numiteration，Infomax算法的循環(huán)數(shù)，缺省值6
%maxnumiteration，fastICA算法的最大循環(huán)數(shù)，缺省值1000

if nargin<5,error('參數(shù)不夠！');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

%先進(jìn)行Infomax，得到一個(gè)W初始值
W1=eye(Dim);%推薦使用單位矩陣
W1old=zeros(size(W1));
I=eye(Dim);
K=eye(Dim);
numofblock=fix(numofsample/Numb) %把整個(gè)樣本分為numofblock個(gè)塊進(jìn)行處理
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

%再進(jìn)行fastICA

%檢查所使用的非線性函數(shù)
if strcmp(lower(g),'pow3'),
    usedNonlinearity=1;
elseif strcmp(lower(g),'tanh'),
    usedNonlinearity=2;
elseif strcmp(lower(g),'gauss'),
    usedNonlinearity=3;
else 
    error('錯(cuò)誤的非線性函數(shù)“g"= %s',g);
end
fprintf('使用非線性函數(shù) %s',g);

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

%固定點(diǎn)算法循環(huán)
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);     
    %計(jì)算混合矩陣
    A = dewhiteningmatrix * B;
    W = B'*whiteningmatrix;
    Y= B'*whitendsig;
    break;
    end
    Bold = B;
    %計(jì)算獨(dú)立分量U1=W*whitendsig,U的列向量為獨(dú)立分量
    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);
    %計(jì)算解混矩陣    
end