Integration of radial functions?

Let $f(|x|)$ be a integrable radial function in $\mathbb{R}^n$ ($|\cdot|$ denotes the euclidean norm as in convention). The following identity is used to simplify computations

$$\int_{\mathbb{R}^n}f(|x|)\mathrm{d}x = \omega_{n-1}\int_0^\infty f(r) r^{n-1} \mathrm{d}r,$$ where $\omega_{n-1}$ denotes the surface area of the $(n-1)$-sphere of radius 1.

But I have never found any full proofs on such topic. Does anyone have some references in mind, better some source which I can actually refer to?

updates

Recently I found this subject has been mentioned by Stein in the appendix from his book "Fourier Analysis", and it further refers to Buck's "Advanced Calculus", Folland's "Advanced Calculus" and Lang's "Undergraduate Analysis". I have gone through the later two. It seems to me that Lang's treatment is pretty self-contained, though it is given in exercise hence details omitted.

Mentioned by TCL, Evan's book on measure theory is decent but I'm not sure if it's too heavy. It let me feel like shooting mosquito with bazooka.

Spherical coordinates in $\mathbb{R}^n$ is treated fully in Karl Stromberg's book "Introduction to Classical Real Analysis" p. 369-370. Note that on p. 370, it shows the Jacobian of the coordinate transformation as $$r^{n-1} \prod_{j=1}^{n-2} (\sin\theta_j)^{n-j-1}$$ So if your function is radial, then \begin{align*}\int_{\mathbb{R}^n}f(|x|)\mathrm{d}x &=\int_{-\pi}^\pi\int_0^\pi\int_0^\pi\cdots\int_0^\infty f(r)r^{n-1}\prod_{j=1}^{n-2} (\sin\theta_j)^{n-j-1}\,dr \, d\theta_1\cdots d\theta_{n-2}\,d\theta_{n-1}\\&= \omega_{n-1}\int_0^\infty f(r) r^{n-1} \mathrm{d}r,\end{align*} where $$\omega_{n-1}=\int_{-\pi}^\pi\int_0^\pi\int_0^\pi\cdots\int_0^\pi \prod_{j=1}^{n-2} (\sin\theta_j)^{n-j-1}\,d\theta_1\cdots d\theta_{n-2}\,d\theta_{n-1}$$

Stromberg, Karl Robert, An introduction to classical real analysis, Providence, RI: AMS Chelsea Publishing (ISBN 978-1-4704-2544-9/hbk). xiv, 577 p. (2015). ZBL1331.00003.

There are a couple good answers already, however I thought I would add a more elementary approach, avoiding the (direct) use of measures or spherical coordinates. The price we pay here is that we must assume $f$ to be continuous, as opposed to merely integrable. (We can get around this, since one can approximate integrable functions with continuous functions in the $L^1$-norm.)

Claim: If $f$ is a continuous real valued function on $\mathbb{R}$, then $$\int_{\mathbb{R}^n}f(|x|)\,dx=\omega_{n-1}\int_0^{\infty}f(r)r^{n-1}\,dr.$$

Proof: Fix $R>0$, set $B_R:=B(0,R)=\{|x|<R\}$. For each fixed integer $N$, one can partition $B_R$ into $N$ sub-annuli $B^i_R\subset B_R$, $1\le i\le N$, where each $B^i_R=\{r_{i-1}\le|x|<r_i\}$, and $r_0=0<r_1<\cdots<r_{N-1}<r_N=R$. We will make the partition size equal, and denote it by $\delta_N:=R/N$. Note that for each $i\in\{1,\dots,N\}$, the mean value theorem furnishes a number $s_i\in(r_{i-1},r_i)$ such that $$r_{i}^n-r_{i-1}^n=ns_i^{n-1}(r_{i}-r_{i-1}).$$ With this in hand (and keeping in mind we are working with the integral of a radial function), we define a new function $f_N:B_R\to\mathbb{R}$ by setting $$f_N(x)=f(s_i),$$ whenever $x\in B^i_R$. We wish to show that $f_N(\cdot)\to f(|\cdot|)$ uniformly on $B_R$ as $N\to\infty$. Fix $\epsilon>0$ and choose $N$ so large that $\delta_N<\delta$, where $\delta>0$ is such that $|f(r)-f(s)|<\epsilon$ whenever $r,s\in B_R$ with $|r-s|<\delta$ (this $\delta$ exists by the uniform continuity of $f$ on $B_R$.) Note then that for any $x\in B_R$, then $x\in B^i_R$ for some $i$ and we have $$|f_N(x)-f(|x|)|=|f(s_i)-f(|x|)|<\epsilon.$$ Therefore $f_N(\cdot)\to f(|\cdot|)$ uniformly on $B_R$ as $N\to\infty$ and hence we have $$\int_{B_R}f_N(x)\,dx\to\int_{B_R}f(|x|)\,dx$$ as $N\to\infty.$ Note that the volume of $B^i_R$ is equal to $\omega_{n-1}(r^n_i-r^n_{i-1})=\omega_{n-1}s_i^{n-1}(r_i-r_{i-1})=\omega_{n-1}s_i^{n-1}\delta_N$. Now, using this and the definition of Riemann integration, we can compute $$\int_{B_R}f_N(x)\,dx=\sum_{i=1}^Nf(s_i)\omega_{n-1}s^{n-1}_i\delta_N\to\omega_{n-1}\int_0^Rf(r)r^{n-1}\,dr,$$ as $N\to\infty$. Hence we have $$\int_{B_R}f(|x|)\,dx=\omega_{n-1}\int_0^Rf(r)r^{n-1}\,dr.$$ Letting $R\to\infty$ gives the desired result.

