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.
[Edit: H. Burkhardt in Trigonometrische Reihen und Integrale, Encykl. Math. Wiss. II A 12, p. 929 confirms the above as the first of several papers where Cauchy uses the complex form. Moreover he goes back even further to Laplace who writes in Theorie analytique des probablilités (Paris, 1812), pp. 83-84:
"Take the equation $u=\sum_{x=0}^\infty y_xt^x$. Substitute on both sides $e^{x\varpi\sqrt{-1}}$ for $t^x$... and write $U$ for what $u$ then becomes. Multiplying the equation by $e^{-x\varpi\sqrt{-1}}$ and integrating... the right-hand side boils down to $2\pi y_x$; one has therefore $y_x = \frac1{2\pi}\int U\ d\varpi\ (\cos x\varpi - \sqrt{-1}\sin x\varpi)$".]
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}$.