Analysis of data generated by Mueller matrix polarimeters using two photoelastic modulators has been evolving with the improvements in data acquisition and digital signal processing (DSP). Historical processing of the temporal data generated by these devices has involved isolating the frequencies via hardware signal processing (e.g., lock-in amplifiers) or the numerical computation of Fourier integrals of recorded temporal data. Both avenues have their advantages, but the DSP aspects of the latter provide greater flexibility in choice of harmonics for processing. While conventional processing uses one harmonic for each desired Mueller matrix element, recent work has demonstrated that theoretical improvements are possible by coherently combining the information in multiple harmonic channels for each element. We demonstrate some recent progress in DSP that enables these polarimeters' data to be more fully exploited by addressing two key issues in the Fourier domain: spectral leakage and phase recovery. Adequately addressing these issues enables numerical analysis of the temporal data in the complex Fourier domain and delivers Mueller matrix results in which spectral phase information is used to recover the matrix elements and determine their signs automatically. We explore the application of this complex analysis and how the precision and accuracy of the results are affected by common experimental and DSP limitations compared to the usual magnitude-only analysis in the Fourier domain. The multi-harmonic method can provide a theoretical factor of 1.3-1.7 improvement in instrumental precision, and our experimental results approach that theoretical range.