We calculate the finite-temperature r-mode spectrum of a superfluid neutron star accounting for both muons in the core and the entrainment between neutrons and protons. We show that the standard perturbation scheme, considering the rotation rate as an expansion parameter, breaks down in this case. We develop an original perturbation scheme which circumvents this problem by treating both the perturbations due to rotation and (weak) entrainment simultaneously. Applying this scheme, we propose a simple method for calculating the superfluid r-mode eigenfrequency in the limit of vanishing rotation rate. We also calculate the r-mode spectrum at finite rotation rate for realistic microphysics input (adopting, however, the Newtonian framework and Cowling approximation when considering perturbed oscillation equations) and show that the normal r-mode exhibits resonances with superfluid r-modes at certain values of temperatures and rotation frequencies in the parameter range relevant to neutron stars in low-mass x-ray binaries (LMXBs). This turns the recently suggested phenomenological model of resonance r-mode stabilization into a quantitative theory, capable of explaining observations. A strong dependence of resonance rotation rates and temperatures on the neutron superfluidity model allows us to constrain the latter by confronting our calculations with the observations of neutron stars in LMXBs.