I've seen Julius' MATLAB code and I know what it does.
Essentially, given an LTI filter with impulse response, $h[n]$, and frequency response:
$$\begin{align}
H(e^{j \omega}) &\triangleq \Big| H(e^{j \omega}) \Big| \, e^{j \phi(\omega) } \\
&= \sum\limits_{n=-\infty}^{\infty} h[n] \, e^{-j \omega n}
\end{align}$$
Then $\Big| H(e^{j \omega}) \Big| > 0$ is the magnitude response and $\phi(\omega)$ is the phase response.
To be a minimum-phase filter, the phase response (expressed in radians) is the negative of the Hilbert Transform of the natural log of the magnitude response.
The natural log of the magnitude response, $\log \big| H(e^{j \omega}) \big|$, is just like decibels (dB) but is expressed in a different dimensionless unit, the neper. 8.685889638 dB is equal to 1 neper. Essentially the nepers is the real part of the complex natural log of the complex frequency response and the phase in radians is the imaginary part.
Just like radians is the mathematically natural unit for angle, so also are nepers the mathematically natural unit for relative change of magnitude. A change of magnitude of 1% is nearly the same as 0.01 neper. But, like dB, going up 0.1 neper followed by going down 0.1 neper will land you at exactly the original magnitude. But going up 10% followed by going down 10% will land you at slightly less than the original magnitude.
So, the (unwrapped) minimum phase of a filter, given it's magnitude response is:
$$\begin{align}
\phi(\omega) &= - \mathscr{H} \Big\{ \log \big| H(e^{j \omega}) \big| \Big\} \\
\\
&= - \int\limits_{-\pi}^{\pi} \log \Big| H(e^{j \theta}) \Big| \, \Big( 2 \pi \tan \big(\tfrac{\omega-\theta}{2} \big) \Big)^{-1} \, \mathrm{d}\theta \\
\end{align}$$
Where $\mathscr{H} \big\{ \cdot \big\}$ is the Hilbert Transform.
If we were to evaluate that integral, we would need to do something called "take the Principal Value" to deal with the division by zero when $\theta = \omega$. But we won't do the Hilbert transform that way.
Okay, so Step 2 is understanding that the Discrete-Time Fourier Transform (DTFT) of this sequence:
$$ g[n] \triangleq \begin{cases}
0 \qquad & n<0 \\
1 \qquad & n=0 \\
2 \qquad & n>0 \\
\end{cases} $$
is
$$\begin{align}
G(e^{j\omega}) &= \sum\limits_{n=-\infty}^{\infty} g[n] \, e^{-j \omega n} \\
\\
&= \frac{e^{j\omega}+1}{e^{j\omega}-1} \\
\\
&= \frac{e^{j\omega/2}+e^{-j\omega/2}}{e^{j\omega/2}-e^{-j\omega/2}} \\
\\
&= \frac{2\cos(\omega/2)}{2j\sin(\omega/2)} \\
\\
&= -j \frac{1}{\tan(\omega/2)} \\
\end{align}$$
i'm running outa time. i gotta return to this.