0
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • $\begingroup$ the difference is small enough that it could be due to a rounding error, caused by you doing a different order of operations than they are. Or your programs could be using different floating point settings (fast math vs precise for instance or disabling denormals). IMO anyways! $\endgroup$ Commented Apr 13, 2015 at 5:42
  • $\begingroup$ I believe it is not a floating point settings issue. if I set f_sampling = 10 Hz and f_cutoff =1 Hz the difference between my formula and Tony Fisher becomes much larger. $\endgroup$ Commented Apr 13, 2015 at 6:35
  • $\begingroup$ If you want to design a discrete-time filter with a certain cut-off frequency, the analog prototype must have a different (pre-warped) cut-off frequency. Did you take this into account? $\endgroup$ Commented Apr 13, 2015 at 7:18
  • $\begingroup$ Matt, I did not know that.. as I wrote, I just more or less followed the recipes, given in wiki - the article called Butterworth filter and the one called Bilinear transformation. Can you tell me please which formula shoul I use in order to convert continuous frequency into discrete one? Or, since you are profi in this field, can you recommend may be a text book on DSP, where I can find it? (by the way, any chance that you know a book which explains how to do high-pass and band-pass filters? wiki explains how to calculate transfer functions only for low-pass Butterworth filters) $\endgroup$ Commented Apr 13, 2015 at 14:03
  • $\begingroup$ The book by Orfanidis is quite good at these things. Find the link in this answer. As soon as I have more time I'll write up an answer with the details of the design process. $\endgroup$ Commented Apr 13, 2015 at 16:56

2 Answers 2

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Matt, how do you put equations in the text? I tried to do fractions and other staff, but did not succeed... $\endgroup$ Commented Apr 17, 2015 at 2:13
1
$\begingroup$

Just wanted to supplement an excellent answer by Matt L., since it was not very clear to me how he calculated the numerical value of $k$ in equation $(3)$. After reading the book "Introduction to signal processing" by Orfandis I have found the formula $$k = \frac{1}{\tan\left(\frac{\omega_c}{2}\right)}$$ where $\omega_c$ is the so called digital cutoff frequency $$\omega_c = \frac{2\pi f_c}{f_s}$$ where $f_c$ is the desired cutoff frequency of the filter and $f_s$ is the sampling frequency. Thus we arrive at the resulting expression $$k = \frac{1}{\tan\left(\frac{\pi f_c}{f_s}\right)}$$ In my question I used sampling frequency $f_s = 30\ Hz$ and the desired cutoff frequency $f_c = 1\ Hz$, which gave $$k = \frac{1}{\tan\left(\pi\frac{1Hz}{30Hz}\right)} = 9.5143646$$ and $$-\frac{1-k}{1+k} = 0.80978403$$

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.