  In this work we derive the covariant and gauge invariant perturbation
equations in general theories of $f(R)$ gravity in the Palatini formalism to
linear order and calculate the cosmic microwave background (CMB) and matter
power spectra for an extensively discussed model, $f(R) = R +
\alpha(-R)^\beta$, which is a possible candidate for the late-time cosmic
accelerating expansion found recently. These spectra are discussed and found to
be sensitively dependent on the value of $\beta$. We are thus able to make
stringent constraints on $\beta$ from cosmological data on CMB and matter power
spectra: The three-year Wilkinson Microwave Anisotropy Probe (WMAP) data alone
gives a constraint $|\beta| \lesssim \mathcal{O}(10^{{-3}})$ while the joint
WMAP, Supernova Lagacy Survey (SNLS) and Sloan Digital Sky Survey (SDSS) data
sets tightens this to $\beta \sim \mathcal{O}(10^{{-6}})$, about an order of
magnitude more stringent than the constraint from SDSS data alone, which makes
this model practically indistinguishable from the standard
$\Lambda\mathrm{CDM}$ paradigm.
