


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.

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{:});