1D Case
The exponential basis for $[0,1]$ is $b_j(x) = \exp(i2\pi j x)$. The Fourier series of $f$ is:
$$S(f)(x) = \sum_{j=-\infty}^{\infty} \langle b_j, f\rangle b_j(x) = \sum_{j=-\infty}^{\infty} c_j \exp(i2\pi j x), \quad c_j = \langle b_j, f\rangle = \int_0^1 f(x)\exp(-i2\pi j x)\,dx.$$
Suppose $f$ is a real function, then $\overline{f(x)} = f(x)$ and thus
\begin{align} c_{-j} &= \int_0^1 f(x)\exp(-i2\pi (-j)x)\,dx = \int_0^1 \overline{f(x)\exp(-i2\pi jx)}\,dx = \overline{c_{j}}. \end{align}
Then decompose $c_j = \alpha_j + i \beta_j$, where $\alpha_j,\beta_j\in \mathbb{R}$ are the real and imaginary parts of $c_j$. Now you can reduce $S(f)$ to the real setting as follows:
\begin{align} S(f)(x) &= \sum_{j=-\infty}^{\infty} c_j \exp(i2\pi j x) \\ &=c_0 + \sum_{j=1}^{\infty}c_j\exp(-i2\pi j x) + c_{-j}\exp(-i2\pi j x) \\ &=c_0 + \sum_{j=1}^{\infty}c_j\exp(i2\pi j x)+ \overline{c_j\exp(i2\pi j x)} \\ &= c_0 + \sum_{j=1}^{\infty} 2\mathfrak{Re}[c_j\exp(-i2\pi j x)] \\ &= c_0 + 2\sum_{j=1}^{\infty}\mathfrak{Re}[(\alpha_j + i\beta_j)(\cos(2\pi jx) + i\sin(2\pi jx))] \\ &= c_0 + 2\sum_{j=1}^{\infty} \alpha_j\cos(2\pi j x) - \beta_j\sin(2\pi j x). \end{align}
In other words $\alpha_j = \frac{c_j+c_{-j}}{2}$ and $\beta_j = \frac{c_j-c_{-j}}{2i}$, and $c_0 = \int_0^1 f(x)\,dx$. Equivalently you can compute those as: $$\alpha_j+i\beta_j = \int_0^1 f(x) \exp(-i2\pi jx)\,dx = \int_0^1 f(x) \cos(2\pi jx)\,dx - \int_0^1 f(x)\sin(2\pi jx)\,dx.$$ In the literature one often finds $A_j = 2\alpha_j$ and $B_j = -2\beta_j$, and $A_0 = c_0$, then: $$S(f)(x) = A_0 + \sum_{j=1}^{\infty}A_j\cos(2\pi j x) + \sum_{j=1}^{\infty}B_j\sin(2\pi j x) .$$
2D Case
The Fourier series of a 2D function is given as follows:
\begin{align}S(f)(x,y) &= \sum_{j=-\infty}^{\infty}\sum_{k=-\infty}^{\infty} c_{j,k} \exp(i2\pi k y)\exp(i2\pi jx), \\ c_{j,k} &= \int_0^1\left(\int_0^1 f(x,y) \exp(-i2\pi ky)\,dy\right)\exp(-i2\pi jx)\,dx.\end{align}
Again, suppose that $f$ is real, then $\overline{f(x,y)}=f(x,y)$ and thus:
\begin{align}c_{-j,-k} &= \int_{[0,1]^2}f(x,y)\exp(-i2\pi (x,y)\cdot (-j,-k))\,dxdy \\ &= \int_{[0,1]^2}\overline{f(x,y)\exp(-i2\pi (x,y)\cdot (j,k))}\,dxdy = \overline{c_{j,k}}.\end{align}
Once again define $c_{j,k} = \alpha_{j,k} + i\beta_{j,k}$ such that $\alpha_{j,k}, \beta_{j,k}\in\mathbb{R}$ are the real and imaginary parts of $c_{j,k}$. Then proceed as for the 1-D case:
\begin{align} S(f)(x,y) &= \sum_{k=-\infty}^{\infty}\sum_{j=-\infty}^{\infty}c_{j,k}\exp(i2\pi (x,y)\cdot(j,k)) \\ &= c_{0,0} \\ &+ \sum_{j=1}^{\infty} c_{j,0} \exp(i2\pi x j) + c_{-j,0} \exp(-i2\pi x j) \\ &+ \sum_{k=1}^{\infty} c_{0,k} \exp(i2\pi y k) + c_{0,-k} \exp(-i2\pi y k)\\ &+ \sum_{j=1}^{\infty}\sum_{k=1}^{\infty} c_{j,k}\exp(i2\pi (x,y)\cdot(j,k)) + \sum_{j=1}^{\infty}\sum_{k=1}^{\infty} c_{-j,-k}\exp(i2\pi (x,y)\cdot(-j,-k)) \\ &+\sum_{j=1}^{\infty}\sum_{k=1}^{\infty} c_{-j,k}\exp(i2\pi (x,y)\cdot(-j,k)) +\sum_{j=1}^{\infty}\sum_{k=1}^{\infty} c_{j,-k}\exp(i2\pi (x,y)\cdot(j,-k)) \\ &= c_{0,0} \\ &+ \sum_{j=1}^{\infty} c_{j,0} \exp(i2\pi x j) + \overline{c_{j,0} \exp(i2\pi x j)} \\ &+ \sum_{k=1}^{\infty} c_{0,k} \exp(i2\pi y k) + \overline{c_{0,k} \exp(i2\pi y k)}\\ &+ \sum_{j=1}^{\infty}\sum_{k=1}^{\infty} c_{j,k}\exp(i2\pi (x,y)\cdot(j,k)) + \overline{c_{j,k}\exp(i2\pi (x,y)\cdot(j,k))} \\ &+\sum_{j=1}^{\infty}\sum_{k=1}^{\infty} c_{-j,k}\exp(i2\pi (x,y)\cdot(-j,k)) +\overline{c_{-j,k}\exp(i2\pi (x,y)\cdot(-j,k))} \\ &= c_{00} + 2\sum_{j=1}^{\infty} \mathfrak{Re}[c_{j,0} \exp(i2\pi x j)] + 2\sum_{k=1}^{\infty} \mathfrak{Re}[c_{0,k} \exp(i2\pi y k)] \\ &+ 2\sum_{j=1}^{\infty}\sum_{k=1}^{\infty} \mathfrak{Re}[c_{j,k}\exp(i2\pi (x,y)\cdot(j,k))] \\ &+ 2\sum_{j=1}^{\infty}\sum_{k=1}^{\infty} \mathfrak{Re}[c_{-j,k}\exp(i2\pi (x,y)\cdot(-j,k))]\\ &= c_{00} +2\sum_{j=1}^{\infty}\alpha_{j,0}\cos(2\pi jx) - \beta_{j,0} \sin(2\pi j x) + 2\sum_{k=1}^{\infty}\alpha_{0,k}\cos(2\pi ky) - \beta_{0,k} \sin(2\pi k xy)\\ &+ 2\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}\mathfrak{Re}[(\alpha_{j,k}+i\beta_{j,k})(\cos(2\pi jx) + i\sin(2\pi j x))(\cos(2\pi ky) + i\sin(2\pi k y))] \\ &+2\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}\mathfrak{Re}[(\alpha_{-j,k}+i\beta_{-j,k})(\cos(2\pi jx) - i\sin(2\pi j x))(\cos(2\pi ky) + i\sin(2\pi k y))] \\ &= c_{00} +2\sum_{j=1}^{\infty}\alpha_{j,0}\cos(2\pi jx) - \beta_{j,0} \sin(2\pi j x) + 2\sum_{k=1}^{\infty}\alpha_{0,k}\cos(2\pi ky) - \beta_{0,k} \sin(2\pi k xy)\\ &+ 2\sum_{j=1}^{\infty}\sum_{j=1}^{\infty} (\alpha_{j,k}+\alpha_{-j,k})\cos(2\pi j x) \cos(2\pi k y) \\ &- 2\sum_{j=1}^{\infty}\sum_{j=1}^{\infty} (\alpha_{j,k}-\alpha_{-j,k}))\sin(2\pi j x) \sin(2\pi k y) \\ &-2\sum_{j=1}^{\infty}\sum_{j=1}^{\infty} (\beta_{j,k}+\beta_{-j,k})\cos(2\pi j x) \sin(2\pi k y)\\ &-2\sum_{j=1}^{\infty}\sum_{j=1}^{\infty} (\beta_{j,k}-\beta_{-j,k})\sin(2\pi j x) \cos(2\pi k y). \end{align}
You can verify the following identities:
\begin{align} a_{j,k}+a_{-j,k} &= \mathfrak{Re}\left[\int_{[0,1]^2} f(x,y)\exp(-i2\pi xj)\exp(-i2\pi yk)\,dxdy\right] \\ &- \mathfrak{Re}\left[\int_{[0,1]^2} f(x,y)\exp(i2\pi xj)\exp(-i2\pi yk)\,dxdy\right] \\ &=\int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\cos(2\pi yk)-\sin(2\pi xj)\sin(2\pi yk))\,dxdy \\ &+ \int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\cos(2\pi yk)+\sin(2\pi xj)\sin(2\pi yk))\,dxdy \\ &= 2\int_{[0,1]^2} f(x,y)\cos(2\pi xj)\cos(2\pi yk)\,dxdy \\ a_{j,k}-a_{-j,k} &= \int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\cos(2\pi yk)-\sin(2\pi xj)\sin(2\pi yk))\,dxdy \\ &- \int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\cos(2\pi yk)+\sin(2\pi xj)\sin(2\pi yk))\,dxdy \\ &= -2\int_{[0,1]^2} f(x,y)\sin(2\pi xj)\cos(2\pi yk)\,dxdy \\ \beta_{j,k}+\beta_{-j,k} &= -\int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\sin(2\pi yk)+\sin(2\pi xj)\cos(2\pi yk))\,dxdy \\ &- \int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\sin(2\pi yk)-\sin(2\pi xj)\cos(2\pi yk))\,dxdy \\ &= -2\int_{[0,1]^2} f(x,y)\cos(2\pi xj)\sin(2\pi yk)\,dxdy \\ \beta_{j,k}-\beta_{-j,k} &= -\int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\sin(2\pi yk)+\sin(2\pi xj)\cos(2\pi yk))\,dxdy \\ &+ \int_{[0,1]^2}f(x,y)(\cos(2\pi xj)\sin(2\pi yk)-\sin(2\pi xj)\cos(2\pi yk))\,dxdy \\ &= -2\int_{[0,1]^2} f(x,y)\sin(2\pi xj)\cos(2\pi yk)\,dxdy. \end{align}
This reproduces the formula I linked in the comments from this math stack answer. You can see it's rather tedious, on the other hand it gives you the general idea of how things should look like in higher dimensions.
$n$-D Case
If you have an $n$-dimensional problem then the $k$-th type of coefficient $k\in\{0,1\}^n$ at $\vec{j} = (j_1,\ldots,j_n)\in\mathbb{N}_0^n$ is defined as:
$$A_{\vec{k},\vec{j}} = 2^{n-z(\vec{j})}\int_{[0,1]^n} f(\vec{x}) b_{\vec{k},\vec{j}}(\vec{x})\,d\vec{x}, \quad b_{\vec{k},\vec{j}}(\vec{x}) = \prod_{l=1}^{n}(\cos(2\pi j_l x_l))^{1-k_l} (\sin(2\pi j_l x_l))^{k_l}.$$
where $z(\vec{j})$ counts the number of zeroes in the vector $\vec{j}$. Then you can write the Fourier series of $f$ as:
$$S(f) = \sum_{k\in\{0,1\}^n}\sum_{j\in\mathbb{N}^n_0}A_{\vec{k},\vec{j}} b_{\vec{k},\vec{j}}(\vec{x}).$$
Fourier Series over $[p,q]$
If you have a function $\tilde{f}:[p,q]\to\mathbb{R}$ then you can remap it to $f:[0,1]\to\mathbb{R}$ as follows: $$t = \varphi(x) = p + (q-p)x \,\,\text{ and } \,\,f(x) := \tilde{f}(\varphi(x)) = (\tilde{f} \circ \varphi)(x).$$ Since $x = \varphi^{-1}(t) = (t-p)/(q-p)$ you have:
$$\cos(2\pi j t) = \cos(2\pi j (x-p)/(q-p)), \quad \sin(2\pi j t) = \sin(2\pi j (x-p)/(q-p)).$$
Performing the change of variables in the integral yields:
\begin{align} A_{0,j} &= 2^{1-z(j)}\int_0^1 f(x) \cos(2\pi j x)\,dx = \frac{2^{1-z(j)}}{q-p}\int_a^b \tilde{f}(t) \cos(2\pi j (t-p)/(q-p))\,dt,\\ A_{1,j} &= 2^{1-z(j)}\int_0^1 f(x) \cos(2\pi j x)\,dx = \frac{2^{1-z(j)}}{q-p}\int_a^b \tilde{f}(t) \sin(2\pi j (t-p)/(q-p))\,dt, \end{align}
and the series is:
$$S(\tilde{f})(t) = A_{0,0} + \sum_{j=1}^{\infty} A_{0,j} \cos(2\pi j (t-p)/(q-p)) + \sum_{j=1}^{\infty} A_{1,j} \sin(2\pi j (t-p)/(q-p)).$$
The shift of the basis functions is rather unaesthetic, and you can actually get rid of it. The functions $\cos(2\pi jt/\lambda + \alpha)$ and $\sin(2\pi jt/\lambda + \beta)$ would work just as well as $\cos(2\pi jt/\lambda)$ and $\sin(2\pi jt/\lambda)$. Thus if you have a function $\tilde{f}:[p,q]\to\mathbb{R}$ you just need $\lambda=q-p$ and you don't need the shift. This results in:
\begin{align} A_{0,j} &= \frac{2^{1-z(j)}}{\lambda}\int_a^b \tilde{f}(t) \cos(2\pi j t/\lambda)\,dt,\\ A_{1,j} &= \frac{2^{1-z(j)}}{\lambda}\int_a^b \tilde{f}(t) \sin(2\pi j t/\lambda)\,dt,\\ S(\tilde{f})(x) &= A_{0,0} + \sum_{j=1}^{\infty} A_{0,j} \cos(2\pi j t/\lambda) + \sum_{j=1}^{\infty} A_{1,j} \sin(2\pi j t/\lambda). \end{align}
Fourier Series over $[p_1,q_1]\times\ldots[p_n,q_n]$
Let $\mathcal{D} = [p_1,q_1]\times\ldots[p_n,q_n]$. Generalising the previous to $n$ dimensions yields
\begin{align} A_{\vec{k},\vec{j}} &= \frac{2^{1-z(j)}}{\prod_{l=1}^n \lambda_l}\int_{\mathcal{D}} f(\vec{x}) b_{\vec{k},\vec{j}}(\vec{x})\,d\vec{x},\\ S(f)(x) &= \sum_{\vec{k}\in\{0,1\}^n}\sum_{\vec{j}\in\mathbb{N}_0}A_{\vec{k},\vec{j}} b_{\vec{k},\vec{j}}(\vec{x}), \\ b_{\vec{k},\vec{j}}(\vec{x})&= \prod_{l=1}^n (\cos(2\pi j_l x_l/\lambda_l))^{1-k_l}(\sin(2\pi j_l x_l/\lambda_l))^{k_l}. \end{align}
You can see that this answer is a special case of the above.
A Bit of Theory
Motivation
What you have seen for the exponential basis is a special case of an orthogonal basis for a Hilbert space. Below I provide a brief discussion of this.
(Hamel) Basis in Finite Dimensions
In finite-dimensional linear algebra you can write any vector $v$, from a vector space $V$, as a linear combination of some basis $B=(b_1,\ldots,b_n)$ of $V$:
$$v = \sum_{j=1}^n [v]_B^j b_j.$$
There are two main reasons we care about a basis. The first is that we want to be able to find coordinates for any vector $v$ (e.g. in order to store its coordinate representation on a computer). The second is that we want a vector to have a unique coordinate representation, i.e., two different coordinate vectors must not refer to the same vector. The fact that you can write any $v$ as a linear combination of the vectors from $B$ is termed as $B$ spanning $V$, i.e., $V = \operatorname{span}(b_1,\ldots,b_n)$. The fact that the linear combination is unique is equivalent to the basis vector $b_i$ being linearly independent. These two conditions are in fact the definition of a basis. You can relax the span condition, and the you get a $k$-frame, that is a basis for some $k$-dimensional subspace of $V$. You could also relax the linear independence condition, and then you get an overcomplete basis (often also termed a frame).
Orthonormal Basis in Finite Dimensions
An inconvenience of just having a basis $B$ in some vector space $V$ is that the above does not prescribe a procedure for computing the coordinates $[v]_B$. One could of course try to find biorthogonal basis $B^* = (b^1,\ldots,b^n)$ of the dual space $V^*$ such that $b^i(b_j) = \delta^i_j$, but the latter is not always easy to do, or nice to compute.
If you have an inner product $\langle\cdot,\cdot\rangle : V\times V \to \mathbb{F}$, however, and the basis $B$ is orthonormal w.r.t. this inner product, i.e., $\langle b_i, b_j\rangle = \delta_{ij}$, then you can easily compute the coordinates of $v$ in the basis $B$ as $[v]_B^i = \langle b_i, v\rangle$, since $$\langle b_i, v\rangle = \left\langle b_i, \sum_{j=1}^n [v]_B^j b_j\right\rangle = \sum_{j=1}^n [v]_B^j\langle b_i,b_j\rangle = \sum_{j=1}^n [v]_B^j\delta_{ij} = [v]_B^i.$$
Typically the latter is considered when we generalise to infinite dimenions (orthonormality is often relaxed to orthogonality).
Schauder Basis in Separable Hilbert Spaces
Now suppose you have an infinite-dimensional vector space. Then you would again like to have a coordinate decomposition of your vectors, i.e., $v = \sum_{j=1}^{\infty} [v]_B^j b_j := \lim_{n\to\infty} \sum_{j=1}^n[b]_B^j b_j$. For the infinite sum to make sense you typically require that limits remain in the space, i.e., you wish for a complete space. Moreover, if you want an orthonormal basis, then you also need the notion of an inner product. A vector space with an inner product that is complete is known as a Hilbert space. Any separable Hilbert space $\mathcal{H}=(V,\mathbb{F},\langle\cdot,\cdot\rangle)$ has an orthonormal basis $B=(b_1,b_2,\ldots)$ such that for any $v\in V$ you have $v = \sum_{j=1}^{\infty} \langle b_j, v\rangle b_j$. You need the separability in order to have that the basis is countable. Do note that this generalisation of a basis is not referring to a Hamel basis but rather to a Schauder basis.
Complex Exponentials Basis for $L_2([0,1],\mathbb{C})$
Now consider the space of square-integrable "functions" $L_2([0,1],\mathbb{C})$ (functions is in quotes because these are really equivalence classes of functions that are equivalent up to sets of measure zero). Then one possible orthonormal basis for $L_2([0,1],\mathbb{C})$ is given by $b_j(x) = \exp(i2\pi j x)$. The inner product w.r.t. which the basis is orthonormal is the standard inner product for complex-valued functions on $[0,1]$: $$\langle f, g\rangle = \int_{0}^1 \overline{f(x)}g(x)\,dx, \,\,\text{ where }\,\, \overline{z} = \overline{a+bi} = a-bi \,\,\text{ is complex conjugation.}$$
You can verify that $b_j$ are orthonormal as follows: \begin{align}\langle b_k, b_k\rangle &= \int_0^1 \overline{\exp(i2\pi kx)}\exp(i2\pi kx)\,dx = \int_0^1 \exp(-i2\pi kx)\exp(i2\pi kx)\,dx \\ &= \int_0^1 \exp(i2\pi (k-k)x)\,dx = \int_0^1 1\,dx = 1\\ \langle b_k, b_j\rangle &= \int_0^1 \overline{\exp(i2\pi kx)}\exp(i2\pi jx)\,dx = \int_0^1 \exp(-i2\pi kx)\exp(i2\pi jx)\,dx \\ &= \int_0^1 \exp(i2\pi (j-k)x)\,dx = \frac{1}{i2\pi(j-k)}\exp(i2\pi(j-k)x)\bigr|^1_0 \\ &=\frac{\exp(i2\pi(j-k))-1}{i2\pi(j-k)} = \frac{\overbrace{\cos(2\pi(j-k))}^{1}+i\overbrace{\sin(2\pi(j-k))}^{0}-1}{i2\pi(j-k)} = 0. \end{align}
The representation of a "function" $f\in L_2([0,1],\mathbb{C})$ is then:
$$S(f) = \sum_{j=-\infty}^{\infty} \langle b_j, f\rangle b_j = \sum_{j=-\infty}^{\infty} c_j b_j, \quad c_j = \langle b_j, f\rangle = \int_{0}^1 \exp(-2i\pi jx)f(x)\,dx.$$
I wrote "function" again, because $f$ is really an equivalence class of functions. If you had an actual square integrable function $g\in\mathcal{L}_2([0,1],\mathbb{C})$ (as opposed to $g\in L_2([0,1],\mathbb{C})$) then $\|g-S(g)\|_2 = 0$ but not necessarily $g(x) = S(g)(x)$ (i.e. the equality doesn't have to hold pointwise, for when the latter holds see this).
Complex Exponentials Basis for $L_2([p,q], \mathbb{R})$
Consider the complex exponential basis $\exp(i2\pi j x/\lambda)$ where $\lambda = q-p$ and the inner product on $[p,q]$: $\langle p, q\rangle = \int_p^q \overline{f(x)}g(x)\,dx$. It is easy to show that the functions are orthogonal:
\begin{align} \langle b_j, b_k\rangle &= \int_p^q \exp(i2\pi(k-j)x/\lambda)\,dx \\ &= \frac{\exp(i2\pi (k-j)x/\lambda)}{i2\pi(k-j)/\lambda}\biggr|^q_p \\ &= \exp(i2\pi (k-j)p/(q-p))\frac{\exp(i2\pi (k-j)(x-p)/(q-p))}{i2\pi(k-j)/(q-p)}\biggr|^q_p \\ &= \exp(i2\pi (k-j)p/(q-p)) \frac{\overbrace{\exp(i2\pi (k-j))}^{1} - \exp(0)}{i2\pi(k-j)/(q-p)} = 0. \\ \langle b_k, b_k\rangle &=\int_p^q 1\,dx = \lambda. \end{align}
The same orthogonality holds for the cosines and sines. But you have to check products of a cosine with a cosine, a sine with a cosine, and a sine with a sine. Moreover you need to use various trigonometric identities such as $\cos(x)\cos(y) = \frac{1}{2}(\cos(x+y)+\cos(x-y))$. And if $p\ne 0$, the last step where I multiplied by $1 = \exp(i2\pi (k-j)p/(q-p)) \exp(-i2\pi (k-j)p/(q-p))$ is quite a bit more tedious in the real case.