delay.m

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

function delaysig = delay(sig,m,step,method);
%delay - delayed the signal
%
% delaysig = delay(sig,m,step,method);
%
% 01.08.03 - initial version
% 01.11.04 - using new delay architecture 
% Peter Gruber

error(nargchk(2,4,nargin))     % check no. of input args
if (nargin<3) step=1; end
if (nargin<4) method='multidim'; end

d=size(sig);

switch method
 case 'multidim'
  flip='flip';
 otherwise
  flip='no';
end

% fast version for 1-dim signals
if (length(d)==2)
 if ((d(1)==1)||(d(2)==1))
  switch flip
   case 'flip'
    sig=[sig(1:prod(d)) sig(prod(d)-1:-1:2)];
  end
  delaysig=[];
  for i=-m:m
   delaysig=[delaysig;circshift(sig',-round(i*step))'];
  end
  return
 end
end

delaysig=delay_multidim(sig,m,step,flip);

return


function delaysig = delay_multidim(sig,lag,step,preprocess);
%dlay - compute delayed/undelayed signal with lag and step
%
% delaysig = delay_multidim(sig,lag,step,preprocess);
%
%  lag  - distance to consider
%  step - step
%
% 01.10.03 - initial version
% 01.10.04 - rewrite
% Peter Gruber

error(nargchk(2,4,nargin))     % check no. of input args
if (nargin<3) step=1; end
if (nargin<4) preprocess='no'; end

d=size(sig);

switch preprocess
 case 'no'
  ndims=d;
 case 'flip'
  ndims=(2*d)-2;
end

% find offsets and weights
mdist=ceil(lag)+2;
tcoord=ones(1,length(d))*(-mdist);
linoffs=[];
while (sum(tcoord)<(length(d)*mdist))
 td=sqrt(sum(tcoord.^2));
 if (td<=lag) linoffs=[linoffs [tcoord']]; end
 tcoord(1)=tcoord(1)+1;
 for i=1:length(d)-1
  if (tcoord(i)>mdist)
   tcoord(i+1)=tcoord(i+1)+1;
   tcoord(i)=-mdist;
  end
 end
end

% delaydimension
lo=size(linoffs,2);

% trivial case
if (lo==0)
 delaysig=[]; return
end

% length of delayed cordinates
nd=prod(ndims);
ld=length(ndims);

delaysig=zeros(nd,lo);
% generate the delayed permutations
switch preprocess
 case 'flip'
  for i=1:ld
   p=prod(ndims(1:i-1));
   q=d(i);
   r=prod(d(i+1:ld));
   t=reshape(sig,[p,q,r]);
   sig=t(:,[[1:q],[q-1:-1:2]],:);
  end
  sig=reshape(sig,ndims);
  t=circshift(sig,-linoffs(:,1));
 case 'no'
  t=circshift(sig,-linoffs(:,1));
end

% shift all dimensions in ns and generate delaysig
delaysig(:,1)=t(:);
for j=2:lo
 t=circshift(t,linoffs(:,[j-1:j])*[1;-1]);
 delaysig(:,j)=t(:);
end

delaysig=delaysig';

return