10
$\begingroup$

I've got a first-order Butterworth filter with the cutoff frequency $\omega_c$. Its transfer function is then

$$H(s) = \frac{\omega_c}{s+\omega_c}$$

Using the bilinear transform to find an $H(z)$ (what is that function called?), I get

$$H(z)=\frac{\omega_c}{\frac{2}{T}\frac{z-1}{z+1} + \omega_c} = \frac{\omega_c z + \omega_c}{\left(\frac{2}{T}+\omega_c\right)z + \omega_c-\frac{2}{T}}$$

However, I can't reconcile this result with what Matlab is doing. It seems wrong, no matter what value of $T$. I assume that B and A below are the coefficients of $H(z)$.

>> [B,A] = butter(1,0.5)
B = 0.5000    0.5000
A = 1.0000   -0.0000
>> [B,A] = butter(1,0.6)
B = 0.5792    0.5792
A = 1.0000    0.1584
>> [B,A] = butter(1,0.7)
B = 0.6625    0.6625
A = 1.0000    0.3249
>> [B,A] = butter(1,0.8)
B = 0.7548    0.7548
A = 1.0000    0.5095

What am I misunderstanding?

Improve this question
$\endgroup$
3
  • $\begingroup$ MATLAB does not use analog-to-digital conversion. It designs the filter digitally, therefore the bilinear transform idea may not be applicable. $\endgroup$
    – Phonon
    Nov 7, 2011 at 21:22
  • 1
    $\begingroup$ @Phonon: This answer seems to indicate that Matlab uses the bilinear transform in some way. $\endgroup$
    – Andreas
    Nov 7, 2011 at 21:33
  • $\begingroup$ Late to the game here but the all the uppercase functions H of z/s/\omega are usually called the transfer function. When the argument is time or samples, it's called the impulse response and it's usually lowercased, h. So the transfer function is the transform (Z, Fourier, Laplace depending the application) of the impulse response. $\endgroup$
    – Emanuel Landeholm
    Feb 19, 2015 at 22:13

2 Answers 2

Reset to default
10
$\begingroup$

A couple things:

Before making the substitution $s = \frac{2}{T} \frac{z-1}{z+1}$, you need to prewarp the cutoff frequency by making the substitution:

$$ \omega_{c,w} = \frac{2}{T} \tan (\omega_c\frac{T}{2}) $$

where $\omega_{c,w}$ is the warped cutoff frequency. This is necessary because the bilinear transform maps the left-half plane in the Laplace domain (used in analog filter design) to the unit circle in the $z$-domain in a nonlinear fashion. Therefore, as you approach the Nyquist rate (digital frequencies of $\pm \pi$), the approximation to the analog filter prototype becomes inaccurate.

Also, the second parameter that you are passing to the butter function is the normalized cutoff frequency, not the sampling interval $T$. The normalized frequency used by that function is in the interval $(0,1)$ and is equal to the ratio of the desired cutoff frequency to the Nyquist rate:

$$ \omega_n = \frac{\omega_c}{2 \pi \frac{f_s}{2}} $$

$$ \omega_n = \frac{\omega_c}{\pi f_s} $$

$$ \omega_n = \frac{\omega_c T}{\pi} $$

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you! Now I get the right coefficients. Here, I substituted $\omega_c$ for $omega_c,w$ in the expression for $H(z)$. This worked because I know where the cutoff frequency affects the Butterworth filter. What if I had a general filter and just knew the poles (and zeroes) of $H(s)$? How would I know which values to substitute? $\endgroup$
    – Andreas
    Nov 8, 2011 at 9:31
  • $\begingroup$ Because I suppose the bilinear transform of a rational $H(s)$ could be done knowing only the sample frequency, not the (normalized) cutoff frequency? $\endgroup$
    – Andreas
    Nov 8, 2011 at 9:53
  • $\begingroup$ You can use the bilinear transform to map any $s$-domain system into an approximation in the $z$ domain. Prewarping is not required, but the caveat that the resulting discrete-time system is just an approximation to the analog system applies. Prewarping any frequencies of interest allows you to "stretch" the mapping so that the region of the frequency band that you care most about has as little distortion from the prototype filter as possible. $\endgroup$
    – Jason R
    Nov 8, 2011 at 14:02
5
$\begingroup$

When opening the code for MATLAB's butter function, we see that it uses frequency pre-warping :

%# step 1: get analog, pre-warped frequencies
if ~analog,
    fs = 2;
    u = 2*fs*tan(pi*Wn/fs);
else
    u = Wn;
end
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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