A multi-configuration mixing approach built on essentially complex, symmetry-projected Hartree-Fock-Bogoliubov (HFB) mean fields is introduced. The mean fields are obtained by variation after projection. The configuration space consists out of the symmetry-projected HFB vacuum and the symmetry-projected two-quasiparticle excitations for even, and the symmetry-projected one-quasiparticle excitations for odd A systems. The underlying complex HFB transformations are assumed to be time-reversal invariant and axially symmetric. The model allows nuclear structure calculations in large model spaces with arbitrary two-body interactions. The approach has been applied to $^{20}$Ne and $^{22}$Ne. Good agreement with the exact shell model results and considerable improvement with respect to older calculations, in which only real HFB transformations were admitted, is obtained.