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.

$\endgroup$
  • $\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$ – Alan Wolfe Apr 13 '15 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$ – John Smith Apr 13 '15 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$ – Matt L. Apr 13 '15 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$ – John Smith Apr 13 '15 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$ – Matt L. Apr 13 '15 at 16:56
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.

$\endgroup$
  • $\begingroup$ Matt, how do you put equations in the text? I tried to do fractions and other staff, but did not succeed... $\endgroup$ – John Smith Apr 17 '15 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$$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.