denoise.m

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

function [Xdenoised, criterion]=denoise(X,varargin)
%
% local projective denoising
%
% Implematation of a denoising algorithms for enhancing noisy signals based on
% LICA, dAMUSE and KPCA. The algorithm LICA relies on applying ICA locally to
% clusters of signals embedded in a high dimensional feature space of delayed
% coordinates. The components resembling the signals can be detected by various
% criteria like estimators of kurtosis or the variance of autocorrelations
% depending on the statistical nature of the signal. The algorithm proposed can
% be applied favorably to the problem of denoising multidimensional data.
%
% Demo code using this can be found in demo.m
%
% usage:
%       [Xdenoised, criterion]=denoise(X [,argID,value [...]])
%
% with input and output arguments ([]'s are optional):
% X                 (matrix) data matrix (any dimension)
% [argID,           (string) Additional parameters are given as argID, value
%   value]          (varies) pairs. See below for list of valid IDs and values.
% Xdenoised         (matrix) denoised data matrix X
% criterion         (matrix) if multiple paramters are given, gives the criterion upon which the selection was based
%
% Here are the valid argument IDs and corresponding values. 
%  'delaydim' or    (integer) delay dimension (if dimension of X>1 only approximatively) 
%  'delaydim'       (vector)  try any of these and use the one yielding the best denoising
%  'cluster' or     (integer) number of clusters to use
%  'cluster'        (vector)  try any of these and use the one yielding the best denoising
%  'search'         (integer) !=0: search for the best M, K up to the maxima given by M, K
%  'method'         (string)  'ica': use ica to separate noise from the signal
%                             'pca': use pca to separate noise from the signal
%  'rateby'         (string)  criterion to identify the signal in the delayed mixture
%                             'var' : strongest most likely to be signal
%                             'kurt': largest (abs) kurtosis most lekely to be signal
%  'gamma'          (scalar)  gamma for the mdl estimators
%  'repeat'         (integer) reapeat the denoising process
%  'delaystep'      (scalar)  delaystep (in timesteps)
%  'verbose'        (integer) report progress of algorithm in text format. (higher number = more verbosity, 0 = no text output)
%
% Author : Peter Gruber <[e-mail]>
% This software is for non-commercial use only.
%
% ------------------------------------------------------------------------------------------------
% REFERENCES:
%
% P. Gruber, K. Stadlthanner, A.M. Tomé, A.R. Teixeira, F.J. Theis, C.G.
% Puntonet, and E.W. Lang. Denoising using local ICA and a generalized
% eigendecomposition with time-delayed signals. In Proc. ICA 2004, volume 3195
% of Lecture Notes in Computer Science, pages 992-1000, Granada, Spain, 2004.
%
% P. Gruber, F.J. Theis, K. Stadlthanner, E.W. Lang, A.M. Tomé, and A.R.
% Teixeira. Denoising using local ICA and kernel-PCA. In Proc. IJCNN 2004, pages
% 2071-2076, Budapest, Hungary, 2004.
%
% P. Gruber, F.J. Theis, A.M. Tomé, and E.W. Lang. Automatic denoising using
% local independent component analysis. In Proc. EIS 2004, Madeira, Portugal,
% 2004.
%

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

M             = 35;
K             = max(10,length(X(:))/5000);
delaymethod   = 'multidim';
delayspeed    = 2;
clustermethod = 'fast';
reducemethod  = 'mpencil';
undelaymethod = 'avg';
sortmethod    = 'comb';
choosemethod  = 'mdl';
ratingmethod  = 'err';
gamma         = 1;
repeat        = 1;
search        = 0;
verbose       = 0;

for i=1:2:length(varargin),
    %% Check that all argument types are strings
    if ~ischar(varargin{i}),
        error('Invalid input identifier names or input argument order.');
    end
  
    %% Lower/uppercase in identifier types doesn't matter: 
    identifier=lower(varargin{i});     % identifier (lowercase)
    value=varargin{i+1};
        
    switch identifier        
    case 'verbose'
	if ( isscalar(value) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        verbose = value;
    case 'delaydim'
	if ( isvector(value) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        M = value(:);
    case 'cluster'
	if ( isvector(value) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        K = value(:);
    case 'method'
	if ( all(isletter(value)) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        reducemethod = lower(value);
	switch reducemethod
	case 'ica'
		reducemethod  = 'pearsonica'
		sortmethod    = 'kurt';
	case 'fastica'
		sortmethod    = 'kurt';
	case 'pearsonica'
		sortmethod    = 'kurt';
	case 'jadeica'
		sortmethod    = 'kurt';
	case 'mpencil'
		sortmethod    = 'kurt';
	case 'pca'
		sortmethod    = 'var';
	end
    case 'sortby'
	if ( all(isletter(value)) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        sortmethod = lower(value);
    case 'chooseby'
	if ( all(isletter(value)) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        choosemethod = lower(value);
    case 'rateby'
	if ( all(isletter(value)) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        ratingmethod = lower(value);
    case 'gamma'
	if ( isscalar(value) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
	gamma = value;
    case 'search'
	if ( isscalar(value) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
        if ( value~=0 )
	    M=max(M(:));
	    K=max(K(:));
	    search = 1;
	end
    case 'repeat'
	if ( isscalar(value) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
	repeat = floor(value)+0.7;
    case 'delaystep'
	if ( isscalar(value) == 0 ) error(sprintf('Invalid value for identifier %s!',identifier)); end
	delayspeed = ceil(value);
    otherwise
        %%% Unknown identifier
        error(sprintf('Invalid argument identifier %s!',identifier));  
    end
end

d=size(X); d=sum(d~=1);
M = 0.5*M.^(1/d);

if (search ~= 0 )
        [Xdenoised,criterion] = denoise_locca_findbest(X,M,K,floor(log(M)),delaymethod,delayspeed,clustermethod,reducemethod,undelaymethod,sortmethod,choosemethod,ratingmethod,gamma,repeat,verbose);
else
    if ( max(length(M(:)),length(K(:))) > 1 )
        [Xdenoised,criterion] = denoise_locca(X,M,K,delaymethod,delayspeed,clustermethod,reducemethod,undelaymethod,sortmethod,choosemethod,ratingmethod,gamma,repeat,verbose);
    else
        Xdenoised = denoise_locca(X,M,K,delaymethod,delayspeed,clustermethod,reducemethod,undelaymethod,sortmethod,choosemethod,ratingmethod,gamma,repeat,verbose);
    end
end

return