A $1$-periodic function, discretized (sampled) at step $h=1/n$, is a collection of $n$ values $y_0,\dots,y_{n-1}$. The functions $f_k = \exp(2\pi i k/n)$, $k=0,\dots,n-1$, form an orthogonal basis among such functions, so there is an expansion $f=\sum c_k f_k$ where
$$
c_k = \frac{1}{n} \sum_{\ell} y_\ell \exp(-2\pi i k\ell/n) \tag{1}
$$
Except for $1/n$, this computation amounts to the multiplication of vector $y$ by the Fourier matrix
$W$ with the entries $\exp(-2\pi i kj)$ (if we index the rows and columns starting with 0). The Discrete Fourier transform is this multiplication (algorithmically, fft arranges it in a clever way, reusing previous intermediate results to speed it up.)
However, one should not jump to the conclusion that
$$
y = \sum_{k=0}^{n-1} c_k \exp(2\pi i k t) \tag{2}
$$
is the right way to interpolate the sampled values.
Although the function on the right does interpolate $y_0,\dots, y_{n-1}$, it has unnatural oscillations in between. One should consider the aliasing of frequencies and fold them appropriately.
For example, take $n=5$, so that (2) says
$$
c_0+ c_1e^{2\pi i t} + c_2e^{4\pi i t} + c_3e^{6\pi i t} + c_4e^{8\pi i t} \tag{3}
$$
The aliasing of frequencies is the fact that $6 \pi t \equiv -4 \pi t \bmod 2\pi$ when $t$ is a grid point (one of $k/5$ points). Similarly, $8\pi t \equiv -2\pi t \bmod 2\pi $. So, we can get smoother curve by folding the frequencies, replacing (3) with
$$
c_0+ c_1e^{2\pi i t} + c_2e^{4\pi i t} + c_3e^{ - 4\pi i t} + c_4e^{-2\pi i t} \tag{4}
$$
The equation (4) is the correct interpolant, and it can be turned into sine-cosine Fourier series (rather, its partial sum) using Euler's formula.
Note that for real functions, the vector of coefficients is conjugate-symmetric: $c_{n-k}=\overline {c_k}$. So one can take the real part of the terms with $c_1, c_2$ and double them; this accounts for $c_3, c_4$. Separate treatment is required for the constant term $c_0$ (which has no counterpart) and, when $n$ is even, for the Nyquist frequency $k=n/2$, which also has no counterpart.
Best Answer
The original motivation for Fourier series was to approximate any periodic function, such as square waves and triangle waves, by sines and cosines. So they naturally come up in acoustics, where the multitude of sounds you hear through your headphones are built up from these sinusoidal functions, which are themselves absurdly simple to generate.
Yet Fourier series are not perfect and their approximation of a square wave is an example. Even after including a lot of sines and cosines, there still remains an overshoot/undershoot at the transition between high and low amplitudes:
This is the Gibbs phenomenon. The study of these overshoots, and how well a Fourier series can approximate any function in general, leads to the most important relation in digital signal processing: the Nyquist–Shannon sampling theorem. This has implications in, for example, capturing audio from a microphone to store in a digital format and then playing it to an audience: how well is that sound?
Fourier series underpin the answers to these questions. While the Fourier transform represents a procedure to convert between the continuous-wave and discrete samples, the series itself remains relevant to the analysis of how well this transformation goes both ways.
Other applications of the actual Fourier series in physics and engineering include the general heat equation (which was the problem Fourier himself solved with the series that got named after him), vibrational modes of structural elements in buildings, quantum harmonic oscillators and generally any place where some function repeats itself over and over.