This a continuation of my first answer. I was trying to edit the previous one but the MathJax processing was freezing my computer. I suppose that answer was getting too long.
@Darij: There has been lots of action on this question since I posted my first answer last night. Although, I must say, there has not been much reading of that answer ... ;)
In your comments you raised two issues.
Issue 1: how does FFT imply part b) of Schur-Weyl duality (the hard part)?
My reply: The proof I gave does in one sweep, with no pause in between, the FFT for $SU(d)$ (and not $U(d)$ as in David's answer) and Schur-Weyl duality. But I could have written it in two steps. In this case the main protagonist $T$ would be an array with entries $$ T_{a_1\ldots a_p,i_1\ldots i_q} $$ corresponding to an invariant of $p$ vectors and $q$ covectors. After eliminating the matrix elements of $g^{-1}$ by Cramer's rule, one has a degree $p+q(d-1)=qd+p-q$ homogeneous polynomial in the entries of $g$ to be hit by a suitable power $r$ of ${\rm det}(\partial g)$. If $d$ does not divide $p-q$ then $T$ must vanish. Otherwise $r=q+\frac{p-q}{d}$ and at the end of the argument above using the elimination of $\epsilon$'s and $\tau$'s in pairs with the identity "${\rm det}(AB)={\rm det}(A){\rm det}(B)$", the number of $\epsilon$'s is $q$ while the number of $\tau$'s is $r$.
- If $p>q$: then $r>q$ and after pairwise elimination of $q$ $\epsilon$-$\tau$ pairs, one is left with $q$ direct contractions of vectors and covectors mediated by a permutation in $S_q$ (this is the Schur-Weyl or $GL$ invariant part) as well as $(p-q)/d$ surviving $\tau$'s giving what the classics called "Klammerfaktoren" involving vectors only.
- If $p<q$: then $r<q$ and after pairwise elimination of $r$ $\epsilon$-$\tau$ pairs, one is left with $p$ direct contractions of vectors and covectors mediated by a permutation in $S_p$ (this is again the Schur-Weyl or $GL$ invariant part) as well as $(q-p)/d$ surviving $\epsilon$'s giving "Klammerfaktoren" involving covectors only.
- What you are interested in is the special case $p=q=n$, or pure Schur-Weyl situation with no "Klammerfaktoren".
Issue 2: This does not count because the proof uses analysis, Haar measure etc.
My reply: The Haar measure is just salad dressing, so the structure of the proof becomes natural. There three steps only: 1) average over the group. 2) figure out what averaging means in algebro/combinatorial terms rather than analytic ones 3) do the calculation and contemplate the result. It would have taken me too long to draw pictures in my answer but this is a purely graphical proof. If you want to see these pictures at least in the case $d=2$, $p$ arbitrary and $q=0$ case, then see pages 16-17 of this article. (If you have access to JKTR, the published version is quite a bit better).
I could have written my first answer in a purely combinatorial way without mentioning the Haar measure at all but then some elements of the proof would have come totally out the blue. In particular, I would have had to all of a sudden say: it $T$ is invariant then $$ T=c_n\ ({\rm det}(\partial g))^n\ gT $$ where $gT$ is the expression obtained after Cramer's rule and the $g$'s are treated as formal variables. So again I insist, my answer is purely combinatorial. But I don't know if it may be useful in nonzero characteristic. This may require understanding the arithmetic of the $c_n$'s.
Additional remark: A more natural proof in the spirit of the one I gave would be to use the FFT for $U(d)$ slash $GL(d)$ rather than $SU(d)$ slash $SL(d)$. This would need an explicit formula for averaging as an infinite order differential operator. I remember seeing something like that in a physics paper, but I would have to find that reference. Usually people resort to Weingarten calculus which goes in the orthogonal direction of the character theory of the symmetric group.
Edit: Yet another rewriting of the proof in order to eliminate analysis and Haar measures etc.
For any $p,q\ge 0$, let $\mathcal{T}_{p,q}$ be the set of arrays $$ T=(T_{a_1,\ldots,a_p,i_1\ldots,i_q})_{a_1,\ldots,a_p,i_1\ldots,i_q\in [d]} $$ of complex numbers. For $g\in GL(d)$ and $T\in\mathcal{T}_{p,q}$ define the new array $gT$ by $$ (gT)_{a_1,\ldots,a_p,i_1\ldots,i_q}=\sum_{b_1,\ldots,b_p,j_1\ldots,j_q\in [d]} T_{j_1,\ldots,j_p,b_1\ldots,b_q} g_{a_1 j_1}\cdots g_{a_p j_p}\ (g^{-1})_{b_1 i_1}\cdots (g^{-1})_{b_q i_q} $$ If $q=0$ then this is well defined for all matrices $g$ and not just invertible ones.
We define two maps $\phi:\mathcal{T}_{p,q}\rightarrow \mathcal{T}_{p+q(d-1),0}$ and $\psi:\mathcal{T}_{p+q(d-1),0} \rightarrow\mathcal{T}_{p,q}$ as follows: $$ (\phi(T))_{a_1\ldots a_p, b_{11}\ldots b_{1 (d-1)},\cdots, b_{q 1}\ldots, b_{q (d-1)}}=\sum_{b_1,\ldots,b_q\in [d]} $$
$$ T_{a_1,\ldots,a_p,b_1\ldots,b_q} \tau_{b_1,b_{11}\ldots b_{1 (d-1)}}\cdots \tau_{b_q, b_{q 1}\ldots b_{q (d-1)}} $$ and $$ (\psi(S))_{a_1,\ldots,a_p,i_1\ldots,i_q}=\sum_{j_{11}\ldots j_{1 (d-1)},\cdots, j_{q 1}\ldots, j_{q (d-1)}\in [d]} $$
$$ S_{a_1\ldots a_p, i_1, j_{11}\ldots j_{1 (d-1)},\cdots, i_q, j_{q 1}\ldots, j_{q (d-1)}} \epsilon_{i_1, j_{11}\ldots j_{1 (d-1)}}\cdots \epsilon_{i_q, j_{q 1}\ldots j_{q (d-1)}} $$ For all invertible $g$ one has $$ g\phi(T)=({\rm det}(g))^{q} \ \phi(gT) $$ and $$ g\psi(S)=({\rm det}(g))^{-q}\ \psi(gS) $$ as trivial consequences of the even more trivial identity $$ \sum_{i_1\ldots i_d\in [d]}\epsilon_{i_1\ldots i_d} M_{i_1 j_1}\cdots M_{i_d j_d}= {\rm det}(M)\epsilon_{j_1\ldots j_d} $$ and also a good pair of glasses.
Moreover, one has $$ \psi(\phi(T))= \frac{1}{(d-1)!^q} \ T $$ as a consequence of Cramer's rule for the identity matrix, namely, $$ \sum_{j_1\ldots j_{d-1}\in [d]}\epsilon_{j,j_1\ldots j_{d-1}}\tau_{i,j_1\ldots j_{d-1}}=(d-1)!\ \delta_{j,i} $$ The above is prep work for a reduction of the FFT for $SU(d)$ to the case $q=0$ and $p$ divisible by $d$.
Let us now assume that $q=0$ and that $T\in\mathcal{T}_{p,0}$ is invariant, i.e., $gT=T$ for all $g\in SU(d)$ or rather $SL(d)$. Since any $g\in GL(d)$ can be written as $g=\lambda h$ with $h\in SL(d)$ and $\lambda$ some $d$-th root of ${\rm det}(g)$ then since $gT$ is a homogeneous polynomial we more generally have $gT=({\rm det}(g))^n\ T$ for all $g\in GL(d)$, and even for all $g\in {\rm Mat}_{d\times d}$. Here $n$ denotes the weight $p/d$. We now have $$ ({\rm det}(\partial g))^n\ gT=({\rm det}(\partial g))^n\ ({\rm det}(g))^n T=\rho_n\ T $$ The heart of the proof (as well as of the understanding of the Haar measure as pure combinatorics), is to show that $$ \rho_n=({\rm det}(\partial g))^n\ ({\rm det}(g))^n \neq 0 $$ The most precise way of doing this is via the "Cayley identity" exercise which gives the explicit formula $$ \rho_n=\frac{1}{c_n} $$ But a cheaper argument goes as follows. When summing over the permutation in $S_{nd}$ which I called a Wick contraction scheme (just the Leibnitz rule) arising in the computation of $\rho_n$, one gets a sum of squares: what I called the two decoupled (rather than connected) components that I mentioned earlier are identical and deliver the same numerical evaluation. One just need choose one of these contraction schemes, for instance, coming from $$ \left[({\rm det}(\partial g))\ ({\rm det}(g))\right]^n $$ and check it is nonzero. The germ of this "mirror symmetry" between the two components is visible in the right-hand side of the basic identity I used $$ \frac{\partial}{\partial g_{cl}}g_{i b}=\delta_{ci}\delta_{lb} $$
I guess that's enough details for you to finish the proof of the FFT for $SU$ or $SL$.
Let me simply say that this is not my proof but Clebsch's in his amazing article "Ueber symbolische Darstellung algebraischer Formen" in Crelle 1861. (Please click on the full text link and read Section 3 of that article).