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}$.)