%extended-Infomax+fastICA
function [A,W,Y]=newICA1(whitendsig,whiteningmatrix,dewhiteningmatrix,...
    elta,g,a1,a2,numIC,epsilon,numiteration,maxnumiteration);
%A，估計的混合矩陣
%W，解混矩陣
%Y,分解出的獨立分量
%whitendsig,白化后的矩陣
%whiteningmatrix,白化矩陣
%dewhiteningmatrix,解白化矩陣
%elta,學習率
%g,使用的非線性函數：pow3,tanh,gauss
%a1,tanh的參數
%a2,gauss的參數
%numIC,估計的獨立分量數，缺省值為觀測矩陣的維數
%epsilon,收斂標準，缺省值0.0001
%numiteration，Infomax算法的循環數，缺省值6
%maxnumiteration，fastICA算法的最大循環數，缺省值1000

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

%先進行Infomax，得到一個W初始值
W1=eye(Dim);%推薦使用單位矩陣
W1old=zeros(size(W1));
I=eye(Dim);
K=eye(Dim);
global U;
for round1 = 1:numiteration ;
    W1old = W1;
    U = W1*whitendsig;
    for i=1:Dim
        ui=U(i,:);
        k=sign(mean((sech(ui)).^2)*mean(ui.^2)-mean(ui.*tanh(ui)));
        K(i,i) = k;
    end
    dW1 = elta*(I-K*tanh(U)*U'-U*U')*W1;
    W1=W1+dW1; 
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