Home > chronux_1_50 > spikesort > 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 for k = 1:D,
0052     if (k <= N)
0053         if ((numel(H{k})>1) && ~isvectord(H{k})), error('All kernels Hk must be vectors.');  end;
0054         if (~isa(H{k},'double')),  H(k) = double(H{k}(:));  end;  % force col/double
0055         Hisreal(k) = isreal(H{k});
0056     else
0057         H{k} = 1;
0058         Hisreal(k) = 1;
0059     end
0060 end
0061 
0062 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Do the Conv %%%%%%%%%%%%%%%%%%%%%%%%%%
0063 C = A;
0064 for k = 1:N,
0065     orient = ones(1,ndims(A));   orient(k) = numel(H{k});
0066     kernel = reshape(H{k}, orient);
0067 
0068     if (Hisreal(k) & isreal(C))
0069         C = convnc(kernel,C);  
0070     elseif (Hisreal(k) & ~isreal(C))
0071         C = convnc(kernel,real(C)) + j*convnc(kernel,imag(C));
0072     elseif (~Hisreal(k) & isreal(C))
0073         C = convnc(real(kernel),C) + j*convnc(imag(kernel),C);  
0074     else
0075         Hr = real(kernel);    Hi = imag(kernel);
0076         Cr = real(C);         Ci = imag(C);
0077         C = convnc(Hr,Cr) - convnc(Hi,Ci) + j*(convnc(Hi,Cr) + convnc(Hr,Ci)); 
0078     end
0079 end
0080 
0081 %%%%%%%%%%%%%%%%%%%%%%%% Get the right shape %%%%%%%%%%%%%%%%%%%%%%
0082 % nothing more to do for 'full' shape
0083 if (strcmp(shape,'full')),  return;  end;
0084  
0085 % but for 'same' or 'valid' we need to crop the conv result
0086 subs = cell(1,ndims(C));
0087 if (strcmp(shape,'same'))  
0088   for k = 1:D
0089       subs{k}  = floor(length(H{k})/2) + [1:size(A,k)];  % central region
0090   end
0091 elseif (strcmp(shape,'valid'))
0092   for k = 1:D
0093       validLen = max(size(A,k)-length(H{k})+1,0);
0094       subs{k}  = length(H{k})-1 + [1:validLen];
0095   end
0096 end
0097 C = C(subs{:});

Generated on Mon 09-Oct-2006 00:54:52 by m2html © 2003