function [S] =dht(s,type)
    [N1, N2] = size(s);
    % Discrete Hartley transform
    for i = 0:N1-1
        for k1 = 0:N1-1
            G1(i+1,k1+1) = cos(2*pi*k1*i/N1)+sin(2*pi*k1*i/N1);
        end
    end
    for i = 0:N2-1
        for k1 = 0:N2-1
            G2(i+1,k1+1) = cos(2*pi*k1*i/N2)+sin(2*pi*k1*i/N2);
        end
    end
    s = double(s);
    
    if type == 1
        S = G1*s*G2';
    elseif type  == 2
        S = 1/N1/N2 * G1' * s * G2;
    end
end