One of the deepest reasons to use the Fourier expansion is the idea of Harmonic Analysis induced by symmetries. The main topic is called Representation Theory and it forms the foundation to large parts of physics.

Here's the idea: consider a (finite - or compact, for which many similar properties hold) group $G$ and a vector space $V$ on which it operates. Assume for simplicity the field of the vector space has characteristic 0 or at least not 2 (the latter case is the nastiest, with many exceptions). If confused, just assume we have an $\mathbb{R}$- or $\mathbb{C}$-based vector space.

This is to mean that group elements $g\in G$ act as vector space homomorphisms on the vector space $V$. I.e. any group element $g:V\to V$ can be considered a linear map over $V$ whose group multiplication respects homomorphism composition, i.e. $(g_1 \cdot g_2)(v) = (g_1 \circ g_2) (v)$ where $\cdot$ denotes the multiplications as group members while $\circ$ denotes the composition as homomorphisms.

The ensuing structure is called a group representation or $G$-space which has a lot of nice properties. Amongst these is the fact that, for a given group $G$, the $G$-space can be decomposed in only very specific ways into $G$-subspaces. They decompose uniquely into so-called irreducible $G$-subspaces which cannot be further decomposed. Of these irreducible subspaces, there exist only very specific ones for the given group $G$ (for finite groups, there are, in fact, only a finite number of inequivalent irreducible $G$-space types).

It means that, no matter how your $G$-space looks, the decomposition will always contain the same irreducible types, albeit with different contributions - think of this as a kind of "prime-factor decomposition" of your $G$-space; its components will only depend on $G$. You can establish that before computing anything else. This is what molecular dynamics scientists used to establish degeneracy of excitations, purely based on symmetry considerations - before even computing the precise frequencies or knowing the dynamical parameters of the system.

If you have additional properties of your $V$, such as a scalar product, and the group operation respects that, you get additional nice properties, such as that elements from distinct irreducible $G$-subspaces of your $G$-space will be orthogonal (Schur's Lemma). In addition, if your group is abelian, i.e. commutative, and your field is $\mathbb{C}$, the irreducible spaces have dimension 1.

And now we come to the Fourier analysis: the Fourier components are nothing else than decomposition of your vector into its irreducible components when your group $G$ are the translations.

For simplicity, consider the group of translation operators $T_k$ (for technical reasons, assume a finite one-dimensional chain and periodic boundary conditions) which transforms $j \in [n]$ to $T_k(j) := (j+k-1) \mod n + 1$ (i.e. shift by $k$ and wrap around to 1 if past $n$). Given a (complex-valued) vector $v = (v_1, v_2, \dots, v_n)$, a translation operator $T_k$ will transform $v$ as per $T_k v := (v_{T_k^{-1}1}, v_{T_k^{-1}2}, \dots, v_{T_k^{-1}n})$, i.e. all coordinates of $v$ are shifted by $k$ and wrapped around. According to Representation Theory, this $G$-space with the $T_k$ forming the group can be decomposed into one-dimensional irreducible subspaces. It turns out that these are precisely the Fourier components, as is easily verifiable.

Namely: consider our problem's $G$-subspace $F_m := \langle (e^{1 \cdot im}, e^{2\cdot im},\dots,e^{n\cdot im}) \rangle$ (the angles indicate the space generated by the vector inside). The $i$ is the imaginary unit here, and $m$ selects the irreducible subspace (in your usual lingo, the harmonic). Clearly, this is irreducible, as it is one-dimensional. It is also indeed a $G$-space. It is enough to apply $T_k$ to $(e^{1 \cdot im}, e^{2\cdot im},\dots,e^{n\cdot im})$. This gives you $(e^{(1-k)im}, e^{(2-k)im},\dots,e^{(n-k)im})$, but that's the same as $e^{-ikm} \cdot (e^{1 \cdot im}, e^{2\cdot im},\dots,e^{n\cdot im})$, thus just a constant times the original vector and therefore in $F_m$. We have shown that $F_m$ is an irreducible $G$-space for the translations.

Since we have a scalar product, we automatically know from representation theory that the $F_m$ are orthogonal to each other for different $m$. This holds also for other relevant groups. So, for instance if you have the rotation groups on the sphere, you get the spherical harmonics as irreducible spaces. Again, no need to run through tedious integrations on the surface of the sphere - it pops out for all of these, automatically, through Schur's Lemma.

TL;DR: If you have a group which operates on a vector space, with a few additional conditions (group is finite or compact, though some results work for non-compact groups, too; field is nicely behaved, no character 2, for instance), you get a large set of guarantees as to how your vector space can be decomposed and how this decomposition will look like. Additional properties, such as orthogonality etc. will lead to additional guarantees. The Fourier decomposition is the decomposition arising with the special case of the group being the translations of the coordinates and can easily be generalised. The relevant field is called "Representation Theory" and is one of the most elegant topics of mathematics with a strong influence on contemporary physics.