The electrical activity of the heart may be modeled by a non-linear system of partial differential equations known as the bidomain model. Due to the rapid variations in the electrical field, accurate simulations require a fine-scale discretization of the equations and consequently the solution of large severely ill-conditioned linear systems at each time step. Solving these systems is a major bottleneck of the whole simulation. We propose a highly effective preconditioning strategy for a general and popular 3D formulation of the problem. A theoretical analysis of the preconditioned matrix ensuring mesh independence of the spectrum is also described. Numerical comparisons with state of the art approaches confirm the effectiveness of our preconditioning technique. Finally, we show that an equivalent but less exercised formulation provides the best performance, in terms of CPU time.