I believe the answer to the question in the title is probably Cauchy, who in Méthode simple et générale pour la détermination numérique des coefficients que renferme le développement de la fonction perturbatrice, C. R. Acad. Sci. Paris 11 (1840), 453-475, writes (page 469):
"La formule $$ Q = \sum Q_{h,h'} e^{h(p-\varpi)\sqrt{-1}} e^{h'(p'-\varpi')\sqrt{-1}} $$ entraîne l'équation $$ Q_{h,h'} = \frac1{4\pi^2} \int_0^{2\pi} \int_0^{2\pi} Q e^{-h(p-\varpi)\sqrt{-1}} e^{-h'(p'-\varpi')\sqrt{-1}} dp\\,dp'." $$ Specializing to $\varpi=0$ one gets the requested complex Fourier series and formula for its coefficients.
Note that Cauchy does it here in two variables -- maybe someone with easy access to his collected works can find an earlier occurrence in one variable.
Note also that specific complex Fourier series were written much earlier. For instance Lagrange in 1766 computes $$ (1-\alpha e^{-i\theta})^{-s} = \sum_{m=0}^{\infty}\binom{-s}{m}(-\alpha)^me^{im\theta} $$ ... except that he still writes everywhere $e^{im\theta}$ in the form $\cos m\theta +\sin m\theta\sqrt{-1}$.