Since $\mathbb{R}^n=\{0\}\cup (0,\infty)\times S_r^{n-1}$, its volume element $dV$ can be written as a the product of $dr$ times the volume element of the $n-1$ sphere of radius $r$. By dimensional grounds, the latter is equal to $r^{n-1}$ times the volume element of the unit $(n-1)$-sphere. For example, if $n=2$ one has $dV=r dr\wedge d\theta$, note that the integral of $d\theta$ over $S^1$ is $2\pi$, i.e. the "volume" of $S^1$; similarly for $n=2$, $dV=r^2 \sin\theta\, dr\wedge d\theta$ and $\int _{S^2} \sin\theta d\theta\wedge d\phi=4\pi$. Since your function is radial, you can factorize the integral into a radial and an angular part. The integral of the angular part, which you denote with $\omega_{n-1}$, is therefore the volume of the unit $(n-1)$-sphere. A nice way of calculating $\omega_{n-1}$ is presented here http://scipp.ucsc.edu/~haber/ph116A/volume_11.pdf

Proof using coarea formula for general $\ell_p$-norm radial functions

Let $p \in [1,\infty]$ and let $\|\cdot\|_p$ denote the $\ell_p$-norm on $\mathbb R^n$. In particular, $p=2$ corresponds to the euclidean case in the OP's question. I will handle the general case.

Let $B_{n,p}(t) := \{x \in \mathbb R^p \mid \|x\|_p \le t\}$ be the $\ell_p$-ball of radius $t \ge 0$ in $\mathbb R^n$. Let $ :[0, \infty) \rightarrow \mathbb R$ be a measurable function, and define measurable functions $g,u:\mathbb R^n \rightarrow \mathbb R$ by $g(y) := f(\|y\|_p)$, $r(y) := \|y\|_p$. The $\omega_{n,p}(t)$ be the surface area of the level-set $r^{-1}(t) := \{y \in \mathbb R^n \mid \|y\|_p = t\}=:S_{n,p}(0,t) = \partial B_{n,p}(t)$. Note that thanks to the triangle inequality for $\ell_p$-norms, $y \mapsto r(y)$ is a Lipschitz function on $\mathbb R^n$ and the level-set $r^{-1}(t)$ is a.e smooth for all $t \ge 0$. Finally, note that $\partial_j r(y) = \dfrac{y_j|y_j|^{p-2}}{\|y\|_p^{p-1}}$ and so for every $y \in \mathbb R^n\setminus\{0\}$, we have $$ \begin{split} \|\nabla r(y)\|_2 &= \frac{1}{\|y\|_p^{p-1}}\sqrt{\sum_{j=1}^n y_j^2y_j^{2(p-2)}}=\frac{1}{\|y\|_p^{p-1}}\sqrt{\sum_{j=1}^n y_j^{2(p-1)}}\\ &=a_{n,p}(y) := \begin{cases}1, &\mbox{ if }p \in \{2,\infty\},\\ \text{ some explicit expression},&\mbox{ else.}\end{cases} \end{split} $$ For $m \in \mathbb N$, let $\mathcal L^m$ be the $m$-dimensional Lebesgue measure (aka $m$-dimensional volume) and let $\mathcal H^m$ be the $m$-dimensional Hausdorff measure (aka $m$-dimensional surface area). Then by the coarea formula, we have

$$ \begin{split} \int_{\mathbb R^n}a_{n,p}(y)f(\|y\|_p)dy &= \int_{\mathbb R^n}g(y)\|\nabla r(y)\|_2dy = \int_{\mathbb R}\left(\int_{r^{-1}(t)}g(y)d\mathcal H^{n-1}(y)\right)dt\\ &= \int_{0}^\infty f(t)\mathcal H^{n-1}(S_{n,p}(0,t))dt = \omega_{n,p}(1)\int_{0}^\infty f(t)t^{n-1}dt, \end{split} $$

Thus, we get the identity

$$ \int_{\mathbb R^n} a_{n,p}(y) f(\|y\|_p)dy = \omega_{n,p}(1)\int_{0}^\infty f(t)t^{n-1}dt. $$