ewasobi.m

ICALAB for Signal Processing Toolbox for BSS, ICA, cICA, ICA-R, SCA and MCA

    end   
    gam1=sum(f2o.*x1(1:k,:),1);
    gam3=sum(f2o.*x3(1:k,:),1)+a0+phi(k,:);
    x4(k+1,:)=ones(1,M);
    b1m=sum(([phi(2:k+1,:); a0]+[zeros(1,M); phi(1:k,:)]).*x4(1:k+1,:));
    b2m=sum(([a0; phi(k+1:-1:2,:)]+phi(k+2:2*k+2,:)).*x4(1:k+1,:));
    latemp=b2m./alm;
    b2m=latemp-lalold; lalold=latemp;
    bom=alm./almold;
    ok=ones(k+1,1);
    x2(1:k+1,:)=x4(1:k+1,:).*(ok*(1./alm));
    x1(1:k+1,:)=[x1(1:k,:); zeros(1,M)]-(ok*gam1).*x2(1:k+1,:);
    x3(1:k+1,:)=[x3(1:k,:); zeros(1,M)]-(ok*gam3).*x2(1:k+1,:);
    x4temp=x4(1:k,:);
    x4(1:k+1,:)=[zeros(1,M); x4(1:k,:)]+[x4(2:k,:); ones(1,M); zeros(1,M)]...
       -(ok*bom).*[x4old; ones(1,M); zeros(1,M)]...
       -(ok*b2m).*x4(1:k+1,:)-(ok*b1m).*x1(1:k+1,:)-(ok*x4(1,:)).*x3(1:k+1,:);
    x4old=x4temp;
    almold=alm;
end
MK=M*K;
G=zeros(K,MK);
G(:,1:K:MK)=x1; clast=zeros(K,M);
f1=[phi(2:K,:); zeros(1,M)]+[zeros(1,M); phi(1:K-1,:)];
f2=[zeros(1,M); phi(K:-1:2,:)]+[phi(K+1:2*K-1,:); zeros(1,M)];
for k=2:K
    ck=G(:,k-1:K:MK);
    G(:,k:K:MK)=[ck(2:K,:); zeros(1,M)]+[zeros(1,M);  ck(1:K-1,:)]...
          -clast-(ok*sum(f1.*ck)).*x1-(ok*sum(f2.*ck)).*x2-(ok*ck(1,:)).*x3...
          -(ok*ck(K,:)).*x4;
    clast=ck;
end 

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   of THinv5

function [W] = ffdiag2(Rx_est,ini)
%
% computes Ziehe's ffdiag with pre-processing for maintaining
% high stability even for badly conditioned mixing matrices
% Next, ffdiag requires matrices stored in 3-dimensional field.
    [d Md]= size(Rx_est);
    M=Md/d;
    [H1,D] = eig(Rx_est(:,1:d));
    W_est = [];
    for k = 1:d
        W_est = [W_est; D(k,k)^(-.5)*H1(:,k)'/(H1(:,k)'*H1(:,k))];
    end
    R_zest=zeros(d,d,M);
    for k = 1:M
        id=(k-1)*d;
        R_zest(:,:,k) = W_est*Rx_est(:,id+1:id+d)*W_est';
    end
    %%%
        [W C stat] = ffdiag(R_zest,ini);  %%% CALLING ORIGINAL FFDIAG
    %%%    
    W=W*W_est;
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   of ffdiag2

function [V,C,stat] = ffdiag(C0,V0)
%FFDIAG  Diagonalizes a set of matrices C_k with a (single) non-orthogonal transformation. 
%
%
% Usage:
%
%  function [V,CD,stat] = ffdiag(C0,V0)  
%
% Input:
%    C0 - N x N x K array containing K symmetric matrices C_k which are to be diagonalized.
%    V0 - initial value for the diagonalizer (default eye(N))
%
% Output:
%    V    - joint diagonalization transformation (dimension N x N)
%    CD   - diagonalized set of matrices (dimension N x N x K).
%    stat - structure with some statistics of the algorithm's performance
%        .etime      - elapsed time
%        .errdiag    - diagonalization error
%  
% code by Andreas Ziehe and Pavel Laskov
%
%  (c) 2004  Fraunhofer FIRST.IDA 
%
% Usage example:  
%  
% gen_mat; [V,C,stat]=ffdiag(C0); imagesc(V*A);


  
%  
% The algorithm is derived in the paper:   
%  A. Ziehe, P. Laskov, G. Nolte and K.-R. Mueller, 
% "A Fast Algorithm for Joint Diagonalization with Non-orthogonal
%  Transformations and its Application to Blind Source Separation"
%  Journal of Machine Learning Research vol 5, pages 777-800, 2004.
%
% An earlier version has been presented at the ICA 2003 Workshop in
% Nara, Japan  
% 
%  A. Ziehe and P. Laskov and K.-R. Mueller and G. Nolte
%  "A Linear Least-Squares Algorithm for Joint Diagonalization",
%  in Proc. of the 4th International Symposium on 
%  Independent Component Analysis and Blind Signal Separation
%  (ICA2003), pages 469--474, Nara, Japan 2003.
% 
%  
% See also
% http://www.first.fraunhofer.de/~ziehe/research/ffdiag.html

  
  eps = 1e-9;

