[proofplan]
The argument is finite-dimensional linear algebra followed by coefficient comparison. First we use orthogonality to compute the coefficient of $f$ along each basis vector $f_j$ by pairing with $f_j$ under the Petersson inner product. Then uniqueness follows because a basis is linearly independent. Finally, since the expansion of $f$ is a finite linear combination of the $q$-expansions of the basis forms, the $n$-th Fourier coefficient is the same finite linear combination of the $n$-th Fourier coefficients.
[/proofplan]
[step:Project $f$ onto each orthogonal basis vector]
Define complex numbers $c_1,\dots,c_d \in \mathbb{C}$ by
\begin{align*}
c_j := \frac{(f,f_j)_{\mathrm{Pet}}}{(f_j,f_j)_{\mathrm{Pet}}}
\end{align*}
for $j \in \{1,\dots,d\}$. The denominator is nonzero because $f_j$ is a basis vector, hence $f_j \neq 0$, and the Petersson inner product is positive definite on $V$.
Set
\begin{align*}
g := f - \sum_{j=1}^{d} c_j f_j \in V.
\end{align*}
For a fixed index $\ell \in \{1,\dots,d\}$, linearity of the Petersson inner product in the first variable gives
\begin{align*}
(g,f_\ell)_{\mathrm{Pet}}
&= (f,f_\ell)_{\mathrm{Pet}} - \sum_{j=1}^{d} c_j (f_j,f_\ell)_{\mathrm{Pet}}.
\end{align*}
Since the basis is orthogonal, $(f_j,f_\ell)_{\mathrm{Pet}} = 0$ whenever $j \neq \ell$. Therefore
\begin{align*}
(g,f_\ell)_{\mathrm{Pet}}
&= (f,f_\ell)_{\mathrm{Pet}} - c_\ell (f_\ell,f_\ell)_{\mathrm{Pet}} \\
&= (f,f_\ell)_{\mathrm{Pet}} - \frac{(f,f_\ell)_{\mathrm{Pet}}}{(f_\ell,f_\ell)_{\mathrm{Pet}}}(f_\ell,f_\ell)_{\mathrm{Pet}} \\
&= 0.
\end{align*}
Thus $g$ is orthogonal to every vector in the basis $(f_1,\dots,f_d)$.
[guided]
We want to determine the coefficient of $f$ in the direction of a fixed basis vector $f_j$. Orthogonality isolates that coefficient: pairing with $f_j$ kills every other basis vector.
Define
\begin{align*}
c_j := \frac{(f,f_j)_{\mathrm{Pet}}}{(f_j,f_j)_{\mathrm{Pet}}}
\end{align*}
for each $j \in \{1,\dots,d\}$. This is well-defined because $f_j \neq 0$, since $f_j$ is a basis vector, and positive definiteness of the Petersson inner product gives $(f_j,f_j)_{\mathrm{Pet}} > 0$.
Now define the error term
\begin{align*}
g := f - \sum_{j=1}^{d} c_j f_j.
\end{align*}
The purpose of $g$ is to measure what remains after subtracting the proposed expansion. We show that this remainder is orthogonal to every basis vector. Fix $\ell \in \{1,\dots,d\}$. By linearity in the first variable,
\begin{align*}
(g,f_\ell)_{\mathrm{Pet}}
&= (f,f_\ell)_{\mathrm{Pet}} - \sum_{j=1}^{d} c_j (f_j,f_\ell)_{\mathrm{Pet}}.
\end{align*}
Orthogonality of the basis means $(f_j,f_\ell)_{\mathrm{Pet}} = 0$ for $j \neq \ell$, so only the term $j=\ell$ remains:
\begin{align*}
(g,f_\ell)_{\mathrm{Pet}}
&= (f,f_\ell)_{\mathrm{Pet}} - c_\ell (f_\ell,f_\ell)_{\mathrm{Pet}} \\
&= (f,f_\ell)_{\mathrm{Pet}} - \frac{(f,f_\ell)_{\mathrm{Pet}}}{(f_\ell,f_\ell)_{\mathrm{Pet}}}(f_\ell,f_\ell)_{\mathrm{Pet}} \\
&= 0.
\end{align*}
Thus the chosen coefficients make the proposed error term orthogonal to every basis vector.
[/guided]
[/step]
[step:Use the basis property to identify the expansion and prove uniqueness]
Since $(f_1,\dots,f_d)$ is a basis of $V$, there exist complex numbers $b_1,\dots,b_d \in \mathbb{C}$ such that
\begin{align*}
g = \sum_{j=1}^{d} b_j f_j.
\end{align*}
Taking the Petersson inner product with $f_\ell$ for a fixed $\ell \in \{1,\dots,d\}$ gives
\begin{align*}
0
= (g,f_\ell)_{\mathrm{Pet}}
= \sum_{j=1}^{d} b_j (f_j,f_\ell)_{\mathrm{Pet}}
= b_\ell (f_\ell,f_\ell)_{\mathrm{Pet}}.
\end{align*}
Since $(f_\ell,f_\ell)_{\mathrm{Pet}} \neq 0$, we get $b_\ell = 0$. This holds for every $\ell$, so $g=0$. Therefore
\begin{align*}
f = \sum_{j=1}^{d} c_j f_j.
\end{align*}
If also
\begin{align*}
f = \sum_{j=1}^{d} c'_j f_j
\end{align*}
for complex numbers $c'_1,\dots,c'_d \in \mathbb{C}$, then
\begin{align*}
0 = \sum_{j=1}^{d} (c_j-c'_j) f_j.
\end{align*}
[Linear independence](/page/Linear%20Independence) of the basis forces $c_j-c'_j=0$ for every $j$, so $c_j=c'_j$ for every $j$. Hence the expansion is unique.
[/step]
[step:Compare the Fourier coefficients term by term]
For each $n \in \mathbb{N}$, let
\begin{align*}
A_n: V &\to \mathbb{C} \\
h &\mapsto a_n(h)
\end{align*}
denote the map taking a cusp form $h \in V$ to its $n$-th Fourier coefficient at the cusp $\infty$. The Fourier expansion is linear in the form: if $h,h' \in V$ and $\lambda,\mu \in \mathbb{C}$, then the $q$-expansion of $\lambda h+\mu h'$ is obtained by adding the corresponding convergent [Fourier series](/page/Fourier%20Series) term by term, so
\begin{align*}
A_n(\lambda h+\mu h') = \lambda A_n(h)+\mu A_n(h').
\end{align*}
Thus $A_n$ is a complex-[linear map](/page/Linear%20Map).
Applying $A_n$ to the already proved finite expansion
\begin{align*}
f = \sum_{j=1}^{d} c_j f_j
\end{align*}
gives
\begin{align*}
a_n(f)
= A_n(f)
= A_n\left(\sum_{j=1}^{d} c_j f_j\right)
= \sum_{j=1}^{d} c_j A_n(f_j)
= \sum_{j=1}^{d} c_j a_n(f_j).
\end{align*}
Because the sum over $j$ is finite, no limiting operation over the basis index is involved. This proves the asserted Fourier coefficient identity for every $n \in \mathbb{N}$.
[/step]