% C-FICA algorithm - time-domain Convolutive extension of FastICA  - v1.2
% 
% J. Thomas, Y. Deville and S. Hosseini
% "Time-domain fast fixed-point algorithms for convolutive ICA",
% IEEE Signal Processing Letters, vol. 13, no. 4, pp. 228-231, April 2006. 
% 
% 
% [Sc,sc,Isc,y,fe,fc,te,tc]=cfica(x,non-lin,epsi,Re,[we1 we2],Rc,[wc1 wc2],itermax,initializ,verbose)
% 
% Inputs (place the deconvolution parameters before the recoloration ones
% themselves before iter_max)
% 
% x         : observation vector
% non-lin   : 'k'->kurtosis, 'g'->Gaussian, 't'->tanh   default 'g'
% initializ : 'i' for initialization with unit filters  optional
% epsi      : stopping criterion (between 0 and 1)      default 10^-12
% Re        : order of the extraction filters           default 80
% [we1 we2] : extraction time window                    optional
% Rc        : order of the recoloration filters         default 3R1
% [wc1 wc2] : recoloration time window                  optional
% itermax   : maximum number of iteration               default 10000
% verbose   : 'v' for verbose                           optional
% 
% Outputs
%
% Sc        : most powerful contributions for each estimated source
% sc        : estimated source contributions  (tensor N*N*nb_samples)
% Isc       : indices of the most powerful contributions
% y         : N-1 estimated innovation processes  (matrix (N-1)*nb_samples)
% fe        : extraction filters
% fc        : recoloration filters
% te        : times of extraction
% tc        : times of recoloration
% -------------------------------------------------------------------------
% History of the software:
% V1.1: original software written by J. Thomas
%
% V1.2: Error messages appeared with recent releases of Matlab in the function
% conv_mix_nc. This version corrects them. 
% Author: Matthieu PUIGT   Contact: matthieu.puigt@univ-littoral.fr
%
% -------------------------------------------------------------------------



function [Sc,sc,Isc,y,fe,fc,te,tc]=cfica(x,varargin)

% N: number of observations, n: number of samples
[N,n]=size(x);

% centering of the observations
x=x-mean(x,2)*ones(1,n);

% ordering of the input arguments
[criterion,initializ,epsi,Re,we1,we2,Rc,wc1,wc2,iter_max,verbose]=args(varargin,n);

% R1 non-causal filter coefficients
% R2 strictly causal filter coefficients
R1=floor(Re/2);Rc1=floor(Rc/2);
clear Re Rc
R2=R1;Rc2=Rc1;

% unit filter
id=[zeros(1,R1),1,zeros(1,R2)];

% deflated observation vector
x_def=x;
for r=1:N-1
    
if verbose,fprintf('convolutive sphering %d .. \n',r),end
tic

% creation of a multi-delayed observation vector
for k=1:N-r+1
    for i=-R1:R2
    x_tilde((R1+R2+1)*(k-1)+i+1+R1,:)=x_def(k,we1-i:we2-i);
    end
end

% whitening process
[x_prime,B]=white(x_tilde,we2-we1+1);
clear x_tilde;

