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}$.