In this paper, we study the problem of multi-band (frequency-variant)covariance interpolation with a particular emphasis towards massive MIMOapplications. In a massive MIMO system, the communication between each BS with$M \gg 1$ antennas and each single-antenna user occurs through a collection ofscatterers in the environment, where the channel vector of each user at BSantennas consists in a weighted linear combination of the array responses ofthe scatterers, where each scatterer has its own angle of arrival (AoA) andcomplex channel gain. The array response at a given AoA depends on thewavelength of the incoming planar wave and is naturally frequency dependent.This results in a frequency-dependent distortion where the second orderstatistics, i.e., the covariance matrix, of the channel vectors varies withfrequency. In this paper, we show that although this effect is generallynegligible for a small number of antennas $M$, it results in a considerabledistortion of the covariance matrix and especially its dominant signal subspacein the massive MIMO regime where $M \to \infty$, and can generally incur aserious degradation of the performance especially in frequency divisionduplexing (FDD) massive MIMO systems where the uplink (UL) and the downlink(DL) communication occur over different frequency bands. We propose a novelUL-DL covariance interpolation technique that is able to recover the covariancematrix in the DL from an estimate of the covariance matrix in the UL under amild reciprocity condition on the angular power spread function (PSF) of theusers. We analyze the performance of our proposed scheme mathematically andprove its robustness under a sufficiently large spatial oversampling of thearray. We also propose several simple off-the-shelf algorithms for UL-DLcovariance interpolation and evaluate their performance via numericalsimulations.