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][1], writes (page [469][2]): "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][3] 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}$. [1]: http://gallica.bnf.fr/ark%3A/12148/bpt6k2970g.f453 [2]: http://gallica.bnf.fr/ark%3A/12148/bpt6k2970g.f469 [3]: http://archive.org/stream/uvresdelagrange01lagr#page/620/