mul_circ_covar.m

% Computes circulant covariance of a and then multiplies x by covariance.
% Both images are of size (m1, m2) with k channels.
%
% size(a) is [m1, m2, k]
% size(x) is [m1, m2, k]
function b = mul_circ_covar(a, x)
    [m1, m2, k] = size(x);

    % compute circulant covariance from shifts of a
    % sf(u1,u2,p,q) = conj(af(u1,u2,q)) * af(u1,u2,p)
    af = fft2(a);
    sf = bsxfun(@times, permute(conj(af), [1, 2, 4, 3]), af);

    % multiply x by covariance
    % bf(u1,u2,p) = sum_q sf(u1,u2,p,q) * xf(u1,u2,q)
    xf = fft2(x);
    bf = sum(bsxfun(@times, sf, permute(xf, [1, 2, 4, 3])), 4);
    b = ifft2(bf, 'symmetric');
end