ewasobi.m

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

function [Wwasobi,Wsobi,ISR,signals]= ewasobi(x,AR_order,rmax,eps0); 
%
% implements algorithm WASOBI for blind source separation of
% AR sources in a fast way, allowing separation up to 100 sources
% in the running time of the order of tens of seconds.
%
% INPUT:    x .... input data matrix d x N
%           d .... dimension of the data
%           N .... length of the data
%           ARmax .. maximum AR order of the separated sources
%           rmax ... a constant that may help to stabilize the algorithm.
%                    it has the meaning of maximum magnitude of poles of the
%                    AR sources. The choice rmax=1 means that no stabilization 
%                    is applied. The choice rmax=0.99 may lead to more stable
%                    results.
%           eps0 ... machine dependent constant to control condition number
%                    of weight matrices
%
% OUTPUT: Wwasobi ...... estimated de-mixing matrix
%         Wsobi ........ initial estimate of the matrix obtained by FFDiag2
%         ISR .......... estimated ISR matrix which represents approximate accuracy 
%                        of the separation provided that there is no additive 
%                        noise in the model.
%         signals....... separated signals
%
% Code by Petr Tichavsky, using inputs from Eran Doron, Pavel Laskov 
% and Andreas Ziehe. A paper that describes the algorithm is submitted 
% to ICASSP 2007.
%
if nargin<4
   eps0=5.0e-7;
end
if nargin<3
   rmax=0.99;
end  
num_of_iterations = 3;
[d N]=size(x);
Xmean=mean(x,2);
x=x-Xmean*ones(1,N);  %%%%%%%%%  removing the sample mean
%[AOL_1 Rx_est] = sobi(x,AR_order); WOL_1=inv(AOL_1);
T=length(x(1,:))-AR_order;
C0=corr_est(x,T,AR_order);
%
Wsobi = ffdiag2(C0,eye(d)); 
Wwasobi=Wsobi; WOL_1=Wsobi;
Rx_est = [];
%
for in = 1:num_of_iterations
    [WOL_1,Rx_est,AR]=newasobi(x,AR_order+1,WOL_1,Rx_est,rmax,eps0);
    Wwasobi=WOL_1*Wwasobi;
end
ISR=CRLB4(AR)/N;
sqdnorm=diag(Wwasobi*C0(1:d,1:d)*Wwasobi');
Wwasobi=Wwasobi./(sqrt(sqdnorm)*ones(1,d));
signals=Wwasobi*x+(Wwasobi*Xmean)*ones(1,N);
end  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of EWASOBI
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [W_est_WASOBI,Rx_est,AR]=newasobi(x,M,W_pre,Rx_est1,rmax,eps0);
%
%   Computes one iteration needed for wasobi
%
%   d ... number of independent components
%   M ... number of lags = AR_order+1
%   Vnorm ... estimated demixing matrix obtained in preprocessing
%   Wnorm ...  estimated demixing matrix obtained by the advanced tuning
%

status=1;
[d,N]=size(x);
d2=d*(d-1)/2;
T=length(x(1,:))-M+1;
if isempty(Rx_est1)
    Rx_est=corr_est(x,T,M-1,W_pre);
else
   for index=1:M
       id=(index-1)*d;
       Rx_est(:,id+1:id+d)=W_pre*Rx_est1(:,id+1:id+d)*W_pre';
%       Rx_est(:,:,index)=W_pre*Rx_est1(:,:,index)*W_pre';
   end
end
R=zeros(M,d);
for index=1:M
    id=(index-1)*d;
    R(index,:)=diag(Rx_est(:,id+1:id+d)).'; 
%     R(index,:)=diag(Rx_est(:,:,index)).'; %%% columns of R will contain 
                           %%% covariance function of the separated components
end
%
[AR,sigmy]=armodel(R,rmax);      %%% compute AR models of estimated components
%
AR3=[]; 
for i=2:d
  for k=1:i-1
    AR3=[AR3 conv(AR(:,i),AR(:,k))];
  end  
end
phi=ar2r(AR3);     %%%%%%%%%% functions phi to evaluate CVinv
CVinv=THinv5(phi,M,d2,eps0*phi(1,:));  %%%% to compute inversions of CV 
                                       %%%% It has dimension zeros(M,M*d2).
im=1; 
for i=2:d
  for k=1:i-1
     fact=1/(sigmy(1,i)*sigmy(1,k));
     CVinv(:,(im-1)*M+1:im*M)=CVinv(:,(im-1)*M+1:im*M)*fact;
     im=im+1;
  end
end  
W_est_WASOBI=serial10(Rx_est, CVinv);  %%% TO MINIMIZE THE WEIGHTED CRITERION

end %%%%%%%%%%%%%%%%%%%%%%%%%%   of newasobi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ r ] = ar2r( a )
%%%%%
%%%%% Computes covariance function of AR processes from 
%%%%% the autoregressive coefficients using an inverse Schur algorithm 
%%%%% and an inverse Levinson algorithm (for one column it is equivalent to  
%%%%%      "rlevinson.m" in matlab)
% 
  if (size(a,1)==1)
      a=a'; % chci to jako sloupce
  end
  
  [p m] = size(a);    % pocet vektoru koef.AR modelu
  alfa = a;
  K=zeros(p,m);
  p = p-1;
  for n=p:-1:1
      K(n,:) = -a(n+1,:);
      for k=1:n-1
          alfa(k+1,:) = (a(k+1,:)+K(n,:).*a(n-k+1,:))./(1-K(n,:).^2);
      end
      a=alfa;
  end
