denoise_delayedcoordinates.m

I developed an algorithm for using local ICA in denoising multidimensional data. It uses delay embed

function [sig,rk,components] = denoise_delayedcoordinates(noisesig,k,clustermethod,reducemethod,choosemethod,sortmethod,gamma,verbose);
%delocalize - denoises delayed data
%
% [sig,rk,components] = denoise_delayedcoordinates(noisesig,k,clustermethod,reducemethod,choosemethod,sortmethod,gamma,verbose);
%
%   clustermethod - kmeans|seq                           (default=kmeans)
%   reducemethod  - pca|ica                              (default=pca)
%   choosemethod  = mdl|mdlk|err                         (default=mdl)
%   sortmethod    = var|kurt                             (default=var)
%
%   k                      - max number of partitions
%   gamma                  - (default=1)
%     (gamma<0)            - number of dimensions to use
%     (gamma>0)
%       (choosemethod=mdl) - mdlgamma
%       (choosemethod=err) - portion of signal to discard (0 none, 1 all)
%   verbose                - 0 no messages, 20 all messages (default=0)
%         
%   sig        - denoised signal
%   rk         - list with number of chosen components
%   components - list of the components
%
% 01.06.03 - initial version
% Peter Gruber

error(nargchk(2,8,nargin))     % check no. of input args

[m,samples]=size(noisesig);

if (nargin<3) clustermethod='kmeans'; end
if (nargin<4) reducemethod='pca'; end
if (nargin<5) choosemethod='mdl'; end
if (nargin<6) sortmethod='kurt'; end
if (nargin<7) gamma=1; end
if (nargin<8) verbose=0; end

[noisesig,meanvalues]=remmean(noisesig);
[localsig,parts,permut]=localize(noisesig,k,m,clustermethod);

