We introduce a multigrid technique for the solution of multilevel circulant linear systems whose coefficient matrix has eigenvalues of the form $f(x_j^{[n]})$, where f is continuous and independent of $n=(n_1,\ldots,n_d)$, and $x_j^{[n]} \equiv 2\pi j/n = (2\pi j_1/n_1, \ldots, 2\pi j_d/n_d)$, $0 \le j_r \le n_r - 1$. The interest of the proposed technique pertains to the multilevel banded case, where the total cost is optimal, i.e., $O(N)$ arithmetic operations (ops), $N=\prod_{r=1}^d n_r$, instead of $O(N\log N)$ ops arising from the use of FFTs. In fact, multilevel banded circulants are used as preconditioners for elliptic and parabolic PDEs (with Dirichlet or periodic boundary conditions) and for some two-dimensional image restoration problems where the point spread function (PSF) is numerically banded, so that the overall cost is reduced from $O(k(\varepsilon,n)N \log N)$ to $O(k(\varepsilon,n)N)$, where $k(\varepsilon,n)$ is the number of PCG iterations to reach the solution within an accuracy of $\varepsilon$. Several numerical experiments concerning one-rank regularized circulant discretization of elliptic $2q$-differential operators over one-dimensional and two-dimensional square domains with mixed boundary conditions are performed and discussed.