This guide omits a fair amount of detail in order to provide a fairly quick guide.
1) Design in the s-domain
Butterworth filters have their poles evenly spaced on a circlular arc in the left hand plane of the [s/Laplace/analogue]-domain.
For a low-pass filter of order \(N\), the complex valued pole locations are:
\[ p_n = \omega_0 e^\frac{j \pi (2 n + N + 1)}{2 N} \]
Where \(0 <= n < N\) and \(\omega_0\) is the cutoff frequency in radians.
The transfer function in the s domain is:
\[ H(s) = \frac{w_0}{\prod\limits_{n=0}^{N-1}(s - p_n)} \]
Because the poles always occur in conjugate pairs for this kind of filter, it can always be expanded out into a polynomial with real-valued coefficients \(a_n\):
\[ H(s) = \frac{w_0}{\sum\limits_{n=0}^{N-1} a_n s^n} \]
2) Convert the s-domain transfer function to the z-domain
There are a few methods to convert an s-domain transfer function this into a [z/digital/discrete]-domain transfer function. The most common is the bilinear transform which makes the substitution:
\[ s = \alpha \frac{z - 1}{z + 1} \]
Where:
\[ \alpha = \frac{2}{T} \]
Or if you’re performing “frequency warping”:
\[ \alpha = \frac{\omega_0}{\tan(\frac{\omega_0 T}{2})} \]
I won’t go into verbose details here, but you probably want frequency warping. The bilinear transform is a linear approximation to the function \(s=\frac{1}{T}\ln(z)\) which maps the s-domain to the digital z-domain. This linear mapping is analogous to the first term of a series expansion of \(s=\frac{1}{T}\ln(z)\) which is most accurate around z=1 (i.e. DC). The frequency warping operation moves the point where most accuracy is preserved to \(\omega_0\). For higher frequency cutoffs, this becomes important. \(\omega_0\) can be set to the cutoff frequency of the filter to improve the frequency response near the transition.
For a 3-rd order system, the conversion yields (see notes at the end of this page about risks):
\[ \frac{b_3 s^3 + b_2 s^2 + b_1 s^1 + b_0}{a_3 s^3 + a_2 s^2 + a_1 s^1 + a_0} \leftrightarrow \frac{b_3 \alpha^3 (z + 1)^3 + b_2 \alpha^2 (z + 1)^2 (z - 1) + b_1 \alpha (z + 1) (z - 1)^2 + b_0 (z - 1)^3}{a_3 \alpha^3 (z + 1)^3 + a_2 \alpha^2 (z + 1)^2 (z - 1) + a_1 \alpha (z + 1) (z - 1)^2 + a_0 (z - 1)^3} \]
Which can be expanded into a rational polynomial transfer function in \(z\).
3) Convert from the z-domain to a difference equation
First, you must make your z-domain expression causal. Positive powers of \(z\) must be removed otherwise your filter would need to know the future! This is simple, just modify the numerator and denominator of the transfer function by the inverse of the highest power of \(z\). For example:
\[ \frac{Y(z)}{X(z)} = H(z) = \frac{z + 4}{2z^2 - 3z + 2} \]
Multiply the numerator and denominator by \(z^{-2}\).
\[ \frac{Y(z)}{X(z)} = \frac{z^{-1} + 4z^{-2}}{2 - 3z^{-1} + 2z^{-2}} \]
Get \(Y(z)\) as a function of \(X(z)\).
\[ Y(z)(2 - 3z^{-1} + 2z^{-2}) = X(z)(z^{-1} + 4z^{-2}) \]
The difference equation can be found by substituting z-powers as offsets to the data they are using. i.e. for the above case:
\[ 2y[n] = x[n-1] + 4x[n-2] + 3y[n-1] - 2y[n-2] \]
That is what you want to implement in your code.
Things to keep in mind
High-order systems can start behaving very poorly when designed in this way. Even on floating point, once a system reaches 4-th order or higher, if there are poles close to the real axis in the z-domain, the output can become severely degraded. There are ways to mitigate this: the simplest of which is to break the implementation into separate 2nd order filters (biquads). Update: for more information see this blog.