How to implement fast algorithm for convolution integrals with space and time variant kernels?

I am trying to implement an algorithm that is described in this paper , which is for convolution integrals with space and time variant kernels. That is, if we have an integral of this form:

,where the $\sigma$ is a function of space, we cannot simply use the multiplication in Fourier space, and so we can approximate the kernel with a series of kernels with different widths, in this way:

So for a 2D field (matrix) $\sigma$, we calculate the approximation coefficient 2D fields $\alpha_l$ for a given set of fixed width of the kernels $\sigma_l$, using this formula:

The code I am using in `matlab` for now is the following:

    function alpha=calc_alpha(sigma,sigma_l)
    alpha = zeros(size(sigma,1),size(sigma,2),size(sigma_l,2));
    for i = (1:size(sigma_l,2))
        pre_alpha=(2*sigma.^2)./(sigma.^2 + sigma_l(i)^2);
        prod = ones(size(sigma));
        for j = (1:size(sigma_l,2))
            if j~=i
                prodA=(sigma.^2 - sigma_l(j)^2)/(sigma_l(i)^2 - sigma_l(j)^2);
                prodB=(sigma_l(i)^2 + sigma_l(j)^2)./(sigma.^2 + sigma_l(j)^2);
                prod = prod*prodA*prodB;
            end
        end
        disp(prod);
        alpha(:,:,i)=pre_alpha*prod;
    end
    end

Is a way to use broadcasting in order to speed up this kind of calculation by avoiding the for loops?