Recent methodological advances have made it possible to measure fine structure on the order of a few hertz in the nuclear magnetic resonance (NMR) spectra of quadrupolar nuclei in polycrystalline samples. Since quadrupolar couplings are often a significant fraction of the Zeeman coupling, a complete analysis of such experimental spectra requires a theoretical treatment beyond first-order. For multiple pulse NMR experiments, which may include sample rotation, the traditional density matrix approaches for treating higher-order effects suffer from the constraint that undesired fast oscillations (i.e., multiples of the Zeeman frequency), which arise from allowed overtone transitions, can only be eliminated in numerical simulations by employing sampling rates greater than 2I times the Zeeman frequency. Here, we present a general theoretical approach for arbitrary spin I that implements an analytical “filtering” of undesired fast oscillations in the rotating tilted frame, while still performing an exact diagonalization. Alternatively, this approach can be applied using a perturbation expansion for the eigenvalues and eigenstates, such that arbitrary levels of theory can be explored. The only constraint in this approach is that the Zeeman interaction remains the dominant interaction. Using this theoretical framework, numerical simulations can be implemented without the need for a high sampling rate of observables and with significantly reduced computation times. Additionally, this approach provides a general procedure for focusing on the excitation and detection of both fundamental and overtone transitions. Using this approach we explore higher-order effects on a number of sensitivity and resolution issues with NMR of quadrupolar nuclei.