Disclaimer. Below, I give a solution for the case of the $d$-dimensional gaussian distribution. I'm not 100% sure of all my arguments. I've hand checked with a few choices of the function $f$ (e.g linear $f(x) \equiv x^\top w$, etc.). Just posting here to start a discussion.
One can always write $\|\nabla f(x)\|^2 d\mu(x) = H(x)^2dx$, where $H(x): = F(x)G(x)$, with $F(x):=\|\nabla f(x)\|$ and $G(x):=\sqrt{p(x)}$.
By the Plancherel Theorem, we may simplify like so
$$
E_\mu(f) = \int_{\mathbb R^d} H(x)^2 dx = \int_{\mathbb R^d} |\hat{H}(z)|^2dz,
$$
where $\hat{H}$ is the Fourier transform of $H$. Now, the convolution property of the Fourier transform, we have $\hat{H} = \hat{F} \star \hat{G}:z \mapsto \int_{\mathbb R^d} \hat{F}(t)\hat{G}(t-z)dt$.
On the other hand, note that $\hat{F}(z) = \widehat{\|\nabla f\|}(z) = -i\|z\| \hat{f}(z)$. It remains to compute the Fourier transform of $G := \sqrt{p}$, and simplify...
Partial answer for gaussian case
In this case, the density of $\mu$ is $p = \gamma_{d,x}(x)$ defined by $\gamma_{d,\sigma}(x) = (2\pi\sigma^2)^{-d/2} e^{-\|x\|^2/(2\sigma^2)}$. This gives $\hat{G}(z) = (2\pi)^{-d/2}e^{-\sigma^2\|z\|^2/2}=\sigma^d\gamma_{d,1/\sigma}(z)$. Thus, $\hat{H} = T_{K_{d,1/\sigma}}\widehat{\|\nabla f\|} = T_{K_{d,1/\sigma}} \tilde{f}$, where $T_{K_{d,1/\sigma}}$ is the kernel integral operator corresponding to the psd radial kernel $$
K_{d,1/\sigma}(z,z'):= \gamma_{d,1/\sigma}(z-z'),
$$
and the function $\tilde{f}:\mathbb R^d \to \mathbb C$ is defined by
$$
\tilde{f}(z) := \|z\| \hat{f}(z).
$$
Because $T_{K_{d,1/\sigma}}$ is an isometry on $L^2(\gamma_{d,1/\sigma}) := L^2(\mathbb R^d,\gamma_{d,1/\sigma})$ (due to spherical symmetry), we deduce that
$$
E_{\gamma_{d,\sigma}}(f) = \sigma^d \|\tilde{f}\|_{L^2(\gamma_{d,1/\sigma})}^2.
$$