cfica.m 7.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288
  1. % C-FICA algorithm - time-domain Convolutive extension of FastICA - v1.2
  2. %
  3. % J. Thomas, Y. Deville and S. Hosseini
  4. % "Time-domain fast fixed-point algorithms for convolutive ICA",
  5. % IEEE Signal Processing Letters, vol. 13, no. 4, pp. 228-231, April 2006.
  6. %
  7. %
  8. % [Sc,sc,Isc,y,fe,fc,te,tc]=cfica(x,non-lin,epsi,Re,[we1 we2],Rc,[wc1 wc2],itermax,initializ,verbose)
  9. %
  10. % Inputs (place the deconvolution parameters before the recoloration ones
  11. % themselves before iter_max)
  12. %
  13. % x : observation vector
  14. % non-lin : 'k'->kurtosis, 'g'->Gaussian, 't'->tanh default 'g'
  15. % initializ : 'i' for initialization with unit filters optional
  16. % epsi : stopping criterion (between 0 and 1) default 10^-12
  17. % Re : order of the extraction filters default 80
  18. % [we1 we2] : extraction time window optional
  19. % Rc : order of the recoloration filters default 3R1
  20. % [wc1 wc2] : recoloration time window optional
  21. % itermax : maximum number of iteration default 10000
  22. % verbose : 'v' for verbose optional
  23. %
  24. % Outputs
  25. %
  26. % Sc : most powerful contributions for each estimated source
  27. % sc : estimated source contributions (tensor N*N*nb_samples)
  28. % Isc : indices of the most powerful contributions
  29. % y : N-1 estimated innovation processes (matrix (N-1)*nb_samples)
  30. % fe : extraction filters
  31. % fc : recoloration filters
  32. % te : times of extraction
  33. % tc : times of recoloration
  34. % -------------------------------------------------------------------------
  35. % History of the software:
  36. % V1.1: original software written by J. Thomas
  37. %
  38. % V1.2: Error messages appeared with recent releases of Matlab in the function
  39. % conv_mix_nc. This version corrects them.
  40. % Author: Matthieu PUIGT Contact: matthieu.puigt@univ-littoral.fr
  41. %
  42. % -------------------------------------------------------------------------
  43. function [Sc,sc,Isc,y,fe,fc,te,tc]=cfica(x,varargin)
  44. % N: number of observations, n: number of samples
  45. [N,n]=size(x);
  46. % centering of the observations
  47. x=x-mean(x,2)*ones(1,n);
  48. % ordering of the input arguments
  49. [criterion,initializ,epsi,Re,we1,we2,Rc,wc1,wc2,iter_max,verbose]=args(varargin,n);
  50. % R1 non-causal filter coefficients
  51. % R2 strictly causal filter coefficients
  52. R1=floor(Re/2);Rc1=floor(Rc/2);
  53. clear Re Rc
  54. R2=R1;Rc2=Rc1;
  55. % unit filter
  56. id=[zeros(1,R1),1,zeros(1,R2)];
  57. % deflated observation vector
  58. x_def=x;
  59. for r=1:N-1
  60. if verbose,fprintf('convolutive sphering %d .. \n',r),end
  61. tic
  62. % creation of a multi-delayed observation vector
  63. for k=1:N-r+1
  64. for i=-R1:R2
  65. x_tilde((R1+R2+1)*(k-1)+i+1+R1,:)=x_def(k,we1-i:we2-i);
  66. end
  67. end
  68. % whitening process
  69. [x_prime,B]=white(x_tilde,we2-we1+1);
  70. clear x_tilde;
  71. % initialization of the extraction process
  72. if initializ
  73. id1=id(ones(1,N-r+1),:);
  74. id1=reshape(id1',1,(N-r+1)*(R1+R2+1));
  75. fe0=id1*B^-1;
  76. else
  77. fe0=ones(1,(N-r+1)*(R1+R2+1))*B^-1;
  78. end
  79. if verbose,fprintf('fixed-point extraction %d .. \n',r),end
  80. % extraction of a non-Gaussian innovation process
  81. [fe{r},iter]=extract(x_prime,fe0,criterion,epsi,iter_max);
  82. te(r)=toc;
  83. clear x_prime;
  84. fe{r}=fe{r}*B;
  85. % reshaping of the extraction vector
  86. fe{r}=reshape1(fe{r},[1,N-r+1,R1+R2+1],[3 2 1]);
  87. % application of the extraction filters
  88. y(r,:)=conv_mix_nc(fe{r},x_def(1:N-r+1,:),R1);
  89. for k=1:N
  90. % coloration of the r^th source for the r^th sensor
  91. if verbose,fprintf('recoloration process %d-%d .. \n',k,r),end
  92. tic
  93. [sc(k,r,:),fc{k,r}]=color(y(r,:),x(k,:),wc1,wc2,Rc1,Rc2);
  94. tc(k,r)=toc;
  95. end
  96. % removing of the last estimated contributions from the sensors
  97. for t=1:N
  98. s0(1,:)=sc(t,r,:);
  99. x_def(t,:)=x_def(t,:)-s0;
  100. end
  101. clear s0;
  102. end
  103. % the last contributions corresponding to the N^th source
  104. sc(:,N,:)=x_def;
  105. % identification of the most powerful contribution
  106. [M,Isc]=max(mean(sc.^2,3));
  107. % output signals
  108. for k=1:N
  109. Sc(k,:)=sc(Isc(k),k,:);
  110. end
  111. %-----------------------------------------
  112. function [x_prime,B]=white(x_tilde,n)
  113. % eigendecomposition of xxt
  114. [E,D]=eig(x_tilde*x_tilde');
  115. % computation of whitening matrix
  116. sqD=sqrt(n)*D^(-1/2);
  117. B=sqD*E';
  118. % whitening process
  119. x_prime=B*x_tilde;
  120. %------------------------------------------------------------
  121. function [w,iter]=extract(x_prime,w,criterion,epsi,iter_max)
  122. [N,n]=size(x_prime);
  123. % extraction vector normalization
  124. w=w./norm(w);
  125. test=1;iter=0;
  126. while test>epsi & iter<iter_max
  127. y=w*x_prime;
  128. switch criterion
  129. case 'k'
  130. % kurtosis fixed-point iteration
  131. w1=1/n*(y.^3)*x_prime'-3*w;
  132. case 'g'
  133. % Gaussian negentropic fixed-point iteration
  134. gauss_y=exp(-y.^2/2);
  135. w1=1/n*(y.*gauss_y)*x_prime' - 1/n*sum(((1-y.^2).*gauss_y))*w;
  136. case 't'
  137. % hyperbolic-tangent negentropic fixed-point iteration
  138. tanh_y=tanh(y);
  139. w1=1/n*tanh_y*x_prime' - 1/n*sum(1-tanh_y.^2)*w;
  140. end
  141. % vector normalization
  142. w1=w1./norm(w1);
  143. % computation of the stopping criterion
  144. test=min([abs(w1-w) abs(w1+w)]);
  145. w=w1;
  146. iter=iter+1;
  147. end
  148. %-------------------------------------------------------
  149. function [sc,c]=color(y,x,w1,w2,R1,R2)
  150. n=w2-w1+1;
  151. % computation of the autocorrelation matrix of y
  152. y2=y(w1+R1:w2+R1);
  153. for r=-R1:R2
  154. y1=y(w1-r:w2-r);
  155. Ry1(r+R1+1)=1/n*y1*y2';
  156. end
  157. for r1=-R1:R2
  158. for r2=-R1:R2
  159. Ry(r1+R1+1,r2+R1+1)=Ry1(abs(r2-r1)+1);
  160. end
  161. % computation of the intercorrelation vector of y and x
  162. y1=y(w1-r1:w2-r1);
  163. ryx(r1+R1+1)=1/n*y1*(x(w1:w2))';
  164. end
  165. % computation of the wiener filter
  166. c=Ry^(-1)*ryx';
  167. % non-causal filtering process
  168. sc=filter_nc(c,y,R1);
  169. %-------------------------------------------------
  170. function x=conv_mix_nc(A,s,R1);
  171. % R1 non-causal coefficients
  172. % R2 strictly causal coefficients
  173. R2=size(A,3)-R1-1;
  174. for i=1:size(A,1)
  175. x(i,:)=zeros(1,size(s,2)+R1+R2);
  176. for j=1:size(s,1)
  177. x(i,:)=x(i,:)+conv(reshape(A(i,j,:),size(A,3),1,1),s(j,:));
  178. end
  179. end
  180. % removing of the out samples
  181. x=x(:,1+R1:end-R2);
  182. %-------------------------------------------------
  183. function s=filter_nc(h,e,R1);
  184. R2=length(h)-R1-1;
  185. % convolution process
  186. s=conv(h,e);
  187. % removing of the out samples
  188. s=s(1+R1:end-R2);
  189. %---------------------------------------------------
  190. function w1=reshape1(w,dims,order)
  191. % ordering of the tensor dimensions
  192. dims1=dims(order);
  193. % reshaping of the vector w
  194. w1=reshape(w,dims1);
  195. % re-ordering of the tensor
  196. w1=ipermute(w1,order);
  197. %-----------------------------------------------------
  198. function [criterion,initializ,epsi,Re,we1,we2,Rc,wc1,wc2,iter_max,verbose]=args(varargin,n)
  199. % ordering of the input arguments
  200. for i=1:length(varargin)
  201. if length(varargin{i})==1
  202. switch varargin{i}
  203. case {'k','g','t'},criterion=varargin{i};
  204. case 'i',initializ=1;
  205. case 'v',verbose=1;
  206. otherwise
  207. if varargin{i}<1 & varargin{i}>0,epsi=varargin{i};
  208. else,if ~exist('Re'),Re=varargin{i};
  209. else,if ~exist('Rc'),Rc=varargin{i};
  210. else iter_max=varargin{i};end,end,end
  211. end
  212. else,if ~exist('we1'),we1=varargin{i}(1);we2=varargin{i}(2);
  213. else,wc1=varargin{i}(1);wc2=varargin{i}(2);end
  214. end
  215. end
  216. if ~exist('criterion'), criterion='g';end
  217. if ~exist('initializ'), initializ=0;end
  218. if ~exist('epsi'), epsi=1e-12;end
  219. if ~exist('Re'), Re=80;end
  220. if ~exist('we1'), we1=Re+1;we2=n-Re-1;end
  221. if ~exist('Rc'), Rc=2*Re;end
  222. if ~exist('wc1'), wc1=Rc+1;wc2=n-Rc-1;end
  223. if ~exist('iter_max'), iter_max=1e4;end
  224. if ~exist('verbose'), verbose=0; end
  225. if verbose,criterion,initializ,epsi,Re,extr_window=[we1,we2],Rc,col_window=[wc1,wc2],iter_max,verbose,end