denoise_locca_findbest.m

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

function [sig,edl,err,snr] = denoise_locca_findbest(noisesig,maxm,maxk,steps,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_findbest(noisesig,maxm,maxk,steps,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)
%   ratingmethod  = var|kurt|autovar                     (default=var)
%   choosemethod  = mdl|err|cluster                      (default=mdl)
%   ratingmethod  = mdl|err                              (default=mdl)
%
%   maxm                   - max of dimensions
%   maxk                   - max of partitions
%   steps                  - use steps steps
%   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(4,16,nargin))     % check no. of input args

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

sigsize=size(noisesig);

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

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

G=Inf;
currentm=[1,maxm]; currentk=[1,maxk];

for step=1:steps
 if (verbose>0) warning(sprintf('Step %d m=%d,%d k=%d,%d',step,currentm(1),currentm(2),currentk(1),currentk(2))); end

 for mm=1:2; m=currentm(mm);

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

  for kk=1:2; k=currentk(kk);
   
   if (err(m,k)==Inf)

    for rcount=1:floor(repeat)
     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
   
     if (verbose>0) & (rcount<2) warning(sprintf('Trying to denoise with parameters (m,k)=(%d,%d) ... ',m,k)); 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
  
    noise=newsig-noisesig;
    M=max(maxm);
    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(m,k)=SNR(orig,newsig);
     if (verbose>1)
      warning(sprintf('Criterion: %e SNR: %f (noisy %f) MDL (real noise: %e)',err(mm,kk),SNR(orig,newsig),SNR(orig,noisesig),origmdl));
      if (verbose>3) cplot('combined',noisesig,newsig,noisesig-newsig,orig-newsig); title('Noise'); pause; end
     end
    else
     if (verbose>1)
      warning(sprintf('Criterion: %e',err(mm,kk));
      if (verbose>3) cplot('combined',noisesig,newsig,noisesig-newsig); title('Noise'); pause; end
     end
    end

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

    if (verbose>0) warning('... done'); end 
   else
    if (verbose>0) warning(sprintf('Trying to denoise with parameters (m,k)=(%d,%d) already done',m,k)); end
   end
  end
 end

 crit=err(currentm,currentk);
 crit=(crit==min(min(crit))).*[1,2;3,4];
 crit=sum(crit(:));
 if (crit>2); currentm(1)=sum(currentm)/2+0.5; crit=crit-2;
 else;        currentm(2)=sum(currentm)/2-0.5; end
 if (crit>1); currentk(1)=sum(currentk)/2+0.5;
 else;        currentk(2)=sum(currentk)/2-0.5; end

 if ((currentm(1)==currentm(2)) && (currentk(1)==currentk(2))) break; 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

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

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

return