% initialization of the extraction process
if initializ
    id1=id(ones(1,N-r+1),:);
    id1=reshape(id1',1,(N-r+1)*(R1+R2+1));
    fe0=id1*B^-1;
else
    fe0=ones(1,(N-r+1)*(R1+R2+1))*B^-1;
end

if verbose,fprintf('fixed-point extraction %d .. \n',r),end
% extraction of a non-Gaussian innovation process
[fe{r},iter]=extract(x_prime,fe0,criterion,epsi,iter_max);

te(r)=toc;
clear x_prime;

fe{r}=fe{r}*B;

% reshaping of the extraction vector
fe{r}=reshape1(fe{r},[1,N-r+1,R1+R2+1],[3 2 1]);

% application of the extraction filters
y(r,:)=conv_mix_nc(fe{r},x_def(1:N-r+1,:),R1);


for k=1:N

% coloration of the r^th source for the r^th sensor
if verbose,fprintf('recoloration process %d-%d .. \n',k,r),end
tic
[sc(k,r,:),fc{k,r}]=color(y(r,:),x(k,:),wc1,wc2,Rc1,Rc2);
tc(k,r)=toc;

end

% removing of the last estimated contributions from the sensors
for t=1:N
    s0(1,:)=sc(t,r,:);
    x_def(t,:)=x_def(t,:)-s0;
end
clear s0;

end

% the last contributions corresponding to the N^th source
sc(:,N,:)=x_def;

% identification of the most powerful contribution
[M,Isc]=max(mean(sc.^2,3));

% output signals
for k=1:N
    Sc(k,:)=sc(Isc(k),k,:);
end

%-----------------------------------------
function [x_prime,B]=white(x_tilde,n)

% eigendecomposition of xxt
[E,D]=eig(x_tilde*x_tilde');

% computation of whitening matrix
sqD=sqrt(n)*D^(-1/2);
B=sqD*E';

% whitening process
x_prime=B*x_tilde;

%------------------------------------------------------------
function [w,iter]=extract(x_prime,w,criterion,epsi,iter_max)

[N,n]=size(x_prime);

% extraction vector normalization
w=w./norm(w);

test=1;iter=0;
while test>epsi & iter<iter_max
    
    y=w*x_prime;
    
    switch criterion
        case 'k'
        % kurtosis fixed-point iteration
        w1=1/n*(y.^3)*x_prime'-3*w;
        case 'g'
        % Gaussian negentropic fixed-point iteration
        gauss_y=exp(-y.^2/2);
        w1=1/n*(y.*gauss_y)*x_prime' - 1/n*sum(((1-y.^2).*gauss_y))*w;
        case 't'
        % hyperbolic-tangent negentropic fixed-point iteration
        tanh_y=tanh(y);
        w1=1/n*tanh_y*x_prime' - 1/n*sum(1-tanh_y.^2)*w;
    end
    
    % vector normalization
    w1=w1./norm(w1);
    
    % computation of the stopping criterion
    test=min([abs(w1-w) abs(w1+w)]);
    w=w1;
    
    iter=iter+1;
end

%-------------------------------------------------------
function [sc,c]=color(y,x,w1,w2,R1,R2)

n=w2-w1+1;

% computation of the autocorrelation matrix of y
y2=y(w1+R1:w2+R1);
for r=-R1:R2     
    y1=y(w1-r:w2-r);      
    Ry1(r+R1+1)=1/n*y1*y2';        
end

for r1=-R1:R2
    for r2=-R1:R2
        Ry(r1+R1+1,r2+R1+1)=Ry1(abs(r2-r1)+1);
    end
    % computation of the intercorrelation vector of y and x
    y1=y(w1-r1:w2-r1);
    ryx(r1+R1+1)=1/n*y1*(x(w1:w2))';
end

% computation of the wiener filter
c=Ry^(-1)*ryx';

% non-causal filtering process
sc=filter_nc(c,y,R1);

%-------------------------------------------------
function x=conv_mix_nc(A,s,R1);
% R1 non-causal coefficients
% R2 strictly causal coefficients

R2=size(A,3)-R1-1;

for i=1:size(A,1)
    x(i,:)=zeros(1,size(s,2)+R1+R2);
    for j=1:size(s,1)
        x(i,:)=x(i,:)+conv(reshape(A(i,j,:),size(A,3),1,1),s(j,:));
    end
end

% removing of the out samples
x=x(:,1+R1:end-R2);

%-------------------------------------------------
function s=filter_nc(h,e,R1);

R2=length(h)-R1-1;

% convolution process
s=conv(h,e);

% removing of the out samples
s=s(1+R1:end-R2);

%---------------------------------------------------
function w1=reshape1(w,dims,order)

% ordering of the tensor dimensions
dims1=dims(order);

% reshaping of the vector w
w1=reshape(w,dims1);

% re-ordering of the tensor
w1=ipermute(w1,order);

%-----------------------------------------------------
function [criterion,initializ,epsi,Re,we1,we2,Rc,wc1,wc2,iter_max,verbose]=args(varargin,n)

% ordering of the input arguments

for i=1:length(varargin)
    if length(varargin{i})==1
    switch varargin{i}
        case {'k','g','t'},criterion=varargin{i};
        case 'i',initializ=1;
        case 'v',verbose=1;
        otherwise
            if varargin{i}<1 & varargin{i}>0,epsi=varargin{i};
            else,if ~exist('Re'),Re=varargin{i};
                else,if ~exist('Rc'),Rc=varargin{i};
                     else iter_max=varargin{i};end,end,end  
    end
    else,if ~exist('we1'),we1=varargin{i}(1);we2=varargin{i}(2);
         else,wc1=varargin{i}(1);wc2=varargin{i}(2);end
    end   
end

if ~exist('criterion'), criterion='g';end
if ~exist('initializ'), initializ=0;end
if ~exist('epsi'), epsi=1e-12;end
if ~exist('Re'), Re=80;end
if ~exist('we1'), we1=Re+1;we2=n-Re-1;end
if ~exist('Rc'), Rc=2*Re;end
if ~exist('wc1'), wc1=Rc+1;wc2=n-Rc-1;end
if ~exist('iter_max'), iter_max=1e4;end
if ~exist('verbose'), verbose=0; end

if verbose,criterion,initializ,epsi,Re,extr_window=[we1,we2],Rc,col_window=[wc1,wc2],iter_max,verbose,end