We give analytical description of generation of harmonics of the undulator radiation (UR) with account for the finite electron beam size, emittance, off-axis beam deviation and electron energy spread, as well as for the constant magnetic components and field harmonics effects. We give exact analytical expressions for the generalized Bessel and Airy functions, which describe the spectrum line shape and intensities in the two-frequency bi-harmonic undulator with account for the above factors. The obtained analytical formulae distinguish contributions of each field component and every undulator and beam parameter on the harmonic radiation in free electron lasers (FEL). The effect of the field on the harmonic radiation is analyzed with account for the beam finite size and its off-axis deviation. The phenomenological model is employed for the FEL modeling; with its help we study the harmonic generation, including even ones, in the experiments LCLS and LEUTL. We demonstrate analytically that strong second FEL harmonic in X-ray FEL at the wavelengths λ = 0.75nm in the LCLS experiment is caused by the deviation of the electron trajectories off the axis in 15 μm on the gain length 1.6 m, which is comparable with the beam size; the strong second FEL harmonic in the LEUTL experiment at the wavelength λ = 192nm can be attributed to interaction of the electrons in wide, ~ 0.2 mm, beam with the photon radiation at the gain length 0.87 m. The modeling results fully agree with the measurements. The developed formalism allows the analysis of projected and built FELs and their radiation, helps minimizing losses and correcting magnetic fields; it also shows physical background and reasons for each harmonic radiated power in the FEL.