Obtaining channel covariance knowledge is of great importance in various Multiple-Input Multiple-Output MIMO communication applications, including channel estimation and covariance-based user grouping. In a massive MIMO system, covariance estimation proves to be challenging due to the large number of antennas ($M\gg 1$) employed in the base station and hence, a high signal dimension. In this case, the number of pilot transmissions $N$ becomes comparable to the number of antennas and standard estimators, such as the sample covariance, yield a poor estimate of the true covariance and are undesirable. In this paper, we propose a Maximum-Likelihood (ML) massive MIMO covariance estimator, based on a parametric representation of the channel angular spread function (ASF). The parametric representation emerges from super-resolving discrete ASF components via the well-known MUltiple SIgnal Classification (MUSIC) method plus approximating its continuous component using suitable limited-support density function. We maximize the likelihood function using a concave-convex procedure, which is initialized via a non-negative least-squares optimization problem. Our simulation results show that the proposed method outperforms the state of the art in various estimation quality metrics and for different sample size to signal dimension ($N/M$) ratios.