How to reduce for loops for Hartley transform?

Question

I am currently working on 2D Hartley transform. The code is shown below:

for u=1:size(img,1)
    for v=1:size(img,2)
        for x=1:size(img,1)
            for y=1:size(img,2)
                a = 2*pi*u*x/size(img,1);
                b = 2*pi*v*y/size(img,2);
                temp= img(x,y)*(cos(a+b) + sin(a+b));
                X(u,v)= X(u,v)+temp;
            end
        end
    end
end

It has 4 for loops and it takes a very long time to execute this. Is there any method to make it more efficient by reducing the amount of for loops? Anything regarding this would be very helpful.

The formula used for this 2-D Hartley transform is shown below:

Reference:Separable two-dimensional discrete Hartley transform by Andrew B. Watson and Allen Poirson.

Answer 1

If you can fit into memory, you can use bsxfun with some additional singleton dimensions:

N = size(img,1);
M = size(img,2);
x = [1:N].'; %' vector of size N x 1 ( x 1 x 1)
y = 1:M;     %  vector of size 1 x M ( x 1 x 1)
u = permute(1:N,[1 3 2]); %vector of size 1 x 1 x N ( x 1)
v = permute(1:M,[1 3 4 2]); %vector of size 1 x 1 x 1 x M

a = 2*pi/N*bsxfun(@times,u,x); % N x 1 x N x 1
b = 2*pi/M*bsxfun(@times,v,y); % 1 x M x 1 x M
apb = bsxfun(@plus,a,b); % N x M x N x M
%img is N x M (x 1 x 1)

X2 = squeeze(sum(sum(bsxfun(@times,img,cos(apb)+sin(apb)),1),2));

Admittedly this is quite brute-force, one could probably come up with a more memory-efficient solution. The solution heavily utilizes that every array implicitly possesses an infinite number of trailing singleton dimensions, which I tried to note in the comments.

Comparison to your original looping version with N=20; M=30; img=rand(N,M);:

>> max(max(abs(X-X2)))
ans =
     1.023181539494544e-12
>> max(max(abs(X)))
ans =
     3.091143465722029e+02

which means they give the same solution within machine precision.

Answer 2

Here's one using bsxfun and fast matrix-multiplication for a vectorized solution -

%// Store size parameters
[m,n] = size(img);

%// Get vectorized versions of a and b
A = 2*pi*(1:m).'*(1:m)/m;
B = 2*pi*(1:n).'*(1:n)/n;

%// Add vectorized a and b's to form a 4D array. Get cos + sin version.
AB = bsxfun(@plus,permute(A,[1 3 2 4]),permute(B,[3 1 4 2]));
cosAB = cos(AB) + sin(AB);

%// Finally bring in magic of matrix-multiplication for final output
Xout = reshape(img(:).'*reshape(cosAB,m*n,[]),m,n);