Fast Hartley Transform Implementation in MATLAB
You can resort to the properties of the DHT with the DFT, given in the wikipedia article. For a real-valued function $f$, its DHT in terms of the DFT is given by

$$
\mathcal{H}f = \Re\{\mathcal{F}f\}-\Im\{\mathcal{F}f\}.
$$

So, a function for fast Hartley transform can be given by

function Hf = fht(f)
F = fft(f);
Hf = real(F) - imag(F);
end


Though, you need to take care of the normalization factors, and tdhis is up to you and how you define the scale factors of the DHT.