We consider the cosmology where some function $f(G)$ of the Gauss-Bonnet term $G$ is added to the gravitational action to account for the late-time accelerating expansion of the universe. The covariant and gauge invariant perturbation equations are derived with a method which could also be applied to general $f(R,{R}^{ab}{R}_{ab},{R}^{abcd}{R}_{abcd})$ gravitational theories. It is pointed out that, despite their fourth-order character, such $f(G)$ gravity models generally cannot reproduce arbitrary background cosmic evolutions; for example, the standard $\ensuremath{\Lambda}\mathrm{CDM}$ paradigm with ${\ensuremath{\Omega}}_{\mathrm{DE}}=0.76$ cannot be realized in $f(G)$ gravity theories unless $f$ is a true cosmological constant because it imposes exclusionary constraints on the form of $f(G)$. We analyze the perturbation equations and find that, as in the $f(R)$ model, the stability of early-time perturbation growth puts some constraints on the functional form of $f(G)$, in this case ${\ensuremath{\partial}}^{2}f/\ensuremath{\partial}{G}^{2}l0$. Furthermore, the stability of small-scale perturbations also requires that $f$ not deviate significantly from a constant. These analyses are illustrated by numerically propagating the perturbation equations with a specific model reproducing a representative $\ensuremath{\Lambda}\mathrm{CDM}$ cosmic history. Our results show how the $f(G)$ models are highly constrained by cosmological data.