Work in progress:
So, my question is: why can't do we use any value (even and odd, spaced by 1 fd), to reduce even more the bandwidth of our signal? Has it relation to interference between symbols in frequency spectrum?
Let's first consider binary or $2$-FSK. We have a carrier signal at fixed frequency $f_c$. Modulation by the binary data changes the frequency of the transmitted signal to $f_c + f_d$ or to $f_c-f_d$. Thus the change in the frequency of the transmitted signal which occurs when the transmitter transmits a $1$ followed by a $0$ is $\pm 2f_d$; the two possible FSK tones at frequencies $f_c + f_d$ and $f_c-f_d$ are spaced $2f_d$$\mathbf 2f_d$ apart in frequency; not $f_d$ apart as the OP apparently believes. (The "$f_d$ apart" model corresponds ta a different system where to send a $0$, a tone at $f_c$ is transmitted whereas to send a $1$, a tone at $f_c+f_d$ is transmitted). The transmitted signal consists of a succession of RF pulses of duration $D$. e assumeLet's assume that $f_c \gg \dfrac 1D$ so that each RF pulse has several periods of the RF sinusoid.
What we would like to arrange is for the two possible tones to be orthogonal over each signaling epoch, that is, the RF pulses
$\cos(2\pi(f_c + f_d)t+\theta_1)$ and $\cos(2\pi(f_c - f_d)t+\theta_0)$ should satisfy
$$\int_{kD}^{(k+1)D}\cos(2\pi(f_c + f_d)t+\theta_1)\cos(2\pi(f_c - f_d)t+\theta_0) \,\mathrm dt = 0$$
for all integers $k$. Why the insistence on orthogonality? Well, ifsuppose that the receiver consists of two matched filters for the two RF pulses and a comparator to determine which filter hasis producing a larger response at the sampling instant. Then, then the output of each matched filter is the cross-correlation function of the input RF pulse with the RF pulse to which the filter is matched, and at the sampling instant, this cross-correlation has value $0$ if the input is the other pulse, and large value (the peak value of the autocorrelation function, in fact) if the input is the one to which the filter is matched (cf. this answer for details.)
Now, assuming without loss of generality that $f_1 > f_0$, we have that
\begin{align}
& ~~~~~~~~\int_{0}^{D}2\cos(2\pi f_1t+\theta_1)\cos(2\pi f_0t+\theta_0) \,\mathrm dt\\
&= \int_{0}^{D}\cos(2\pi (f_1-f_0)t+\theta_1-\theta_0) +\cos(2\pi (f_1+f_0)t+\theta_1+\theta_0) \,\mathrm dt\\
&= \left.\frac{\sin(2\pi (f_1-f_0)t+\theta_1-\theta_0)}{2\pi (f_1-f_0)}
+ \frac{\sin(2\pi (f_1+f_0)t+\theta_1+\theta_0)}{2\pi (f_1+f_0)}\right\vert_{0}^{)D}\\
&= ~~~~~\frac{\sin(2\pi (f_1-f_0)D+\theta_1-\theta_0)-\sin(\theta_1-\theta_0)}{2\pi (f_1-f_0)}\\
&~~~~+\frac{\sin(2\pi (f_1+f_0)D+\theta_1+\theta_0)-\sin(\theta_1+\theta_0)}{2\pi (f_1+f_0)}
\end{align}\begin{align}
& ~~~~~~~~\int_{kD}^{(k+1)D}2\cos(2\pi f_1t+\theta_1)\cos(2\pi f_0t+\theta_0) \,\mathrm dt\\
&= \int_{kD}^{(k+1)D}\cos(2\pi (f_1-f_0)t+\theta_1-\theta_0) +\cos(2\pi (f_1+f_0)t+\theta_1+\theta_0) \,\mathrm dt\\
&= \left.\frac{\sin(2\pi (f_1-f_0)t+\theta_1-\theta_0)}{2\pi (f_1-f_0)}
+ \frac{\sin(2\pi (f_1+f_0)t+\theta_1+\theta_0)}{2\pi (f_1+f_0)}\right\vert_{kD}^{(k+1)D}\\
&= ~~~~~\frac{\sin(2\pi (f_1-f_0)(k+1)D+\theta_1-\theta_0)-\sin(2\pi (f_1-f_0)kD+\theta_1-\theta_0)}{2\pi (f_1-f_0)}\\
&~~~~+\frac{\sin(2\pi (f_1+f_0)(k+1)D+\theta_1+\theta_0)-\sin(2\pi (f_1+f_0)kD+\theta_1+\theta_0)}{2\pi (f_1+f_0)}
\end{align}
The numerators of boththe two fractions in the last RHS above are bounded byboth have values in $2$$[-2,+2]$ while in the denominators we have that $2\pi (f_1+f_0)\gg 2\pi (f_1-f_0)$ and so. So the magnitude of the second fraction is considerably smaller thansmall in comparison to the magnitude of the first. ButNow, bearing in mind that $k$ and $k+1$ are integers, let us note that if $f_1$ and $f_0$ are such that $(f_1-f_0)D$ is an integer, then the first fraction has value $0$ while. Similarly, if $f_1$ and $f_0$ are such that $(f_1+f_0)D$ is an integer, then the second fraction has value $0$. Note that both conditions hold wheneverNow, if $f_0D$$f_1D$ and $f_1D$$f_0D$ both are (distinct) integers ($n$$m$ and $m$$n$, say, with $m> n$) or both are half-integers ($\frac{2m-1}{2}$ and $\frac{2n-1}{2}$, say), then both $(f_1-f_0)D$ and give us$(f_1+f_0)D$ are integers, and the desired orthogonalitytwo possible RF pulses are orthogonal over the interval $[kD, (k+1)D]$. We will restrict attention toNote that in all symbol epochs $[kD,(k+1)D]$, it is either the case that both RF pulses have an integer number of periods of the sinusoid, or it is the case onlythat both RF pulses have a half-integer number of periods of the sinusoid.
Applying the above notion to
$$\int_{0}^{D}\cos(2\pi(f_c + f_d)t+\theta_1)\cos(2\pi(f_c - f_d)t+\theta_0) \,\mathrm dt = 0$$$$\int_{kD}^{(k+1)D}\cos(2\pi(f_c + f_d)t+\theta_1)\cos(2\pi(f_c - f_d)t+\theta_0) \,\mathrm dt = 0$$
we see that it must be that both $(f_c+f_d)D$ and $(f_c-f_d)D$ are integers, or they both are half-integers.