**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$ apart in frequency; _not_ $f_d$ apart as the OP apparently believes. The transmitted signal consists of a succession of RF pulses of duration $D$. e 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, if the receiver consists of two _matched filters_ for the two RF pulses and a comparator to determine which filter has larger response at the sampling instant, 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](https://dsp.stackexchange.com/a/9389/235) 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}
The numerators of both fractions in the last RHS above are bounded by $2$ while $2\pi (f_1+f_0)\gg 2\pi (f_1-f_0)$ and so the second fraction is considerably smaller than the first. But, 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 if $(f_1+f_0)D$ is an integer, the second fraction has value $0$. Note that both conditions hold whenever $f_0D$ and $f_1D$ both are integers ($n$ and $m$, say) or both are _half-integers_ ($\frac{2m-1}{2}$ and $\frac{2n-1}{2}$, say), and give us the desired orthogonality. We will restrict attention to the integer case only.

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$$
we see that