The problem of inverting the very large system of equations implied by the variational principle defining a statistical objective analysis scheme is normally tackled by an iterative numerical method, such as the suitably preconditioned conjugate-gradient or quasi-Newton techniques. At each iteration, there is a requirement to convolve the background error covariance with a gridded field in order to recover another gridded field. Essentially, the background covariance becomes a filtering operator. This particular step, while it is only one of the several algebraic steps that constitute a single iteration cycle, is typically the one that dominates the computational cost in 3D-variational analysis. It is therefore very important to be able to execute it efficiently. It is also a prerequisite for the numerical stability of the solvers that the representation of the covariance operator retains the properties of self-adjointness and positive-definiteness possessed by actual covariances. The synthesis of covariance operators by carefully constructed combinations of simpler filters can be an efficient and versatile approach, particularly if the basic filters act in one dimension at a time. In order to make this approach effective without the principal grid directions imposing a spurious anisotropic imprint on the morphology of the synthesized covariance function, the basic one-dimensional filters need to approximate a Gaussian profile closely. Recursive filters are able to mimic this profile efficiently and accurately, as we shall demonstrate. The covariance profile itself is not restricted to be of Gaussian type since, by superposition of a few Gaussian components of different scales, it is possible to produce a range of `fat-tailed' covariance profiles that are arguably better suited to practical data assimilation. Also, the application of the negative-Laplacian operator to an originally bell-shaped function yields a profile possessing negative side-lobes, which are desirable in the analysis of some variables. Recent developments of the recursive filtering technique have allowed us the freedom to consistently extend the synthesis of covariances to incorporate geographically adaptive modifications of covariance scale and amplitude, and general degrees of local stretching and compression in both two and three dimensions at arbitrary orientations. We present some of these new techniques and describe the numerical approaches adopted to implement these schemes on a massively parallel computer.