Here is a sketch. Feel free to ask for clarifications if my writing gets too terse or confusing in places :-).
First recall the Bernstein (a.k.a. Hoeffding, Chernov, etc.) bound. If $Y_m$ are mean $0$
independent random variables bounded by $s$, then for $Y=\sum_{m=1}^n Y_m$, we have for every
positive $t$,
$$
P(|Y|\ge t)\le Ce^{-c\frac{t^2}{ns^2}}\,.
$$
with some $C,c>0$. The proof goes via the consideration of $Ee^{\beta Y}$ with appropriately
chosen $\beta$, as usual.
We want a small refinement of this bound. Suppose that we know in addition that each $Y_m$ is
non-zero with probability at most $p\in(0,1)$. Then, conditioning upon the events that some
$k$ of $Y_m$ have any chance to be non-zero, we get the bound
$$
P(|Y|\ge t)\le C\sum_{k=1}^n e^{-c\frac{t^2}{ks^2}}{n\choose k}p^k(1-p)^{n-k}\,.
$$
Using the inequality $\frac 1k\ge 2\beta-\beta^2 k$, we can estimate the RHS by
$$
Ce^{-2c\frac{t^2}{s^2}\beta}\sum_{k=1}^n {n\choose k}e^{c\beta^2\frac{t^2}{s^2}k}p^k(1-p)^{n-k}
=Ce^{-2c\frac{t^2}{s^2}\beta}\left[1+p(e^{c\beta^2\frac{t^2}{s^2}}-1)\right]^n\,.
$$
for any $\beta>0$ we want.
Now let us look at the distribution of our random variable $X$ assuming that $E|X|^q=1$. For
every $p$, we can take the set $F$ of probability $p$ on which it attains the largest values
and split $X$ as $X'+X''$ where $X'=X$ outside $F$ and $X'=\frac 1p E(X\chi_F)$ on $F$. If
$EX=0$, then $EX'=EX''=0$, $E|X'|^q\le E|X|^q$, $E|X''|^q\le CE|X|^q$ but $X''$ is not zero
only with probability $\le p$ and $|X'|\le s$ where $s^qp=1$ (by Jensen). We can apply this
trick successively with $p=2^{-r}, \log_2 n\ge r\ge 0$ and get the decomposition of $X$ into
the sum of mean $0$ random variables $Z+\sum_{r=0}^{\log_2n}{X_r}$ where $Z$ is different from
$0$ with probability about $\frac 1n$ and $E|Z|^q\le C$, while for each $r$, we have $|X_r|\le
s_r=2^{r/q}$ and $X_r$ is not zero only with probability $p_r=2^{-r+1}$.
Now it will be enough to treat each $X_r$ and $Y$ separately. The exact nature of the discrete
Fourier transform does not matter. All we need to know is that we are interested in the
maximum of $n$ linear forms of $n$ iid copies of our random variables with coefficients not
exceeding $\frac 1{\sqrt n}$.
Let's start with $n$ iid copies of $Z$. Let $N$ be the number of non-zero values among $Z_1,
\dots, Z_n$. Notice that it is the sum of $n$ iid Bernoulli random variables each of which is
$1$ with probability $\frac 1n$, so $Ee^N=(1+\frac en)^n\le e^e$. In particular, any fixed
moment of $N$ is bounded by some constant.
We have for $q>2$,
$$
\frac 1n \left(\sum_m|Z_m|\right)^2\le \frac 1n
N^{2-\frac 2q}\left[\sum_m |Z_m|^q\right]^{2/q}
$$
so, using this crude bound and the Holder inequality, we get
$$
E(\max(Z-\text{forms})^2)\le\frac 1n (EN^{\text{something}})^{1-\frac 2q}\left[E\sum_m |Z_m|
^q\right]^{2/q}\le Cn^{\frac 2q-1}\to 0
$$
as $n\to\infty$.
Thus this part is negligible for large $n$.
Now let us fix $r$ and consider $X_r$. Notice that we can find $\delta>0$ such that $\frac 1q
+\delta<\frac 12$. Then if we change the notation $X_r$ to $2^{-\delta r}X$, $p_r=2^{-r+1}$ to
$p$, and $s_r=2^{r/q}$ to $2^{-\delta r}s$, we shall still have $s\le n^{\frac 12-
\varepsilon}$ with some $\varepsilon>0$ and $p\le s^{-2}$ for all $r\le\log_2n$. The extra
exponential factor $2^{-\delta r}$ we introduced is strong enough to enable us to consider
each such $X$ separately and just to get a uniform bound of order $\log n$ for the
expectations of the squared maximum of $n$ linear $X$-forms $L_j$.
Now comes the trick: in order to show that $E(\max_{1\le j\le n}|L_j|)^2\le C^2\log n+O(1)$,
it suffices to show that for each individual $j$, we have
$$
E(|L_j|^2-C^2\log n)_+=\int_{C\sqrt{\log n}}^\infty 2tP(|L_j|>t)\,dt\le \frac Cn\,.
$$
However, we have the refined bound for the probability in question and, taking into account
the $\sqrt n$ in the denominator, changing $\beta$ to $\beta/n$, and using the bound $p\le s^
{-2}$, we can rewrite it as
$$
P(|L_j|>t)\le Ce^{-2c\frac{t^2}{s^2}\beta}\left[1+s^{-2}(e^{c\beta^2\frac{t^2}{ns^2}}-
1)\right]^n
$$
We have to choose the optimal $\beta=\beta(t)$. There are two cases to consider:
Case 1: $t^2s^2\le n$.
In this case we can take $\beta$ a small multiple of $s^2$ (the difference of the exponent and $1$ can be treated like a linear function in this range) and get the estimate $e^{-c't^2}$, which is as good as if we were dealing with (sub)Gaussian variables.
Case 2: $t^2s^2>n$. In this case we still want to stay in the linear range for the exponent in the parentheses, so we are forced to take $\beta$ a small multiple of $\frac{s\sqrt n}t$. Fortunately, it still gives the bound $e^{-c'\frac{\sqrt n}{s}t}$ and even if we integrate it against $t$ from $1$, we still obtain something like $e^{-c'\frac{\sqrt n}s}\le e^{-c'n^\varepsilon}$, which is much smaller than what we need for large $n$.
