How does hard clipping change the frequency of a pure sinusoidal signal? Suppose that we decide to limit the magnitude of a real-valued signal $f(t)$ by maximum cutoff $V_s$. Thus, if $|f(t)| > V_s$, a transformed signal $g(t) = V_s$ or $g(t) = -V_s$ depending on the sign of $f(t)$. If $|f(t)| \leq V_s$, $f(t) = g(t)$.
Now suppose $f(t) = A\cos (\omega t)$, with some real $\omega$ and real $A > V_s$, which is a pure cosine signal. Now clipping is done at cutoff $V_s$ as in above.
What happens to Fourier series of $g(t)$?
 A: In general, the Fourier series of a signal $f(t)$ (given all the conditions for such a Fourier transformation are satisfied) is represented as
$$ f_0 + \sum_{n=1}^{\infty}{A_n\cos(nt) + B_n \sin(nt)}$$
where $A_n = \frac{1}{\pi} \int_{0}^{2\pi}{f(t) \cos(nt) dt}$ and $B_n = \frac{1}{\pi} \int_{0}^{2\pi}{f(t) \sin(nt) dt}$
If an LTI (linear time-invariant) system is applied on a cosine wave, the output is also a cosine wave, possibly with different phase and amplitude, but no higher order harmonics. 
Now for the cosine wave passed through this clipping device, the top and bottoms of the output wave are clipped to $\pm V_s$, when $A > V_s$, we no longer have an LTI system so we must do the whole Fourier analysis on the output waveform to understand its harmonic content. This distorted wave has a Fourier series in which the $A_n$ coefficients are zero because the output wave is odd in $\omega t$. Thus, to find the Fourier transform for this nonlinearity, we only have to evaluate $B_n$ as follows:
$$  \begin{aligned}
B_n = \frac{1}{\pi} \left[ \int_{0}^{\beta}{V_s \sin(nt) dt} + 
      \int_{\beta}^{\pi-\beta}{\cos(\omega t) \sin(nt) dt} - \\
      \int_{\pi-\beta}^{\pi+\beta}{V_s \sin(nt) dt} + 
      \int_{\pi+\beta}^{2\pi-\beta}{\cos(\omega t) \sin(nt) dt} +
      \int_{2\pi-\beta}^{2\pi}{V_s \sin(nt) dt} \right]  
\end{aligned}
$$
where $\beta\omega = \arccos\left( \frac{V_s}{V}\right)$. I didn't try to simplify these terms to make a point that when you have to deal with a nonlinear system, you need to proceed with the calculation of Fourier coefficients and figure out how the nonlinear distortion modifies the integrand terms and integral boundaries. Unless some properties of the system impose particular Fourier coefficients to be zero, higher order harmonics are not zero in general.
