I thank everyone for their help.
I've studied the impact of non-linear transformations to harmonic and bi-harmonic signals.
I'm not sure if i use the proper terms because i've got them from the books in other language. So you can fix me in the comment and i'll fix them here.
So.
Any non-linear transformation may be approximated by a polynomial therefore we can put our signal to this one and checks the output harmonics. For example if we use the polynomial of 3 we can calculate the harmonics:
$$ i = a_0 + a_1(u - U_0) + a_2(u - U_0)^2 + a_3(u - U_0)^3 $$
Where $U_0$ - the offset voltage and $u$ - the input signal.
Let's firstly check how this non-linear transformation affects to a harmonic signal.
Using input signal $u(t) = U_0 + U_mcos(\omega t)$ we get:
$$ i(t) = a_0 + a_1U_mcos(\omega t) + a_2U_m^2cos^2(\omega t) + a_3U_m^3cos^3(\omega t) = a_0 + a_1U_mcos(\omega t) + a_2U_m^2(\frac{1}{2}+\frac{1}{2}cos(2\omega t)) + a_3U_m^3(\frac{3}{4}cos(\omega t) + \frac{1}{4}cos(3\omega t)) = (a_0 + \frac{1}{2}a_2U_m^2)+(a_1U_m+\frac{3}{4}a_3U_m^3)cos(\omega t)+\frac{1}{2}a_2U_m^2cos(2\omega t) + \frac{1}{4}a_3U_m^3cos(3\omega t)=I_0 + I_{m1}cos(\omega t)+ I_{m2}cos(2\omega t)+ I_{m3}cos(3\omega t) $$
, where $I_0$ - DC offset, $I_{m1} - I_{m3}$ - harmonics.
As we can see there is 3 harmonics in the output signal. By adding degrees to the polynomial, we get additional harmonics that will decrease infinitely in amplitude. For example let's get the diode approximation at -1 to 1 with the following coefficients:
$$ \frac{1}{3} x^3 + \frac{1}{2} x^2 + \frac{1}{6} x $$
$\hskip2in$
Substituting the coefficients in the formula and use $U_m = 1$ we get:
$I_0 = 0.25$, $I_1 = \frac{5}{12} \approx 0.416$, $I_2 = 0.25$, $I_3 = \frac{1}{12} \approx 0.083$.
Okay let's test the result experimentally using clipping and FFT:
$\hskip1in$
Okay, that's good. If we add the highest degrees we can approximate it more accurate. The other approximation we can use here is piecewise-linear one. This allows to calculate the harmonics of a cut-offed harmonic signal based on a cut-off angle. The animated graph of the amplitudes of the harmonics versus cut-off angle:
This functions is called Berg's functions. The smaller cut-off angle (the impulse is narrowed), the smaller harmonics attenuation.
This is was for a harmonic signal. The other behavior we can see if we use a biharmonic signal (or multiharmonic one). A signal which consists the sum of two harmonic oscillation with different frequencies $\omega_1$, $\omega_2$ and amplitudes $U_{m1}$, $U_{m2}$ is called biharmonic:
$$U_0 + U_{m1}cos(\omega_1t)+ U_{m2}cos(\omega_2t)$$
Substituting this signal to an approximation polynomial we can calculate the harmonics of this signal. To simplify the calculation we use the polynomial of degree 2:
$$ i = a_0 + a_1(u - U_0) + a_2(u - U_0)^2 \Rightarrow $$
$$ i(t) = a_0 + a_1U_{m1}cos(\omega_1t) + a_2U_{m2}cos(\omega_2t) + a_2U_{m1}^2cos^2(\omega_1t) + 2a_2U_{m1}U_{m2}cos(\omega_1t)cos(\omega_2t)+ a_2U_{m2}^2cos^2(\omega_2t) $$
Using following trigonometric functions:
$$cos^2(\psi)=\frac{1}{2}(1 + cos(2\psi))$$
$$cos(\psi_1)cos(\psi_2)=\frac{1}{2}(cos(\psi_1 + \psi_2)+cos(\psi_1 - \psi_2)) \Rightarrow$$
$$(a_0 + \frac{a_2}{2}(U_{m1}^2 + U_{m2}^2)) + a_1U_{m1}cos(\omega_1t) + a_1U_{m2}cos(\omega_2t)+\frac{a_2U_{m1}^2}{2}cos(2\omega_1t)+\frac{a_2U_{m2}^2}{2}cos(2\omega_2t) + a_2U_{m1}U_{m2}cos((\omega_1+\omega_2)t)+a_2U_{m1}U_{m2}cos((\omega_1-\omega_2)t)$$
As we can see there is the new harmonics which weren't there when we processed a harmonic signal with frequencies $\omega_1\pm\omega_2$. Those frequencies are called combined frequencies. Let's check the spectrum using FFT experimentally. $U_{m1} = 0.5; U_{m2} = 0.5$ Using the following polynomial:
$$ \frac{1}{2} x^2 + \frac{1}{2} x \Rightarrow$$
$$I_0 = 0.125; I_{\omega_1} = 0.25; I_{\omega_2} = 0.25; I_{2\omega_1} = 0.0625; I_{2\omega_2} = 0.0625; I_{\omega_1+\omega_2} = 0.125; I_{\omega_1-\omega_2} = 0.125$$
$\hskip1in$
Using polynomials with the higher degree $N$ the result contains the combined frequencies $p\omega_1 \pm q\omega_2$, where $p + q = N$; $p, q\in \mathbb{N}$.
For example the polynomial of 3 degree contains the following frequencies:
┌─────┬─────┬────────┬────────┬─────┐
│ q\p │ 0 │ 1 │ 2 │ 3 │
├─────┼─────┼────────┼────────┼─────┤
│ 0 │ I0 │ ω1 │ 2ω1 │ 3ω1 │
│ 1 │ ω2 │ ω1±ω2 │ 2ω1±ω2 │ - │
│ 2 │ 2ω2 │ ω1±2ω2 │ - │ - │
│ 3 │ 3ω2 │ - │ - │ - │
└─────┴─────┴────────┴────────┴─────┘
The spectrum and the waveform are more accurate:
$\hskip1in$
When we have the bunch of the input harmonics we can represents them as:
$$u(t) = U_0 + \sum_{k = 1}^{\infty}U_{mk}cos(\omega_kt -\varphi_k)$$
Depending on the degree $N$ of the approximation polynomial we have the following combined frequencies:
$$p\omega_1\pm q\omega_2\pm s\omega_3\pm ... \pm k\omega_k\pm ...$$
$$p + q + s + ... + k + .. = N$$
$$p, q, s, k \in \mathbb{N}$$
For example if we use the polynomial of degree 2 with a bi-harmonic input the result spectrum contains DC, two first harmonics of each input frequencies and combined frequencies $\omega_1 \pm \omega_2, \omega_1 \pm \omega_3, \omega_2 \pm \omega_3$. Using polynomial of degree 3 third harmonics appears $3\omega_1, 3\omega_2, 3\omega_3$ with the combined frequencies $\omega_1 \pm \omega_2 \pm \omega_3; 2\omega_1 \pm \omega_3; \omega_1 \pm 2\omega_3$ etc.
So the answer to my question is:
Cutting a signal in the time domain adds the higher harmonics with the combined frequencies in the frequency domain. Particularly, the combined frequencies allows to demodulate a signal in the crystal detector, i.e. shift them down (and up).
Thanks to all again!