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

 1. 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.
 2. 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.
 3. 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 are 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 <a href="https://arxiv.org/abs/0904.1734">this article</a>. (If you have access to JKTR, the <a href="http://www.worldscientific.com/doi/abs/10.1142/S0218216511009522">published version</a> 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
<a href="https://eudml.org/doc/147825">"Ueber symbolische Darstellung algebraischer Formen"</a> in Crelle 1861.
(Please click on the full text link and read Section 3 of that article and in particular pages 12 and 13 which contain the sum of squares non-vanishing argument).