denoise_locca.m

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

function [sig,err,snr] = denoise_locca(noisesig,mrange,krange,delaymethod,delayspeed,clustermethod,reducemethod,undelaymethod,sortmethod,choosemethod,ratingmethod,gamma,repeat,verbose,orig);
%denoise_locca - denoise with local component analysis
%
% [sig,edl,err,snr] = denoise_locca(noisesig,mrange,krange,delaymethod,delayspeed,clustermethod,reducemethod,undelaymethod,choosemethod,ratingmethod,gamma,repeat,verbose,orig)
%
%   delaymethod   = reshape|multidim|multidimflip        (default=multidim)
%   clustermethod = kmeans|seq                           (default=kmeans)
%   reducemethod  = pca|ica|mdlica                       (default=pca)
%   undelaymethod = stat|avg|quad                        (default=avg)
%   sortmethod    = var|kurt|autovar                     (default=var)
%   choosemethod  = mdl|err|cluster                      (default=mdl)
%   ratingmethod  = mdl|err                              (default=mdl)
%
%   mrange                 - array of dimensions
%   krange                 - array of partitions
%   delayspeed             - lag per dimension (default=1)
%   gamma                  - (default=1)
%     (gamma<0)            - number of dimensions to use
%     (gamma>0)
%       (choosemethod=mdl) - mdlgamma (1 parsimonious, 64 less restrictive)
%       (choosemethod=err) - portion of signal to discard (0 none, 1 all)
%   repeat                 - how many time to repeat the reduction step (default=1)
%                            (repeat-floor(repeat)) mixing of the signals (0 only old signal, 1 only new signal)
%   verbose                - 0 no messages, 20 all messages (default=0)
%   orig                   - original signal to compare (for SNR)
%         
%   sig - denoise signal
%   edl - MDLs of the computed noises
%   err - strength of the computed noises
%   snr - SNRs of the computed signals
%  
% 01.08.03 - initial version
% 01.11.04 - reorganistaion
% Peter Gruber

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

if (nargin<4)  delaymethod='multidim'; end
if (nargin<5)  delayspeed=1; end
if (nargin<6)  clustermethod='kmeans'; end
if (nargin<7)  reducemethod='pca'; end
if (nargin<8)  undelaymethod='avg'; end
if (nargin<9)  sortmethod='var'; end
if (nargin<10) choosemethod='mdl'; end
if (nargin<11) ratingmethod='err'; end
if (nargin<12) gamma=1; end
if (nargin<13) repeat=1; end
if (nargin<14) verbose=0; end
if (nargin<15) orig=0; end

sigsize=size(noisesig);

switch delaymethod
 case 'reshape'
  orgsigsize=sigsize;
  noisesig=reshape(noisesig,1,prod(orgsigsize));
  sigsize=size(noisesig);
 otherwise
  orgsigsize=sigsize;
end

if (verbose>0) warning('Signal %dx%d\n',sigsize(1),sigsize(2)); end

err=zeros(size(mrange,1),size(krange,1));
snr=zeros(size(mrange,1),size(krange,1));

G=Inf;

M=max(mrange);

if (prod(size(orig))==prod(size(noisesig)))
 origmdl=getmdl(orig-noisesig,M,delaymethod,delayspeed);
end

for mm=1:length(mrange) m=mrange(mm);

 [savedsig,savedmv]=remmean(delay(noisesig,m,delayspeed,delaymethod));

 for kk=1:length(krange) k=krange(kk);

  for rcount=1:floor(repeat)
   if (verbose>0) & (rcount<2) warning(sprintf('Denoising with parameters (m,k)=(%3.1f,%d) ...',m,k)); end

   if (rcount>1)
    rmix=repeat-floor(repeat);
    [delayedsig,mv]=remmean(delay((1-rmix)*noisesig+(rmix)*newsig,m,delayspeed,delaymethod));
   else
    delayedsig=savedsig; mv=savedmv;
   end
   
   [delayedsig,comp]=denoise_delayedcoordinates(delayedsig,k,clustermethod,reducemethod,choosemethod,sortmethod,gamma,verbose);
   if (verbose>1) warning([sprintf('Number of chosen components: ') sprintf('%d ',comp)]); end

   newsig=undelay(delayedsig+mv*ones(1,size(delayedsig,2)),sigsize,m,delayspeed,delaymethod,undelaymethod);
  end
  
  switch ratingmethod
   case 'err'
    err(mm,kk)=-dist(noisesig,newsig);
   case 'mdl'
    err(mm,kk)=getmdl(newsig-noisesig,M,delaymethod,delayspeed);
  end
  
  if (prod(size(orig))==prod(size(noisesig)))
   snr(mm,kk)=SNR(orig,newsig);
   if (verbose>1)
    warning(sprintf('Criterion: %e SNR: %f (noisy %f) MDL (real noise: %e)',err(mm,kk),snr(mm,kk),SNR(orig,noisesig),origmdl));
    if (verbose>8) cplot('combined',noisesig,newsig,noisesig-newsig,orig-newsig); title('Criterion'); pause; end
   end
  else
   if (verbose>1)
    warning(sprintf('Criterion: %e',err(mm,kk)));
    if (verbose>8) cplot('combined',noisesig,newsig,noisesig-newsig); title('Criterion'); pause; end
   end
  end

  if (err(mm,kk)<=G) G=err(mm,kk); result=newsig; end
 end
end

criterion=err;

if (size(criterion,1)==1)
 [t,k]=min(criterion); m=1;
else
 if (size(criterion,2)==1) [t,m]=min(criterion); k=1;
 else [t,m]=min(criterion); [dummy,k]=min(t); m=m(k);
 end
end
mm=mrange(m); kk=krange(k);

if (verbose>0) warning(sprintf('Optimal for m=%f k=%d',mm,kk)); end

switch delaymethod
 case 'reshape'
  sig=reshape(result,orgsigsize);
 otherwise
  sig=result;
end

return