In classical physics, the canonical ensemble is defined over six-dimensional phase space $(\mathbf{x},\mathbf{p})$ only, i.e. three dimensions for coordinates and three dimensions for momentum. To answer your question of finding marginal energy density from full phase (state) space, it is sufficient to convert momentum density distribution $f(p)dp$ to the corresponding energy density distribution $f(E)dE$ because Maxwell-Boltzmann distribution has no dependence on $\textbf{x}$ space. Note that momentum density distribution is derived first, followed by energy density derivation.
Changing momemtum vector integral to scalar integral is similar to conversion from Cartesian coordinate integral to spherical coordinate integral:
$$
f(\textbf{p})dp_xdp_ydp_z = 4\pi p^2f(p)dp
$$
where the three-dimensional integral is reduced to one dimension integral by integrating away the two angular dimensions. The remaining work is just a change of variable using $E = p^2/2m$ (c.f. wikipedia for full expression).
In quantum mechanics, the continuous $\textbf{p}$ space is replaced by discretized momentum vector $\textbf{k}$ space. In one dimension, $k_x = 2\pi n/L_x$ as a consequence of periodic boundary condition. So, $\textbf{k}$-density in one dimension $$\rho(k)_x dk_x=\frac{n}{k_x}dk_x=\frac{L_x}{2\pi}dk_x$$. In three dimension,$$\rho(k)dk_xdk_ydk_z=\frac{L_xL_yL_z}{(2\pi)^3}dk_xdk_ydk_z=\frac{V}{(2\pi)^3}dk_xdk_ydk_z$$.
We can convert the "Cartesian" k-space into "spherical" form. Integrating over all angles, the magnitude of $\textbf{k}$-density becomes
$$\rho(k)dk_xdk_ydk_z=\frac{V}{(2\pi)^3}4\pi k^2dk$$
The last equation above gives the momentum vector $\textbf{k}$-density in quantum mechanical form. Finally, conversion to energy density can be obtained if you have the equation that relates $E$ and $k$:
$$\rho(E) = \frac{V}{(2\pi)^3}4\pi k^2\left(\frac{dk}{dE}\right)dE$$