if (verbose>8) plot(localsig'); title('delayed signal'); pause; end

ncomp=size(parts,1);
if (verbose>1) warning(sprintf('clusterd into %d parts',ncomp)); end

rk=zeros(1,ncomp);
if (nargout>2) orglocalsig=localsig; end

for i=1:ncomp
 [lsig,lmean]=remmean(localsig(:,parts(i,1):parts(i,2)));

 [V,W,D]=sortcomponents(lsig,reducemethod,sortmethod,verbose);

 p=choosenumberofcomponents(V*lsig,D,choosemethod,gamma,verbose);

 r=size(V,1);
 if (verbose>1) warning(sprintf('Reduction Step (%d/%d/%d)',p,r,m)); end
 
 R=eye(r);
 R(:,p+1:r)=0;

 if (verbose>5) scplot(V*lsig,R*V*lsig); title('reduced Signal'); pause; end

 rk(i)=p;
 
 if (verbose>10) scplot(lsig); title('local delayed Signal'); pause; end

 if (nargout>2)
  for j=1:r
   tempR=eye(r);
   tempR(j,j)=0;
   tempdelayedsig=orglocalsig;
   tempdelayedsig(:,parts(i,1):parts(i,2))=W*tempR*V*(lsig+lmean*ones(1,size(lsig,2)));
   tempdelayedsig=orglocalsig-tempdelayedsig;
   components{(i-1)*m+j}=delocalize(tempdelayedsig,permut);
  end
 end
 lsig=W*R*V*lsig;

 if (verbose>10) scplot(lsig); title('denoised delayed Signal'); pause; end
  
 localsig(:,parts(i,1):parts(i,2))=lsig+lmean*ones(1,size(lsig,2));
 
 if (verbose>8) plot(delocalize(localsig,permut)'); title('denoised Signal'); pause; end
end

sig=delocalize(localsig,permut)+meanvalues*ones(1,samples);

return


function [V,W,D] = sortcomponents(signal,reducemethod,ratingmethod,verbose);
%sortcomponents - find the components
%
% [V,W,D] = sortcomponents(signal,reducemethod,ratingmathod,verbose); 
%
%   reducemethod  = pca|{min,max,comb,auto,var}{rca,ica}    (default=pca)
%
% 01.08.03 - initial version
% 01.11.04 - rewrite
% Peter Gruber

error(nargchk(1,4,nargin))     % check no. of input args

[m,samples]=size(signal);

if (nargin<2) reducemethod='pca'; end
if (nargin<3) ratingmethod='var'; end
if (nargin<4) verbose=0; end

switch reducemethod
 case 'pca'
  [W,D]=eig(cov(signal'));
  V=W';
  ksignal=V*signal;
 case 'mpencil'
  [W,D]=eig(cov(signal'));
  E=D;E(E>0)=1./sqrt(E(E>0));
  V=W*E*W';
  c=max(10,ceil(samples/1000));
  c=exp(-[-c:1:c].^2)/sqrt(2*pi);
  [P,Q]=eig(cov(matconv(V*signal,c)'));
  E=D;E(E>0)=sqrt(E(E>0));
  W=W*E*W'*P;
  W=W./(ones(m,1)*sqrt(sum(W.^2)));
  W=ort(W);
  V=W^-1;
  ksignal=V*signal;
 case 'fastica'
  ni=max(100,size(signal,2)/10);
  if (verbose > 2)
   [ksignal,W,V]=fastica(signal,'approach','symm','verbose','on','g','tanh','maxNumIterations',ni,'stabilization','on');
  else
   [ksignal,W,V]=fastica(signal,'approach','symm','verbose','off','g','tanh','maxNumIterations',ni,'stabilization','on');
  end
  W=W./(ones(m,1)*sqrt(sum(W.^2)));
  W=ort(W);
  V=W^-1;
  ksignal=V*signal;
 case 'pearsonica'
  [ksignal,W,V]=pearson_ica(signal);
  W=W./(ones(m,1)*sqrt(sum(W.^2)));
  W=ort(W);
  V=W^-1;
  ksignal=V*signal;
 case 'jadeica'
  [W,ksignal]=jade(signal);
  W=W./(ones(m,1)*sqrt(sum(W.^2)));
  W=ort(W);
  V=W^-1;
  ksignal=V*signal;
 case 'mdlica'
  [E,D]=eig(cov(signal'));
  D=sum(D,1); D=D(m:-1:1); E=E(:,m:-1:1);
  p=mdl(D,samples,1);
  ksignal=E'*signal;
  if (p==1)
   ksignal=ksignal(1,:);
   V=1; W=1;
  else
   [ksignal,W,V]=pearson_ica(ksignal(1:p,:));
  end
  V=[V,zeros(p,m-p)]*E';
  W=E*[W;zeros(p,m-p)];
  W=W./(ones(m,1)*sqrt(sum(W.^2)));
  V=W^-1;
  ksignal=V*signal;
end

k=size(ksignal,1);
R=eye(k);
R=R(k:-1:1,:);

switch ratingmethod
 case 'kurt'
  [dummy,kp]=sort(kurtosis(ksignal').^2);
  [dummy2,kp]=sort(kp);
 case 'maxkurt'
  [dummy,kp]=sort(-kurtosis(ksignal'));
  [dummy2,kp]=sort(kp);
 case 'minkurt'
  [dummy,kp]=sort(kurtosis(ksignal'));
  [dummy2,kp]=sort(kp);
 case 'comb'
  [dummy,kp]=sort((var(ksignal').*kurtosis(ksignal')).^2);
  [dummy2,kp]=sort(kp);
 case 'auto'
  [dummy,kp]=sort(var(autocor(ksignal')));
  [dummy2,kp]=sort(kp);
 otherwise
  [dummy,kp]=sort(var(ksignal'));
  [dummy2,kp]=sort(kp);
end
if (verbose>2) plot(dummy); title(['Sortcriterion: ',ratingmethod]); pause; end
R=R(:,kp);
V=R*V;
W=W*R';

% the rows contain eigenvalues of the covariances of the
% signal and the estimated noise
switch reducemethod
 case 'pca'
  D=diag(var((R*ksignal)'));
  D=ones(m,1)*diag(D)';
 otherwise
  T=W;
  C=cov((R*ksignal)');
  D(1,:)=diag(C)';
  for i=1:m
   S=ort(V(i+1:m,:));
   S=ort(S(m-i+1:m,:));
   d=S*T*C*T'*S';
   D(i,1:i)=eig(d(1:i,1:i))';
   D(i,i+1:m)=eig(d(i+1:m,i+1:m))';
  end
end

return

function N=matconv(M,C)
[d,m]=size(M);
c=length(C);
N=zeros(d,m+c-1);
for i=1:size(M,1)
 N(i,:)=conv(M(i,:),C(:));
end
return
endfunction



function [localdata,parts,permutation] = localize(sig,k,m,method);
%localize - generate localized data
%
% [localdata,parts,permutation] = localize(sig,k,m,method);
%
%  method - {kmeans,seq}
%
% 1.8.03
% Peter Gruber

error(nargchk(2,4,nargin))     % check no. of input args

if (nargin<3) m=0; end
if (nargin<4) method='kmeans'; end

[dim,samples]=size(sig);

if (k==1) localdata=sig; parts=[1,samples]; permutation=[1:samples]; return; end

if (m>samples) m=samples; end

switch method
 case 'kmeans'; [cent,cl]=kmeans_cluster('batch',sig',k);cent=cent';
 case 'fast';   [cent,cl]=cluster_fast(sig,k,max(500,ceil(samples/500)));
 case 'seq';    [cent,cl]=cluster(sig,k);
 case 'diadic'; [cent,cl]=cluster_diadic(sig,k);
 case 'simple'; [cent,cl]=cluster(sig,k,0);
end

p=zeros(size(cent,2),1); for i=1:size(cent,2) p(i)=sum(cl==i); end

while ((min(p)<m) || (sum(isfinite(p))>k))
 [i,im]=min(p);
 q=zeros(size(cent,2),1);
 for i=1:size(cent,2) 
  if (isfinite(p(i))) q(i)=dist(cent(:,im),cent(:,i));
  else q(i)=+Inf; end
 end
 [i,j]=sort(q);
 cl=cl+(cl==j(1))*(j(2)-j(1));
 p(j(2))=p(j(2))+p(j(1));
 p(j(1))=+Inf;
end

k=sum(isfinite(p));
pt=zeros(k,3);
[i,j]=sort(p);
pt(:,3)=cumsum(i(1:k));
pt(:,1)=j(1:k);
pt(1,2)=1; pt(2:k,2)=pt(1:k-1,3)+1;

i=zeros(samples,1);
for j=1:k i=i+(cl==pt(j,1))*(j-pt(j,1)); end
cl=cl+i; clear i

localdata=zeros(dim,samples);
permutation=zeros(samples,1);

nn=[1:samples];
for i=1:k
 localdata(:,pt(i,2):pt(i,3))=sig(:,cl==i);
 permutation(pt(i,2):pt(i,3))=nn(cl==i);
end

parts=pt(:,2:3);
return



function sig = delocalize(localsig,permutation);
%delocalize - regenerate localized data
%
% sig = delocalize(localsig,permutation);
%
% 1.6.03
% Peter Gruber

error(nargchk(2,2,nargin))     % check no. of input args

[dummy,p]=sort(permutation);
sig=localsig(:,p);
return



function [centroids,clusters,err] = cluster(D, n, epochs, verbose);
% cluster - find clusters
%
% [centroids,clusters,err] = cluster(D, n, epochs, verbose)
%
% 1.6.03
% Peter Gruber

error(nargchk(2,4,nargin))     % check no. of input args

if (nargin<3) epochs=100; end
if (nargin<4) verbose=0; end

[dim,samples] = size(D);

if (n==1)
 centroids=sum(D,2)/samples;
 clusters=ones(samples,1);
 err=sum(sqrt(sum((D-centroids*ones(1,samples)).^2,1)))/samples;
 return
end

partition=zeros(1,n+1); centers=zeros(1,n);

for i=2:n partition(i)=floor(samples/n)*(i-1); end
partition(n+1)=samples;

for t=1:epochs
 oldpart=partition;
 for i=1:n centers(i)=partition(i)+ceil((partition(i+1)-partition(i))/2); end
 for i=1:n-1
  num=centers(i+1)-centers(i)+1;
  if (num==2) partition(i+1)=centers(i);
  else
   dists=zeros(1,num);
   for j=1:num-1 dists(j+1)=dists(j)+sqrt(sum((D(:,centers(i)+j-1)-D(:,centers(i)+j)).^2)); end
   partition(i+1)=centers(i)+sum(dists<(dists(num)/2));
  end
 end
 if all(partition==oldpart) break; end
end

centroids=zeros(dim,n); clusters=zeros(samples,1); err=0;
for i=2:n+1
 s=partition(i-1)+1;
 e=partition(i); n=e-s+1;
 centroids(:,i-1)=sum(D(:,s:e),2)/n;
 clusters(s:e)=ones(1,n)*(i-1);
 err=err+sum(sqrt(sum((D(:,s:e)-centroids(:,i-1)*ones(1,n)).^2,1)));
end

err=err/samples;
return



function [centroids,clusters,err] = cluster_diadic(D, n, verbose);
% cluster_diadic - find clusters
%
% [centroids,clusters,err] = cluster(D, n, verbose)
%
% 1.10.03
% Peter Gruber

error(nargchk(2,3,nargin))     % check no. of input args

if (nargin<3) verbose=0; end

[dim,samples] = size(D);

if (n==1)
 centroids=sum(D,2)/samples;
 clusters=ones(samples,1);
 err=sum(sqrt(sum((D-centroids*ones(1,samples)).^2,1)))/samples;
 return
end

clusters=ones(samples,1);

m=1; while ((2^m<n) && (m<dim)) m=m+1; end

Dm = D - sum(D,2) * ones (1,samples) / samples;
[E,C]=eig(cov(Dm'));
C=ones(1,dim)*C;
Dm=E'*Dm;

[s,p]=sort(C);

for j=1:m
 d=Dm(p(j),:);
 clusters(d>0)=clusters(d>0)+2^(j-1);
end

centroids=zeros(dim,2^m); err=0;
for i=1:2^m
 if (sum(clusters==i)==0)
  centroids(:,i)=zeros(dim,1);
 else
  centroids(:,i)=sum(D(:,clusters==i),2)/sum(clusters==i);
  err=err+sum(sqrt(sum((D(:,clusters==i)-centroids(:,i)*ones(1,sum(clusters==i))).^2,1)));
 end
end

err=err/samples;
return



function [centroids,clusters,err] = cluster_fast(D, n, numkmeans, verbose);
% cluster_fast - find clusters
%
% [centroids,clusters,err] = cluster(D, n, numkmeans, verbose)
%
% 1.10.03
% Peter Gruber

error(nargchk(2,4,nargin))     % check no. of input args

if (nargin<3) numkmeans=500; end
if (nargin<4) verbose=0; end

[dim,samples] = size(D);

numkmeans=max(numkmeans,n);

sig=randperm(samples);
sig=D(:,sig);
sig=sig(:,1:numkmeans);
[centroids,cl]=kmeans_cluster('batch',sig',n);centroids=centroids';

distcl=ones(1,samples)*Inf;
clusters=zeros(samples,1);
for i=1:size(centroids,2)
 d=sum((D-centroids(:,i)*ones(1,samples)).^2,1);
 t=d-distcl;
 clusters(t<0)=i;
 distcl=min(d,distcl);
end

err=0;
for i=1:size(centroids,2)
 d=D(:,clusters==i);
 err=err+sum(sqrt(sum((d-centroids(:,i)*ones(1,size(d,2))).^2,1)));
end

err=err/samples;

return