%  
  r  = 1./prod(1-K.^2);
  f=[r(1,:); zeros(p,m)];
  b=f;
  for k=1:p 
      for n=k:-1:1
          f(n,:)=f(n+1,:)+K(n,:).*b(k-n+1,:);
          b(k-n+1,:)=-K(n,:).*f(n+1,:)+(1-K(n,:).^2).*b(k-n+1,:);
      end
      b(k+1,:)=f(1,:);
      r=[r; f(1,:)]; 
  end       
end %%%%%%%%%%%%%%%%%%%%%%%%%%%  of ar2r
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function R_est=corr_est(x,T,q,W_pre)
%
[K,KK]=size(x);
R_est=[];
if ~exist('W_pre')
   for index=1:q+1
       R_est=[R_est 1/T*(x(:,1:T)*x(:,index:T+index-1)')];
   end    
else
   for index=1:q+1
       R_est=[R_est 1/T*W_pre*(x(:,1:T)*x(:,index:T+index-1)')*W_pre'];
   end    
end

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of corr_est
%
%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


function [AR,sigmy]=armodel(R,rmax)
%
% to compute AR coefficients of the sources given covarience functions 
% but if the zeros have magnitude > rmax, the zeros are pushed back.
%
[M,d]=size(R); v1=zeros(1,d); v2=v1; Rs=R;
for id=1:d
    AR(:,id)=[1; -toeplitz(R(1:M-1,id),R(1:M-1,id)')\R(2:M,id)];
    v=roots(AR(:,id)); %%% mimicks the matlab function "polystab"
%    v1(1,id)=max(abs(v));
    vs=0.5*(sign(abs(v)-1)+1);
    v=(1-vs).*v+vs./conj(v);
    vmax=max(abs(v));
%    v2(1,id)=max(abs(v));
    if vmax>rmax
       v=v*rmax/vmax;
    end   
    AR(:,id)=real(poly(v)'); %%% reconstructs back the covariance function
end 
Rs=ar2r(AR);
sigmy=R(1,:)./Rs(1,:);
% [v1; v2]
end %%%%%%%%%%%%%%%%%%%%%%%  of armodel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function W=serial10(Rx_est,CVinv);
%
[d Md]=size(Rx_est);
M=Md/d; dd=d^2-d; dd2=dd/2;
W=eye(d);
%
num_of_iteration=5;
for iteration=1:num_of_iteration
    A0=eye(d);
    im=1; 
    for i=2:d
      for k=1:i-1
        Y(im,:)=(Rx_est(i,k:d:Md)+Rx_est(k,i:d:Md))/2;  
%        Y(im,:)=squeeze(Rx_est(i,k,:)+Rx_est(k,i,:))/2; 
        im=im+1;
      end
    end
    for i=1:d  
      LaMat(i,:)=Rx_est(i,i:d:Md);
    end 
    b11=zeros(dd2,1); b12=b11; b22=b11; c1=b11; c2=c1;
    m=0; 
    for id=2:d        
      for id2=1:id-1
          m=m+1; im=(m-1)*M;
          Wm=CVinv(:,im+1:im+M);
          Wlam1=Wm*LaMat(id,:)';
          Wlam2=Wm*LaMat(id2,:)';
          b11(m,1)=LaMat(id2,:)*Wlam2;
          b12(m,1)=LaMat(id,:)*Wlam2;
          b22(m,1)=LaMat(id,:)*Wlam1;
          c1(m,1)=c1(m,1)+Wlam2'*Y(m,:)';
          c2(m,1)=c2(m,1)+Wlam1'*Y(m,:)';
       end
    end
    det0=b11.*b22-b12.^2; 
    d1=(c1.*b22-b12.*c2)./det0; 
    d2=(b11.*c2-b12.*c1)./det0;
%    value=norm([d1; d2])
    m=0;
    for id=2:d
        for id2=1:id-1
            m=m+1;
            A0(id,id2)=d1(m,1);
            A0(id2,id)=d2(m,1);
        end
    end  
    Ainv=inv(A0);  
    if iteration<num_of_iteration       
       for m=1:M
 %          Rx_est(:,:,m)=Ainv*Rx_est(:,:,m)*Ainv';
           id=(m-1)*d;
           Rx_est(:,id+1:id+d)=Ainv*Rx_est(:,id+1:id+d)*Ainv';
       end
    end
    W=Ainv*W;
end     
end  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  of serial10


function ISR = CRLB4(AR);
%
% CRLB4(AR) generates the CRLB for gain matrix elements (in term 
% of ISR) for blind separation of K Gaussian autoregressive sources 
% whose AR coefficients (of the length M, where M-1 is the AR order)
% are stored as columns in matrix AR.

[M K]=size(AR);

Rs=ar2r(AR);

sum_Rs_s=zeros(K,K);

for s=0:M-1
    for t=0:M-1
        sum_Rs_s=sum_Rs_s+(AR(s+1,:).*AR(t+1,:))'*Rs(abs(s-t)+1,:);
    end
end

denom=sum_Rs_s'.*sum_Rs_s+eye(K)-1;
ISR=sum_Rs_s'./denom.*(ones(K,1)*Rs(1,:))./(Rs(1,:)'*ones(1,K));
ISR(eye(K)==1)=0;

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% of CRLB4

function G=THinv5(phi,K,M,eps);
%
%%%% Implements fast (complexity O(M*K^2))
%%%% computation of the following piece of code:
%
%C=[];
%for im=1:M 
%  A=toeplitz(phi(1:K,im),phi(1:K,im)')+hankel(phi(1:K,im),phi(K:2*K-1,im)')+eps(im)*eye(K);
%  C=[C inv(A)];
%end  
%
% DEFAULT PARAMETERS: M=2; phi=randn(2*K-1,M); eps=randn(1,2);
%   SIZE of phi SHOULD BE (2*K-1,M).
%   SIZE of eps SHOULD BE (1,M).

phi(2*K,1:M)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
almold=2*phi(1,:)+eps;
C0=1./almold;
x1=zeros(K,M); x2=x1; x3=x1; x4=x1;
x1(1,:)=C0; x2(1,:)=C0; 
x3(1,:)=-C0.*phi(2,:); 
x4(1,:)=-2*C0.*phi(2,:); 
x4old=[];
lalold=2*phi(2,:)./almold;
for k=1:K-1
    f2o=phi(k+1:-1:2,:)+phi(k+1:2*k,:);
    alm=sum(f2o.*x4(1:k,:),1)+phi(1,:)+eps+phi(2*k+1,:);
    a0=zeros(1,M); 
    if k<K-1
       a0=phi(k+2,:);