Given the unitary Fourier matrix \(F\) and its inverse \(F^{\dagger}\) you can construct the two Hartley matrices \(P\) and \(Q\) via:
\begin{align*}
P &\equiv \frac{1}{2}(1 + i)F + \frac{1}{2}(1 - i)F^{\dagger} & Q &\equiv \frac{1}{2}(1 - i)F + \frac{1}{2}(1 + i)F^{\dagger}
\end{align*}
Note that both \(P\) and \(Q\) are hermitian:
\begin{align*}
P^{\dagger} &= P & Q^{\dagger} &= Q
\end{align*}
Both of the Hartley matrices square to the identity:
\begin{align*}
P^{2} &= I & Q^{2} &= I
\end{align*}
The product of the Hartley matrices gives \(S\):
\begin{equation*}
P Q = S
\end{equation*}
Note that \(S^{2} = I\) also. Thus, the four matrices \(I\), \(P\), \(Q\), and \(S\) are a representation of the group \(\mathbb{Z}_{2} \times \mathbb{Z}_{2}\). On the other hand, the four matrices \(F\), \(F^{2} = S\), \(F^{3} = F^{\dagger}\), and \(F^{4} = I\) are a representation of the group \(\mathbb{Z}_{4}\).
The Fourier transform matrix and its inverse can be re-written in terms of the Hartley matrices:
\begin{align*}
F &= \frac{1}{2}(1 - i) P + \frac{1}{2} (1+i) Q \\
F^{\dagger} &= \frac{1}{2}(1 + i) P + \frac{1}{2} (1-i) Q
\end{align*}
In terms of the Hartley matrices, you have involutions. These allow you to use familiar results:
\begin{align*}
\exp{(i a P)} &= I\cos{(a)} + i P \sin{(a)} \\
\exp{(i b Q)} &= I\cos{(b)} + i Q \sin{(b)} \\
\exp{(i c S)} &= I\cos{(c)} + i S \sin{(c)}
\end{align*}
Then
\begin{equation*}
\exp{(i a P)} \exp{(i b Q)} \exp{(i c S)} = I C_{I}(a, b, c) - P C_{P}(a, b, c) - Q C_{Q}(a, b, c) - S C_{S}(a, b, c)
\end{equation*}
with
\begin{align*}
C_{I}(a, b, c) &= \cos{(a)} \cos{(b)} \cos{(c)} - i \sin{(a)} \sin{(b)} \sin{(c)} \\
C_{P}(a, b, c) &= \cos{(a)} \sin{(b)} \cos{(c)} - i \sin{(a)} \cos{(b)} \sin{(c)} \\
C_{Q}(a, b, c) &= \sin{(a)} \cos{(b)} \cos{(c)} - i \cos{(a)} \sin{(b)} \sin{(c)} \\
C_{S}(a, b, c) &= \sin{(a)} \sin{(b)} \cos{(c)} - i \cos{(a)} \cos{(b)} \sin{(c)}
\end{align*}
You need
\begin{align*}
C_{I}(a, b, c) &= 0 & C_{P}(a, b, c) &\neq 0 & C_{Q}(a, b, c) &\neq 0 & C_{S}(a, b, c) &= 0
\end{align*}
Use
\begin{align*}
a &= x + y & b &= x - y
\end{align*}
Then
\begin{align*}
2 C_{I}(x, y, c) &= \left[ \cos{(2y)} + \cos{(2x)} \right] \cos{(c)} - i \left[ \cos{(2y)} - \cos{(2x)} \right] \sin{(c)} \\
2 C_{P}(x, y, c) &= \left[ \sin{(2x)} - \sin{(2y)} \right] \cos{(c)} - i \left[ \sin{(2x)} + \sin{(2y)} \right] \sin{(c)} \\
2 C_{Q}(x, y, c) &= \left[ \sin{(2x)} + \sin{(2y)} \right] \cos{(c)} - i \left[ \sin{(2x)} - \sin{(2y)} \right] \sin{(c)} \\
2 C_{S}(x, y, c) &= \left[ \cos{(2y)} - \cos{(2x)} \right] \cos{(c)} - i \left[ \cos{(2y)} + \cos{(2x)} \right] \sin{(c)}
\end{align*}
The following solutions guarantee \(C_{I} = C_{S} = 0\):
\begin{align*}
x_{k} &= \frac{(4 k + 1)\pi}{4} & x_{l} &= \frac{(4 l - 1)\pi}{4} & y_{m} &= \frac{(4 m + 1)\pi}{4} & y_{n} &= \frac{(4 n - 1)\pi}{4}
\end{align*}
Note that
\begin{align*}
\sin{(2x_{k})} &= 1 & \sin{(2x_{l})} &= -1 & \sin{(2y_{m})} &= 1 & \sin{(2y_{n})} &= -1
\end{align*}
The \((x_{k}, y_{m})\) solution gives
\begin{align*}
C_{P} &= -i \sin{(c)} & C_{Q} &= \cos{(c)}
\end{align*}
Note that
\begin{align*}
a_{km} &= \frac{\pi}{2} + (k + m)\pi & b_{km} &= (k - m) \pi
\end{align*}
The \((x_{k}, y_{n})\) solution gives
\begin{align*}
C_{P} &= \cos{(c)} & C_{Q} &= -i\sin{(c)}
\end{align*}
Note that
\begin{align*}
a_{kn} &= (k + n) \pi & b_{kn} &= \frac{\pi}{2} + (k - n) \pi
\end{align*}
The \((x_{l}, y_{m})\) solution gives
\begin{align*}
C_{P} &= -\cos{(c)} & C_{Q} &= i\sin{(c)}
\end{align*}
Note that
\begin{align*}
a_{lm} &= (l + m) \pi & b_{lm} &= -\frac{\pi}{2} + (l - m) \pi
\end{align*}
The \((x_{l}, y_{n})\) solution gives
\begin{align*}
C_{P} &= i \sin{(c)} & C_{Q} &= -\cos{(c)}
\end{align*}
Note that
\begin{align*}
a_{ln} &= -\frac{\pi}{2} + (l + n) \pi & b_{ln} &= (l - n) \pi
\end{align*}
With \(c = \pm \pi / 4\) and some more work, you can recover \(F\) as an exponential of an anti-hermitian matrix.
