In the covariant cosmological perturbation theory, a 1+3 decomposition ensures that all variables in the frame-independent equations are covariant, gauge-invariant and have clear physical interpretations. We develop this formalism in the case of Brans-Dicke gravity, and apply this method to the calculation of cosmic microwave background (CMB) anisotropy and large scale structures (LSS). We modify the publicly available Boltzmann code CAMB to calculate numerically the evolution of the background and adiabatic perturbations, and obtain the temperature and polarization spectra of the Brans-Dicke theory for both scalar and tensor mode, the tensor mode result for the Brans-Dicke gravity are obtained numerically for the first time. We first present our theoretical formalism in detail, then explicitly describe the techniques used in modifying the CAMB code. These techniques are also very useful to other gravity models. Next we compare the CMB and LSS spectra in Brans-Dicke theory with those in the standard general relativity theory. At last, we investigate the ISW effect and the CMB lensing effect in the Brans-Dicke theory. Constraints on Brans-Dicke model with current observational data is presented in a companion paper (paper II).