  The set $X$ of $k$-subsets of an $n$-set has a natural graph structure where
two $k$-subsets are connected if and only if the size of their intersection is
$k-1$. This is known as the Johnson graph. The symmetric group $S_n$ acts on
the space of complex functions on $X$ and this space has a multiplicity-free
decomposition as sum of irreducible representations of $S_n$, so it has a
well-defined Gelfand-Tsetlin basis up to scalars. The Fourier transform on the
Johnson graph is defined as the change of basis matrix from the delta function
basis to the Gelfand-Tsetlin basis.
  The direct application of this matrix to a generic vector requires
$\binom{n}{k}^2$ arithmetic operations. We show that this matrix can be
factorized as a product of $n-1$ orthogonal matrices, each one with at most two
nonzero elements in each column. The factorization is based on the construction
of $n-1$ intermediate bases which are parametrized via the Robinson-Schensted
insertion algorithm. This factorization shows that the number of arithmetic
operations required to apply this matrix to a generic vector is bounded above
by $2(n-1) \binom{n}{k}$.
  We show that each one of these sparse matrices can be constructed using
$O(\binom{n}{k})$ arithmetic operations. Our construction does not depend on
numerical methods. Instead, they are obtained by solving small linear systems
with integer coefficients derived from the Jucys-Murphy operators. Then both
the construction and the succesive application of all these $n-1$ matrices can
be performed using $O(n \binom{n}{k})$ operations.
  As a consequence, we show that the problem of computing all the weights of
the isotypic components of a given function can be solved in $O(n
\binom{n}{k})$ operations, improving the previous bound $O(k^2 \binom{n}{k})$
when $k$ asymptotically dominates $\sqrt{n}$.
