undelay.m

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

function sig = undelay(delaysig,d,m,step,method,recon);
%undelay - regenerate the undelayed signal
%
% sig = undelay(delaysig,d,m,step,method,recon);
%
%  recon='avg'  average over possible values 
%  recon='quad' average over (quadratic) weighted possible values 
%
% 01.08.03 - initial version
% 01.11.04 - using undelay architecture
% Peter Gruber

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

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))
  q=m^2+1-[-m:m].^2; q=(2*m+1)*q/sum(q);
  for i=-m:m
   j=i+m+1;
   b=round(i*step);
   switch recon
    case 'avg'
     delaysig(j,:)=circshift(delaysig(j,:)',b)';
    case 'quad'
     delaysig(j,:)=q(j)*circshift(delaysig(j,:)',b)';
   end
  end
  sig=sum(delaysig,1)/size(delaysig,1);
  switch flip
   case 'flip'
    l=length(sig);
    c=ceil(l/2)+1;
    sig=sum([sig(1:c);[sig(1),sig(l:-1:c)]],1)/2;
  end
  return
 end
end

pot=0;
switch recon
 case 'avg'; pot=0;
 case 'quad'; pot=2;
 otherwise; pot=0;
end

sig=undelay_multidim(delaysig,d,m,step,pot,flip);

return


function sig = undelay_multidim(delaysig,d,lag,step,pot,preprocess);
%delay - compute delayed/undelayed signal with lag and step
%
% sig = undelay_multidim(delaysig,d,lag,step,pot,preprocess);
%
%  lag  - distance to consider
%  step - step
%
% 1.10.03 - initial version
% 1.10.04 - rewrite
% Peter Gruber

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

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=[];
weights=[];
while (sum(tcoord)<(length(d)*mdist))
 td=sqrt(sum(tcoord.^2));
 if (td<=lag)
  linoffs=[linoffs tcoord']; 
  weights=[weights (td^pot)]; 
 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);

if ((lo~=size(delaysig,1))||(nd~=size(delaysig,2)))
 error('wrong input (does not match dalayed data)')
end

% normalise weight
weights=weights-(min(weights)-1);
weights=weights/sum(weights);

% shift the dimensions in ns and recontruct sig
sig=zeros(ndims);

for j=1:lo
 tns=reshape(delaysig(j,:),ndims);
 tns=circshift(tns,linoffs(:,j));
 sig=sig+weights(j)*tns;
end

switch preprocess
 case 'flip'
  for i=1:ld
   p=prod(d(1:i-1));
   q=ndims(i);
   s=d(i);
   r=prod(ndims(i+1:ld));
   t=reshape(sig,[p,q,r]);
   sig=(t(:,[1:s],:)+t(:,[1,[q:-1:s]],:))/2;
  end
  sig=reshape(sig,d);
end

return