Mise à jour de 'cfica.m'

cfica.m

@@ -1,288 +1,288 @@
-% 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@lisic.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
+% 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