A very nice way to see all of this is to look start with the Poisson semi-group $f \mapsto e^{-t\sqrt{-\Delta}}f$ for $t>0$. These operators are defined by Fourier Transform as $$\widehat{e^{-t\sqrt{-\Delta}}f}(p)=e^{-t|p|}\widehat{f}(p),$$ for any function $f$ so that the right hand side makes sense. Then for $f$, say in $ L^1$, we have $$e^{-t\sqrt{-\Delta}}f(x)=\int_{\mathbb{R}^n} P_t(x-y)f(y) dy$$ where $P_t(x)$ is the Poisson kernel $$P_t(x)=\frac{1}{(2\pi)^n} \int_{\mathbb{R}^n} e^{ix\cdot p} e^{-t|p|} dp = C_n \frac{t}{(t^2 +|x|^2)^\frac{n+1}{2}}. $$ (The computation of $P_t(x)$ is carried out in the first chapter of Stein and Weiss, Fourier Analysis on Euclidean Spaces.) The formula for $\sqrt{-\Delta}$ follows by taking the limit $$\sqrt{-\Delta}f = \lim_{t\downarrow 0} \frac{f- e^{-t\sqrt{-\Delta}}f}{t}$$ whenever this limit exists in a suitable sense. Because $\int P_t(x)dx=1$ we have $$\frac{1}{t} \left (f(x)-e^{-t\sqrt{-\Delta}}f(x) \right ) = C_n \int_{\mathbb{R}^n} \frac{f(x)-f(y)}{(t^2 +|x-y|^2)^{\frac{n+1}{2}}}dy.$$ Your identity now follows whenever $f$ is smooth enough and decays fast enough for the integral to make sense.
