I think the constant in the equality you wrote down depends on how exactly you normalize the Fourier transform. Up to worrying about the constant, here's a proof that should work and may satisfy your criterion:

In dimension $d$, write $$ \int \frac{\widehat{f(\xi)}}{ |\xi|^{2}} d\xi = \int_0^\infty ds \int \hat{f}(\xi) e^{-s |\xi|^2} d\xi $$ Up to the normalization constant involving $(2 \pi)^d$ and justifying the convergence, the right hand side equals (by Plancharel) $$ \int_0^\infty ds \int f(x) e^{- \frac{|x|^2}{4 s} } (2s)^{-d/2} dx $$ Interchange the integrals and change the s variable such that $\tilde{s}^2 = (4s)^{-1}|x|^2$; the above equals $$ C \int f(x) |x|^{-d + 2} dx \int_0^\infty ~ e^{-\tilde{s}^2} \tilde{s}^{d - 3} d\tilde{s} $$ for an explicit constant $C$. Note that the power $d - 3 > -1$ since we are in dimension $d > 2$. The value of the $d\tilde{s}$ integral is then known and that will give your formula.

I think the constant in the equality you wrote down depends on how exactly you normalize the Fourier transform. Up to worrying about the constant, here's a proof that should work and may satisfy your criterion:

In dimension $d$, write $$ \int \frac{\widehat{f(\xi)}}{ |\xi|^{2}} d\xi = \int_0^\infty ds \int \hat{f}(\xi) e^{-s |\xi|^2} d\xi $$ Up to the normalization constant involving $(2 \pi)^d$ and justifying the convergence, the right hand side equals (by Plancharel) $$ \int_0^\infty ds \int f(x) e^{- \frac{|x|^2}{4 s} } (2s)^{-d/2} dx $$ Interchange the integrals and change the s variable such that $\tilde{s}^2 = (4s)^{-1}|x|^2$; the above equals $$ C \int \frac{f(x)}{|x|^{d - 2}} dx \int_0^\infty ~ e^{-\tilde{s}^2} \tilde{s}^{d - 3} d\tilde{s} $$ for an explicit constant $C$. Note that the power $d - 3 > -1$ since we are in dimension $d > 2$. The value of the $d\tilde{s}$ integral is then known and that will give your formula. (Equivalently, this answer gives a way to compute the Fourier transform of $|\xi|^{-2}$.)