A generalized approach derived from Bloch's equation of motion of nuclear magnetic moments is presented to model the frequency, magnetic field, spin density, and temperature dependencies in the electromagnetic permeability tensor for materials with magnetic resonances. The resulting tensor model predicts characteristic polarization signatures which can be observed, for example, in Mueller matrix element spectra measured. When augmented with thermodynamic considerations and suitable Hamiltonian description of the magnetic eigenvalue spectrum, important parameters such as density, spectral amplitude distribution, relaxation time constants, and geometrical orientation parameters of the magnetic moments can be obtained from comparing the generalized model approach to experimental data. We demonstrate our approach by comparing model calculations with full Mueller matrix element spectra measured at an oblique angle of incidence in the terahertz spectral range, across electron spin resonance quintuplet transitions observed in wurtzite-structure GaN doped with iron. Our model correctly predicts the complexity of the polarization signatures observed in the 15 independent elements of the normalized Mueller matrix for both positive and negative magnetic fields and will become useful for future analysis of frequency and magnetic field-dependent magnetic resonance measurements.