My approach:
In general we have:
$$\delta(x-x_0)=\frac{1}{2\pi}\int_{\mathbb{R}} e^{-it(x-x_0)}\text{d}t$$
Then ($y \ge 0$):
$$I=\int_{\mathbb{R}^n} \delta(x_1^2+\dots+x_n^2-y)\text{d}x_1\dots\text{d}x_n= \\ =\int_{\mathbb{R}^n} \left[\frac{1}{2\pi}\int_{\mathbb{R}} e^{-it(x_1^2+\dots+x_n^2-y)}\text{d}t\right]\text{d}x_1\dots\text{d}x_n= \\ =\frac{1}{2\pi}\int_{\mathbb{R}} e^{ity}\left[\int_{\mathbb{R}^n} e^{-it(x_1^2+\dots+x_n^2)}\text{d}x_1\dots\text{d}x_n\right]\text{d}t= \\ =\frac{1}{2\pi}\int_{\mathbb{R}} e^{ity}\left[\prod_{i=1}^n \int_{-\infty}^{+\infty} e^{-itx_i^2}\text{d}x_i\right]\text{d}t$$
Now I use the fact that:
$$\int_{-\infty}^{+\infty} e^{-Ax^2+Bx}\text{d}x=\sqrt{\frac{\pi}{A}}e^{\frac{B^2}{4A}}$$
So we can write:
$$I=\frac{1}{2\pi}\int_{\mathbb{R}} e^{ity}\prod_{i=1}^n \sqrt{\frac{\pi}{it}}\text{d}t= \\ =\frac{1}{2\pi}\pi^{n/2}i^{-n/2}\int_{\mathbb{R}} t^{-n/2}e^{ity}\text{d}t$$
Here I got stuck. Is my approach correct? Thank you!
Note that
$$ \delta(x) = - \frac{\mathrm{d}}{\mathrm{d}\epsilon}\biggr|_{\epsilon=0} \mathbf{1}[x+\epsilon \leq 0]$$
in distribution sense. Using this, we first compute
\begin{align*} F(\epsilon) &= - \int_{\mathbb{R}^n} \mathbf{1}[x_1^2 + \cdots + x_n^2 - y + \epsilon \leq 0] \, \mathrm{d}x_1\cdots\mathrm{d}x_n. \end{align*}
Note that this is the negative of the volume of an $n$-dimensional ball of radius $\sqrt{y-\epsilon}$. So, using the formula for the volume of an $n$-ball, we get
\begin{align*} F(\epsilon) &= - \frac{\pi^{n/2}}{\Gamma(\frac{n}{2}+1)}(y - \epsilon)^{n/2}, \end{align*}
provided $\epsilon < y$. So, if $y > 0$, then taking derivative at $\epsilon = 0$ gives
$$ I = F'(0) = \frac{\pi^{n/2}}{\Gamma(\frac{n}{2})} y^{\frac{n}{2}-1}. $$