Detailed knowledge of the primordial power spectrum of curvature perturbations is essential both in order to elucidate the physical mechanism (`inflation) which generated it, and for estimating the cosmological parameters from observations of the cosmic microwave background and large-scale structure. Hence it ought to be extracted from such data in a model-independent manner, however this is difficult because relevant cosmological observables are given by a convolution of the primordial perturbations with some smoothing kernel which depends on both the assumed world model and the matter content of the universe. Moreover the deconvolution problem is ill-conditioned so a regularisation scheme must be employed to control error propagation. We demonstrate that `Tikhonov regularisation can robustly reconstruct the primordial spectrum from multiple cosmological data sets, a significant advantage being that both its uncertainty and resolution are then quantified. Using Monte Carlo simulations we investigate several regularisation parameter selection methods and find that generalised cross-validation and Mallows $C_p$ method give optimal results. We apply our inversion procedure to data from the Wilkinson Microwave Anisotropy Probe, other ground-based small angular scale CMB experiments, and the Sloan Digital Sky Survey. The reconstructed spectrum (assuming the standard $Lambda$CDM cosmology) is emph{not} scale-free but has an infrared cutoff at $k lesssim 5 times 10^{-4}; mathrm{Mpc}^{-1}$ (due to the anomalously low CMB quadrupole) and several features with $sim 2 sigma$ significance at $k/mathrm{Mpc}^{-1} sim$ 0.0013--0.0025, 0.0362--0.0402 and 0.051--0.056, reflecting the `WMAP glitches. To test whether these are indeed real will require more accurate data, such as from the Planck satellite and new ground-based experiments.