The calculation of spectra of two coupled spin systems within the stochastic Liouville equation formalism is systematized by the use of the coupled multipole-operator representation. The Liouville superoperator of the system is easily calculated in this representation in a unified way for all interactions by using a generalized Wigner-Eckart theorem in operator space. The spectrum is defined as the Fourier-Laplace transform of a single generalized correlation function, which, besides the autocorrelation of the observed spin can also contain contributions of the unobserved spin degrees of freedom, as well as all counterrotating parts when needed. The expression for the spectrum, which is an average of the appropriate first-order multipoles, is the same irrespective of the particular system; i.e., the multiplicity of the spin system, its states, and the relevant transition probabilities are all treated implicitly. This definition of the spectrum also introduces, through the preparation stage of the spin system prior to the free induction decay, another first-order multipole which is important for the case of strongly correlated spins. Two spins are strongly correlated when the difference in their Zeeman resonance frequency is comparable to their broadening. Accordingly, strong correlation is most likely to be important in low field and/or slow motion. The present method imposes no limitations on the relative strength of the relevant interactions or the rate of their modulation and is consequently valid in the slow-motion regime. A fast matrix-tridiagonalization method; based on the Lanczos algorithm, is used to calculate the spectrum, working directly in the frequency domain.