The FFT relies on decomposing the DFT of vectors $\mathbf x\in\mathbb{R}^N$, $N=2^n$ into a sequence of DFTs of lower-dimension. On the other hand, finding a "quantum circuit", at a purely algebraic level, amounts to expressing the procedure operating at the level of individual bits.
I'll try to show how this can be done in a relatively straightforward manner.
FFT
Here's a very rough rundown of the FFT algorithm, focusing on the decomposition it provides for $F_N \mathbf x$.
Suppose we want to compute the DFT of a vector $\mathbf x$ of length $N=2^n$. A standard approach is to decompose $\mathbf x$ into even and odd components, $\mathbf x_e$ and $\mathbf x_o$, respectively, and then express the DFT as
$$F_N \mathbf x = F_N \mathbf x_e+F_N \mathbf x_o
\\= \binom{F_{N/2}\mathbf x_e }{F_{N/2}\mathbf x_e}
+ \binom{\Omega F_{N/2}\mathbf x_o }{-\Omega F_{N/2}\mathbf x_o}
= \binom{F_{N/2} \mathbf x_e + \Omega F_{N/2} \mathbf x_0 }{F_{N/2} \mathbf x_e - \Omega F_{N/2} \mathbf x_0 },\tag 1$$
where $\Omega\equiv \operatorname{diag}(\omega_N^{0},\omega_N^1,...,\omega_N^{N/2})$, and $\omega_N\equiv e^{2\pi i/N}$.
This is the core idea behind the FFT: instead of applying $F_N$ to $\mathbf x$, which has cost $O(N^2)$, we compute $F_{N/2}\mathbf x_e$ and $F_{N/2}\mathbf x_o$, for a total cost of $\sim 2(N/2)^2=N^2/2$ (plus overheads).
This recursive formula is the cornerstone of the FFT algorithm. Iterating this procedure, we manage to compute $F_N\mathbf x$ using $2^k$ applications of $F_{N/2^k}$, and thus using $N$ applications of $F_1$, for a cost of $N$ (times logarithmic overheads).
QFT vs FFT
The first thing to stress when discussing FFT and QFT together is that when talking about DFT/FFT we're usually operating on vectors $\mathbf x\in\mathbb{R}^N$ represented as length-$N$ vectors (as you'd naturally do). However, when $N=2^n$, one can also represent each index for such vectors in binary notation.
Thus, we can express $\mathbf{x}$ as:
$$\mathbf x = \sum_{k=1}^N x_k \mathbf e_k
= \sum_{b_1,...,b_n\in\{0,1\} } x_{b_1,...,b_n} (\mathbf e_{b_1}\otimes\cdots\otimes\mathbf e_{b_n})
\equiv \sum_{\mathbf b} x_{\mathbf b} |\mathbf b \rangle,\tag2$$
where $\mathbf{e}_k$ are the standard basis vectors, and $\vert \mathbf{b} \rangle$ denotes the tensor product of basis vectors corresponding to the binary string $\mathbf{b} = b_1 b_2 \dots b_n$.
This representation highlights that quantum gates operate at the bit level; they manipulate individual bits (qubits) rather than entire vector components. By reformulating the FFT steps to operate on bits, we move closer to constructing a quantum circuit that implements the FFT.
FFT-inspired QFT circuit decomposition
Our objective is now to express the FFT decomposition in terms of (qu)bit operations.
Firstly, note that the even and odd components of the vector $\mathbf{x}$ correspond to indices where the least significant bit (LSB) is 0 and 1, respectively. In our notation, this is the last bit.
This allows us to write:
$$\mathbf x = \mathbf x_e\otimes \mathbf e_0 + \mathbf x_o\otimes \mathbf e_1.\tag3$$
So if we apply the linear operation $F_{N/2}$ "locally" on the first $n-1$ bits, we have
$$(F_{N/2}\otimes I)\mathbf x
= (F_{N/2}\mathbf x_e)\otimes\mathbf e_0
+ (F_{N/2}\mathbf x_o)\otimes\mathbf e_1.\tag4$$
Imagine we now implement a conditional operation: we apply the linear operator $\Omega$ to the first $n-1$ bits, iff the last one is $1$ (in other words, we apply it to the odd indices only). Denote this operation with $C(\Omega)$. This gets us to
$$C(\Omega)(F_{N/2}\otimes I)\mathbf x
= (F_{N/2}\mathbf x_e)\otimes\mathbf e_0
+ (\Omega F_{N/2}\mathbf x_o)\otimes\mathbf e_1.\tag5$$
We're almost there. To get (1) from (5), we need to put together the even and odd terms, with suitable signs. We can do it with a Hadamard operation $H$ acting only on the last bit:
$$\sqrt2
(I_{n-1}\otimes H)C(\Omega)(F_{N/2}\otimes I)\mathbf x
\\= (F_{N/2}\mathbf x_e + \Omega F_{N/2}\mathbf x_o)\otimes\mathbf e_0
+ (F_{N/2}\mathbf x_o - \Omega F_{N/2}\mathbf x_o)\otimes\mathbf e_1.
\tag 6
$$
The remaining difference between equations (1) and (6) is that in equation (1), the components correspond to the first and last $N/2$ indices, whereas in equation (6), they are associated with the even and odd indices.
This discrepancy is trivially resolved by performing a bitwise permutation: rotate bits so that the last one becomes the first one. This leaves us with precisely (1).
You'll notice that the sequence of operations I just outlined is precisely the sequence of gates used to decompose the QFT, with $C(\Omega)$ corresponding for example to the sequence of controlled-phase gates you typically find in these decompositions.
Toy example
To be a bit more explicitly, let's consider how this translates in the case of $N=4$. We have
$$\mathbf x = \begin{pmatrix}x_0\\x_1\\x_2\\x_3\end{pmatrix}
= \begin{pmatrix}x_0 \\ x_2\end{pmatrix}\otimes\binom10
+ \begin{pmatrix}x_1 \\ x_3\end{pmatrix}\otimes\binom01.$$
The DFT decomposes, following (1) as,
$$F_4 \mathbf x
= \begin{pmatrix}F_2 \mathbf x_e\\ F_2\mathbf x_e\end{pmatrix}
+ \begin{pmatrix}\Omega F_2 \mathbf x_o\\ -\Omega F_2\mathbf x_o\end{pmatrix}
= \begin{pmatrix}x_0+x_2 \\ x_0-x_2 \\ x_0+x_2 \\ x_0-x_2\end{pmatrix}
+
\begin{pmatrix}i^0(x_1+x_3) \\ i(x_1-x_3) \\ i^2(x_1+x_3) \\ i^3(x_1-x_3)\end{pmatrix},$$
where $\Omega=S\equiv\operatorname{diag}(1,i)$.
The "QFT-like decomposition" is then
$$F_4 \mathbf x
= (I\otimes H)C(\Omega)(F_2\otimes I)\mathbf x.$$
