How Is the Formula for the Wiener Deconvolution Derived?
The Wiener Filter can also be derived by another (Easier) way.

Let's assume the following model:

$$ y = h \ast x + n $$

Namely the data is a result of a linear combination (Convolution) of $ x $ with Additive Noise.  

If we assume the noise model is Gaussian and our data is also formed by a Gaussian distribution then we should try to minimize (MAP Estimator):

$$ \hat{x} = \frac{1}{2 {\sigma}_{n}^{2}} \left\| h \ast x - y \right\|_{2}^{2} + \frac{1}{ 2 {\sigma}_{x}^{2} } \left\| x \right\|_{2}^{2} $$

Since the convolution is Linear Operation and assuming we in finite dimension then it can be rewritten as:

$$ \hat{x} = \frac{1}{2} \left\| H x - y \right\|_{2}^{2} + \frac{{\sigma}_{n}^{2}}{ 2 {\sigma}_{x}^{2} } \left\| x \right\|_{2}^{2} $$

This is just a Tikhonov Regulaized Least Squares with the solution given by:

$$ \hat{x} = \left( {H}^{T} H + \frac{{\sigma}_{n}^{2}}{{\sigma}_{x}^{2}} I \right)^{-1} {H}^{T} y $$

Now, if we use the discrete form of the circular convolution then $ H $ is a circulant matrix which means all can be solved in the Frequency Domain which results in the same equation as above.