discrete signals - my Butterworth lowpass formulas do not agree with Fisher webpage

I want to implement Butterworth low-pass filter. Thanks to this question, I have found out that the filter coefficients can be generated using Tony Fisher web-site or using his code. But the problem arose when I had tried to verify his formulas myself.

Wikipedia says that the derivation of low-pass formulas is simple: we start with the Butterworth polynomial of order n. I took n=1 $$B_1\left(\frac{s}{\omega_{cutoff}}\right)=1+\frac{s}{\omega_{cutoff}}$$

(note that $\omega_{cutoff}$ is angular frequency, not usual one) then do bilinear transform $$s = 2f_{sampling}\cdot\frac{z-1}{z+1}$$

($f_{sampling}$ is usual frequency, not angular) and rewrite the resulting fraction in the form with non-positive powers of $z$.

To make the story short, my final formula for transfer function is $$H(z) = \frac{Y(z)}{X(z)}$$ where $$Y(z) = \frac{a}{1+a}+\frac{a}{1+a}z^{-1}$$ $$X(z) = 1-\frac{1-a}{1+a}z^{-1}$$ and $$a = \frac{\omega_{cutoff}}{2f_{sampling}}$$ and the resulting formula is $$y[n] = \frac{a}{1+a}x[n]+\frac{a}{1+a}x[n-1]+\frac{1-a}{1+a}y[n-1]$$ For testing I use $f_{cutoff}=1$ (which is $\omega_{cutoff}=2\pi$) and $f_{sampling}=30$.

We should not worry about coefficients in front of $x[]$ since Fisher multiplies them by gain factor anyway, but the coefficient in front of $y[n-1]$ equals

$$\frac{1-a}{1+a} = 0.8104139027$$

while Fisher's web-page for $f_{cutoff}=1$ and $f_{sampling}=30$ gives $0.8097840332$.

If you had a patience to finish this reading, may be you can explain, where I (or Fisher?) am wrong.

Answer

The problem is in the way you apply the bilinear transform. You have to use the appropriate (pre-)warping of the frequencies. Since the bilinear transform warps the frequency axis, you have to make sure that the corner frequency of the discrete-time filter is correct. One way to do that is as follows. The bilinear transform is defined as

$$s=k\frac{z-1}{z+1}\tag{1}$$

with some constant $k$ yet to be defined. If we denote the analog frequency by $\Omega$, and the normalized discrete-time frequency by $\omega$, Eq. (1) becomes for $s=j\Omega$ and $z=e^{j\omega}$

$$j\Omega=k\frac{e^{j\omega}-1}{e^{j\omega}+1}=k\frac{e^{j\omega/2}-e^{-j\omega/2}}{e^{j\omega/2}+e^{-j\omega/2}}=jk\tan(\omega/2)\tag{2}$$

Eq. (2) describes the frequency warping caused by the bilinear transform. If we use an analog lowpass filter with a normalized corner frequency $\Omega_0=1$, we must choose the constant $k$ such that for the desired discrete-time corner frequency $\omega_0$ the term on the right-hand side of Eq. (2) becomes $1$:

$$k=\frac{1}{\tan(\omega_0/2)}\tag{3}$$

where $\omega_0$ is the desired corner frequency of the digital filter. Eq. (3) and Eq. (1) define the appropriately normalized bilinear transform that you must use.

So for your example, the normalized first-order analog Butterworth lowpass transfer function is given by

$$H(s)=\frac{1}{1+s}\tag{4}$$

Applying the bilinear transform gives

$$H(z)=\frac{1}{1+k\frac{z-1}{z+1}}=\frac{z+1}{z(1+k)+1-k}=\frac{1}{1+k}\cdot\frac{1+z^{-1}}{1+\frac{1-k}{1+k}z^{-1}}\tag{5}$$

Ignoring the gain term $1/(1+k)$, the corresponding difference equation is

$$y[n]=x[n]+x[n-1]-\frac{1-k}{1+k}y[n-1]\tag{6}$$

With a desired corner frequency $\omega_0=2\pi/30$ you get from (3) $k=9.5144$ and $(1-k)/(1+k)=-0.80978$, just like on Fisher's website.