%theta = 0.1;    % threshold for stepsize heuristics 

[m,n,k] = size(C0);
C = C0;
Id=eye(n);
V = Id;

%more clever initialization ???
%[V,D]=eig(C0(:,:,1),C0(:,:,2));
%or
%[V,D]=eig(sum(C0,3));
%V=V';

inum = 1;
df = 1;

stat.method='FFDIAG';
tic;
while (df > eps & inum < 100)
  
  % Compute W
  W = getW(C);

  %old normalization  no longer recommended 
  %as it deteriorates  convergence
  % 
  %if (norm(W,'fro') > theta)
  %  W = theta/norm(W,'fro')*W;
  %end

  % A much better way is to 
  % scale W by power of 2 so that its norm is <1 .
  % necessary to make approximation assumptions hold
  % i.e. W should be small
  % 
  
  %if 0,
  [f,e] = log2(norm(W,'inf'));
  % s = max(0,e/2);

  s = max(0,e-1);
  W = W/(2^s );    
  
  % Compute update
  V = (Id+W)*V;

  %re-normalization
  V=diag(1./sqrt(diag(V*V')))*V;  %norm(V)=1 
  
  for i=1:k
    C(:,:,i) = V*C0(:,:,i)*V';
  end
   
  % Save stats
  stat.W(:,:,inum) = W;
  [f] = get_off(V,C0);

  stat.f(inum) = f;
  stat.errdiag(inum)=cost_off(C0,normit(V')'); 
  stat.nw(inum) = norm(W(:));
  
  
  if inum > 2
    df = abs(stat.f(inum-1)-stat.f(inum));

  end
  
%  fprintf(1,'itr %d :  %6.8f  \n ',inum, f);
  inum = inum+1;
end

stat.etime=toc;   %elapsed time
stat.niter=inum;  %number of iterations
stat.V=V;         %estimated diagonalizer

%subplot(2,1,1);
%semilogy(stat.f);
%title('objective function');

%subplot(2,1,2);
%semilogy(stat.nw');
%title('norm W');

end   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of ffdiag

function W = getW(C)
%getW   computes the update for FFDIAG.
% 
% Usage:
%    function W = getW(C)
%
% Input: 
%    C -  3-d array of matrices C_k
%
%  Output:
%    W - update matrix with zeros on the main diagonal
%
% See also:
%          FFDIAG
%

  
[m,n,K] = size(C);
W = zeros(n);
z = zeros(n);
y = zeros(n);

for i=1:n
  for j=1:n
    
    for k=1:K
      z(i,j) = z(i,j)+C(i,i,k)*C(j,j,k);
      y(i,j) = y(i,j)+0.5*C(j,j,k)*(C(i,j,k)+conj(C(j,i,k)));
    end
  end
end

for i=1:n-1
  for j=i+1:n
    W(i,j) = (z(j,i)*y(j,i)-z(i,i)*y(i,j))/(z(j,j)*z(i,i)-z(i,j)^2);
    W(j,i) = ((z(i,j)*y(i,j)-z(j,j)*y(j,i))/(z(j,j)*z(i,i)-z(i,j)^2));
  end
end
    
end         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of getW

function V=normit(V)
%NORMIT  Normalize the Frobenius norm of all columns of a matrix to 1. 
[n,m]=size(V);

for k=1:n,
  nn=norm(V(:,k),'fro');
  if nn<eps,warning('division may cause numerical errors');end
  V(:,k)=V(:,k)/nn;
end
end         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of normit

function [cost]=cost_off(C,V)
%COST_OFF computes diagonalization error as the norm of the off-diagonal elements.
%  C set of matrices in a 3-d array
%  V diagonalizing matrix


[N,N,K]=size(C);

if nargin>1,
  for k=1:K,
    C(:,:,k)=V*C(:,:,k)*V';
  end
end

cost=0;
for k=1:K
  Ck=C(:,:,k);Ck=Ck(:);
  Ck(1:N+1:N*N)=0;
  cost=cost+norm(Ck)^2;
end

end   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of cost_off

function [f,g] = get_off(V,C)
%GET_OFF  a wrapper for OFF 
%
% This function loops over 'off' and  adds up things over multiple matrices C.

[m,n,K] = size(C);
f = 0;
g = sparse(n,n);
h = sparse(n^2,n^2);
h1 = sparse(n^2,n^2);
h2 = sparse(n^2,n^2);

for k=1:K
  if nargout <= 1
    f = f + off(V,C(:,:,k));
  elseif nargout == 2
    [ff,gg] = off(V,C(:,:,k));
    f = f+ff;
    g = g+gg;
    gg
  end
end
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of get_off

function [f,g] = off(V,C)
%OFF - computes the value and the gradient of the "off" diagonality measure 
% 
% Usage:
%        [f,g] = off(V,C)
%
% code by Pavel Laskov and Andreas Ziehe
% 
% (c) 2004 Fraunfofer FIRST.IDA

F = V*C*V';
f = trace(F'*F) - trace(F.*F);

if nargout > 1
  g = 4*(F-diag(diag(F)))*V*C;
end
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of off  %%% and of the whole ffdiag