This paper introduces the formalism which connects between rotation measure ( { RM } ) measurements for extragalactic sources and the cosmological magnetic field power spectrum . It is shown that the amplitude and shape of the cosmological magnetic field power spectrum can be constrained by using a few hundred radio sources , for which Faraday { RM } s are available . This constraint is of the form { { { { B _ { rms } \mathrel { \mathchoice { \lower 0.5 pt \vbox { \halign { \cr } $ \displaystyle \hfil < % $ \cr$ \displaystyle \hfil \sim$ } } } { \lower 0.5 pt \vbox { \halign { \cr } $ \textstyle \hfil% < $ \cr$ \textstyle \hfil \sim$ } } } { \lower 0.5 pt \vbox { \halign { \cr } $ \scriptstyle \hfil% < $ \cr$ \scriptstyle \hfil \sim$ } } } { \lower 0.5 pt \vbox { \halign { \cr } $% \scriptscriptstyle \hfil < $ \cr$ \scriptscriptstyle \hfil \sim$ } } } } 1 \times [ 2.6 \times 1 % 0 ^ { -7 } { cm } ^ { -3 } / \bar { n } _ { b } ] h nano-Gauss ( nG ) on \sim 10 - 50 { h ^ { -1 } Mpc } scales , with \bar { n } _ { b } the average baryon density and h the Hubble parameter in 100 km s ^ { -1 } Mpc ^ { -1 } units . The constraint is superior to and supersedes any other constraint which come from either CMB fluctuations , Baryonic nucleosynthesis , or the first two multipoles of the magnetic field expansion . The most adequate method for the constraint calculation uses the Bayesian approach to the maximum likelihood function . I demonstrate the ability to detect such magnetic fields by constructing simulations of the field and mimicking observations . This procedure also provides error estimates for the derived quantities . The two main noise contributions due to the Galactic RM and the internal RM are treated in a statistical way following an evaluation of their distribution . For a range of magnetic field power spectra with power indices -1 \leq n \leq 1 in a flat cosmology ( \Omega _ { m } =1 ) we estimate the signal-to-noise ratio , Q , for limits on the magnetic field B _ { rms } on \sim 50 { h ^ { -1 } Mpc } scale . Employing one patch of a few square degrees on the sky with source number density n _ { src } , an approximate estimate yields Q \simeq 3 \times ( B _ { rms } / 1 { nG } ) ( n _ { src } / 50 { deg } ^ { -2 } ) ( 2.6 \times 10 ^ { - % 7 } { cm } ^ { -3 } / \bar { n } _ { b } ) h . An all sky coverage , with much sparser , but carefully tailored sample of \sim 500 sources , yields Q \simeq 1 with the same scaling . An ideal combination of small densely sampled patches and sparse all-sky coverage yields Q \simeq 3 with better constraints for the power index . All of these estimates are corroborated by the simulations .