A positive semidefinite Toeplitz matrix, which often arises as the finite covariance matrix of a stationary random process, can be decomposed as the sum of a nonnegative multiple of the identity corresponding to a white noise, and a singular term corresponding to a purely deterministic process. Moreover, the singular nonnegative Toeplitz matrix admits a unique characterization in terms of spectral lines which are associated to an oscillatory signal. This is the content of the famous Carathéodory-Fejér theorem. Its importance lies in the practice of extracting the signal component from noise, providing insights in modeling, filtering, and estimation. The multivariate counterpart of the theorem concerning block-Toeplitz matrices is less well understood, and in this paper, we aim to partially address this issue. To this end, we first establish an existence result of the line spectrum representation for a finite covariance multisequence of some underlying random vector field. Then, we give a sufficient condition for the uniqueness of the representation, which indeed holds true in the special case of bivariate time series. Equivalently, we obtain the Vandermonde decomposition for positive semidefinite block-Toeplitz matrices with $2\times 2$ blocks. The theory is applied to the problem of frequency estimation with two measurement channels within the recently developed framework of atomic norm minimization. It is shown that exact frequency recovery can be guaranteed in the noiseless case under suitable conditions, while in the noisy case, extensive numerical simulations are performed showing that the method performs well in a wide range of signal-to-noise ratios.