Home > chronux_2_00 > spikesort > utility > datatools > convnsep.m

convnsep

PURPOSE ^

CONVNSEP N-dimensional convolution with separable kernels.

SYNOPSIS ^

function C = convnsep(varargin)

DESCRIPTION ^

CONVNSEP          N-dimensional convolution with separable kernels.
   C = CONVNSEP(H1,H2,...,A) performs an N-dimensional convolution of A
   with the separable kernel given by the outer product of the H1,H2,...
   such that the input vector Hk is convolved along the kth dimension of
   the real N-D array A.  If the number of kernels Hk is equal to M-1,
   then the (HM..HN) kernels are taken to be 1.  No convolution occurs
   for dimensions k with corresponding kernel Hk = 1.

   C = CONVN(H1,H2,...,A,'shape') controls the size of the answer C:
     'full'   - (default) returns the full N-D convolution
     'same'   - returns the central part of the convolution that
                is the same size as A.
     'valid'  - returns only the part of the result that can be
                computed without assuming zero-padded arrays.  The
                size of the result is max(size(A,k)-size(Hk,k)+1,0)
                in the kth dimension.

   See also CONVN, CONV2.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function C = convnsep(varargin)
0002 %CONVNSEP          N-dimensional convolution with separable kernels.
0003 %   C = CONVNSEP(H1,H2,...,A) performs an N-dimensional convolution of A
0004 %   with the separable kernel given by the outer product of the H1,H2,...
0005 %   such that the input vector Hk is convolved along the kth dimension of
0006 %   the real N-D array A.  If the number of kernels Hk is equal to M-1,
0007 %   then the (HM..HN) kernels are taken to be 1.  No convolution occurs
0008 %   for dimensions k with corresponding kernel Hk = 1.
0009 %
0010 %   C = CONVN(H1,H2,...,A,'shape') controls the size of the answer C:
0011 %     'full'   - (default) returns the full N-D convolution
0012 %     'same'   - returns the central part of the convolution that
0013 %                is the same size as A.
0014 %     'valid'  - returns only the part of the result that can be
0015 %                computed without assuming zero-padded arrays.  The
0016 %                size of the result is max(size(A,k)-size(Hk,k)+1,0)
0017 %                in the kth dimension.
0018 %
0019 %   See also CONVN, CONV2.
0020 
0021 % Modified from TMW's CONVN for efficiency.
0022 %  The following commented code tests CONVNSEP against CONVN:
0023 % D = 50;  R = 3;  sd = 1;
0024 % data = randn([D,D,D]) + i*randn([D,D,D]);
0025 % gauss1 = gausskernel(R,sd);    gauss3 = gausskernel([R,R,R],sd);
0026 % tic; orig = convn(data,gauss3,'same');  toc
0027 % tic; modf = convnsep(gauss1,gauss1,gauss1,data,'same'); toc
0028 % mse = sqrt(mean(abs(orig(:)-modf(:)).^2))   % mean squared error
0029 
0030 
0031 %%%%%%%%%%%%%%%%%%%%%%%%%%% Parse Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%
0032 % determine output shape
0033 if (nargin < 2), error('At least two arguments are required.'); end;
0034 if (ischar(varargin{end})),  
0035     shape = varargin{end};   varargin = varargin(1:end-1);
0036     if (length(varargin) == 1), 
0037         C = varargin{1};   return;     % If no kernels specified, no work to do
0038     end;
0039 else
0040     shape = 'full';
0041 end
0042 
0043 % get target matrix
0044 A = varargin{end};   varargin = varargin(1:end-1);
0045 if (~isa(A,'double')),   A = double(A);  end;
0046 D = ndims(A);
0047 
0048 % get kernels
0049 H = varargin;   N = length(H);
0050 if (N > D),  error('Can not have more kernels than the number of dimensions in A.'); end;
0051 Hisreal = ones(D,1);
0052 for k = 1:D,
0053     if (k <= N)
0054         if ((numel(H{k})>1) && ~isvectord(H{k})), error('All kernels Hk must be vectors.');  end;
0055         if (~isa(H{k},'double')),  H{k} = double(H{k}(:));  end;  % force col/double
0056         Hisreal(k) = isreal(H{k});
0057     else
0058         H{k} = 1;
0059         Hisreal(k) = 1;
0060     end
0061 end
0062 
0063 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Do the Conv %%%%%%%%%%%%%%%%%%%%%%%%%%
0064 C = A;
0065 for k = 1:N,
0066     orient = ones(1,ndims(A));   orient(k) = numel(H{k});
0067     kernel = reshape(H{k}, orient);
0068 
0069     if (Hisreal(k) && isreal(C))
0070         C = convnc(kernel,C);  
0071     elseif (Hisreal(k) && ~isreal(C))
0072         C = convnc(kernel,real(C)) + j*convnc(kernel,imag(C));
0073     elseif (~Hisreal(k) && isreal(C))
0074         C = convnc(real(kernel),C) + j*convnc(imag(kernel),C);  
0075     else
0076         Hr = real(kernel);    Hi = imag(kernel);
0077         Cr = real(C);         Ci = imag(C);
0078         C = convnc(Hr,Cr) - convnc(Hi,Ci) + j*(convnc(Hi,Cr) + convnc(Hr,Ci)); 
0079     end
0080 end
0081 
0082 %%%%%%%%%%%%%%%%%%%%%%%% Get the right shape %%%%%%%%%%%%%%%%%%%%%%
0083 % nothing more to do for 'full' shape
0084 if (strcmp(shape,'full')),  return;  end;
0085  
0086 % but for 'same' or 'valid' we need to crop the conv result
0087 subs = cell(1,ndims(C));
0088 if (strcmp(shape,'same'))  
0089   for k = 1:D
0090       subs{k}  = floor(length(H{k})/2) + (1:size(A,k));  % central region
0091   end
0092 elseif (strcmp(shape,'valid'))
0093   for k = 1:D
0094       validLen = max(size(A,k)-length(H{k})+1,0);
0095       subs{k}  = length(H{k})-1 + (1:validLen);
0096   end
0097 end
0098 C = C(subs{:});

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003