The seismic response of inhomogeneous soils is explored analytically by means of one-dimensional viscoelastic wave propagation theory. The system under investigation comprises of a continuously inhomogeneous layer over a homogeneous one of higher stiffness. The excitation is specified at the bottom of the base layer in the form of vertically propagating harmonic S waves. Shear wave propagation velocity in the inhomogeneous layer is described by a generalized parabolic function, which allows modeling of soil having vanishing shear modulus at the ground surface. The problem is treated analytically leading to an exact solution of the Bessel type for the natural frequencies, mode shapes and base-to-surface response transfer function. The model is validated using available theoretical solutions and finite-element analyses. The exact analytical solution is compared with energybased Rayleigh techniques and equivalent homogeneous soil approximations. The latter are defined by means of alternative definitions for the representative shear wave velocity in the inhomogeneous layer. Results are presented in the form of normalized graphs demonstrating the effect of salient model parameters such as layer thickness, impedance contrast between surface and base layer, surface-to-base shear wave velocity ratio in the inhomogeneous layer, rate of inhomogeneity and hysteretic damping ratio. Harmonic response of inhomogeneous soils with vanishing shear wave velocity near soil surface is explored by asymptotic analyses.