# Design IIR Bandpass Filters

In this post, I present a method to design Butterworth IIR bandpass filters. My previous post [1] covered lowpass IIR filter design, and provided a Matlab function to design them. Here, we’ll do the same thing for IIR bandpass filters, with a Matlab function bp_synth.m. Here is an example function call for a bandpass filter based on a 3^{rd} order lowpass prototype:

N= 3; % order of prototype LPF fcenter= 22.5; % Hz center frequency, Hz bw= 5; % Hz -3 dB bandwidth, Hz fs= 100; % Hz sample frequency [b,a]= bp_synth(N,fcenter,bw,fs) b = 0.0029 0 -0.0087 0 0.0087 0 -0.0029 a = 1.0000 -0.8512 2.6169 -1.3864 2.1258 -0.5584 0.5321

There are seven “b” (numerator) and seven “a” (denominator) coefficients, so H(z) is 6^{th} order. To find the frequency response:

[h,f]= freqz(b,a,512,fs); H= 20*log10(abs(h));

The magnitude response of the filter is shown in Figure 1, along with the response of the equivalent analog bandpass filter.

This article is available in PDF format for easy printing.

Figure 1. Magnitude response of bandpass filter based on N= 3 lowpass prototype.

fcenter= 22.5 Hz, bw= 5 Hz, and fs= 100 Hz. The equivalent analog filter’s response is also shown.

## Filter Synthesis

Here is a summary of the steps for computing the bandpass filter coefficients. Recall from the previous post that F is continuous (analog) frequency in Hz and Ω is continuous radian frequency.

- Find the poles of a lowpass analog prototype filter with Ω
_{c}= 1 rad/s. - Given upper and lower -3 dB frequencies of the digital bandpass filter, find the corresponding frequencies of the analog bandpass filter (pre-warping).
- Transform the analog lowpass poles to analog bandpass poles.
- Transform the poles from the s-plane to the z-plane, using the bilinear transform.
- Add N zeros at z= -1 and N zeros at z= +1, where N is the order of the lowpass prototype.
- Convert poles and zeros to polynomials with coefficients a
_{n}and b_{n}.

These steps are similar to the lowpass design procedure, with steps 2, 3, 5, and 6 modified for the bandpass case. Now let’s look at the design procedure in detail. A Matlab function bp_synth that performs the filter synthesis is provided in Appendix A.

1. Poles of the analog lowpass prototype filter. For a Butterworth filter of order N with Ω_{c} = 1 rad/s, the poles are given by [2,3]:

$$p'_{ak}= -sin(\theta)+jcos(\theta)$$

where $$\theta=\frac{(2k-1)\pi}{2N},\quad k=1:N$$

Here we use a prime superscript on p to distinguish the lowpass prototype poles from the yet to be calculated bandpass poles.

2. Given upper and lower -3 dB frequencies of the digital bandpass filter, find the corresponding frequencies of the analog bandpass filter. We define f_{center} as the center frequency in Hz, and bw_{Hz} as the -3 dB bandwidth in Hz. Then the -3 dB discrete frequencies are :

$$f_1=f_{center}-bw_{Hz}/2$$

$$f_2=f_{center}+bw_{Hz}/2$$

As before, we’ll adjust (pre-warp) the analog frequencies to
take the nonlinearity of the bilinear transform into account:

$$F_1=\frac{f_s}{\pi}tan\left(\frac{\pi
f_1}{f_s}\right)$$

$$F_2=\frac{f_s}{\pi}tan\left(\frac{\pi f_2}{f_s}\right)$$

We need to define two more quantities:

BW_{Hz} = F_{2} – F_{1 }Hz is the pre-warped -3 dB bandwidth, and

$F_0=\sqrt{F_1F_1}$ is the geometric mean of F1 and F2.

3. Transform the analog lowpass poles to analog bandpass poles. See Appendix B for a derivation of the transformation. For each lowpass pole p_{a}' , we get two complex-conjugate bandpass poles:

$$p_a=2\pi F_0\left[\frac{BW_{Hz}}{2F_0}p'_a \pm j\sqrt{1-\left(\frac{BW_{Hz}}{2F_0}p'_a\right)^2}\ \right]$$

4. Transform the
poles from the s-plane to the z-plane, using the bilinear transform. This is the same as for the IIR lowpass,
except there are 2N instead of N poles:

$$p_k=\frac{1+p_{ak}/(2f_s)}{1-p_{ak}/(2f_s)},\quad
k=1:2N$$

5. Add N zeros at z=
-1 and N zeros at z= +1. See Appendix B
for details. We can now write H(z) as:

$$H(z)=K\frac{(z+1)^N(z-1)^N}{(z-p_1)(z-p_2)...(z-p_{2N})}\qquad(1)$$

In bp_synth, we represent the N zeros at -1 and N zeros at +1 as a vector:

q= [-ones(1,N) ones(1,N)]

6. Convert poles and zeros to polynomials with coefficients a_{n} and b_{n}. If we expand the numerator and denominator of equation 1 and divide numerator and denominator by z^{2N}, we get polynomials in z^{-n}:

$$H(z)=K\frac{b_0+b_1z^{-1}+...+b_{2N}z^{-2N}}{1+a_1z^{-1}+...+a_{2N}z^{-2N}}$$

The Matlab code to perform the expansion is:

a= poly(p) a= real(a) b= poly(q)

We want H(z) to have a gain of 1 at f = f_{0}. We can do this by evaluating H(z) at f_{0} and setting K = 1/|H(f_{0})|. The Matlab code is:

f0= sqrt(f1*f2); h= freqz(b,a,[f0 f0],fs); K= 1/abs(h(1));

## Example

Here’s an example IIR bandpass filter based on a 2^{nd} order lowpass. We’ll show the different poles and zeros computed by bp_synth.m. (Note the pole and zero values are not printed by the code listed in Appendix A). The function call is as follows:

N= 2; % order of prototype LPF fcenter= 20; % Hz center frequency, Hz bw= 4; % Hz -3 dB bandwidth, Hz fs= 100; % Hz sample frequency [b,a]= bp_synth(N,fcenter,bw,fs)

Here are the pole and zero computed values. Note Matlab uses an a+bi format for complex numbers.

Two poles of analog lowpass butterworth prototype with Ω’ = 1 rad/s

(p_lp is the same as p_{a}’):

p_lp = -0.7071 + 0.7071i -0.7071 - 0.7071i

Four poles of analog bandpass:

pa = 1.0e+002 * -0.1490+ 1.5854i -0.1234+ 1.3130i -0.1234- 1.3130i -0.1490- 1.5854i

Four poles of IIR bandpass:

p= 0.2053+ 0.8892i 0.3627+ 0.8426i 0.3627- 0.8426i 0.2053- 0.8892i

Two zeros at z= -1 and two zeros at z= 1:

q = -1 -1 1 1

The IIR bandpass poles and zeros are shown in Figure 2. Now let’s look at the filter coefficients.

Numerator coefficients before scaling:

b = 1 0 -2 0 1

Numerator scale factor:

K = 0.0134

Filter coefficients b and a:

b = 0.0134 0 -0.0267 0 0.0134 a = 1.0000 -1.1361 1.9723 -0.9498 0.7009

From K, b, and a we can write the filter’s transfer function:

$$H(z)=.0134*\frac{1-2z^{-2}+z^{-4}}{1-1.1361z^{-1}+1.9723z^{-2}-.9498z^{-3}+.7009z^{-4}}$$

The filter’s magnitude response is plotted as the
N=2 curve in Figure 3.

Figure 2. Pole-zero plot of bandpass filter based on 2^{nd} order lowpass prototype.

fcenter= 20 Hz, bw= 4 Hz, and fs= 100 Hz. Zeros at z= +/-1 are 2^{nd} order.

## Filter Performance

As was the case for lowpass IIR filters, the bandpass IIR response has more rejection than that of the analog prototype as you approach f_{s}/2 (see Figure 1). This is due to the bilinear transform, which compresses the entire bandwidth into 0 to f_{s}/2. Thus, the analog zeros at infinity occur at f_{s}/2 for the discrete-time filter.

Figure 3 compares the magnitude response of bandpass filters based on 1^{st}, 2^{nd}, and 3^{rd} order Butterworth lowpass filters. They all have the same -3 dB bandwidth of 4 Hz. Figure 4 shows the group delay response of the same filters: as you would expect for IIR filters, the group delay is not flat.

## Caveats

Performance of IIR filters can be hurt by coefficient quantization and roundoff noise [4,5]. In particular, when the poles are close to the unit circle, quantization of the denominator coefficients can cause a pole to move outside the unit circle, making the filter unstable. Bandpass filters with narrow bandwidth and higher order are most susceptible.

Figure 3. Magnitude response of bandpass filters vs. order of lowpass prototype.

fcenter= 20 Hz, bw= 4 Hz, and fs= 100 Hz.

Figure 4. Group Delay of bandpass filters vs. order of lowpass prototype.

fcenter= 20 Hz, bw= 4 Hz, and fs= 100 Hz. [gd,f]= grpdelay(b,a,512,fs);

## Appendix A. Matlab Function bp_synth

This program is provided as-is without any guarantees or warranty. The author is not responsible for any damage or losses of any kind caused by the use or misuse of the program.

%function [b,a]= bp_synth(N,fcenter,bw,fs) 12/29/17 neil robertson % Synthesize IIR Butterworth Bandpass Filters % % N= order of prototype LPF % fcenter= center frequency, Hz % bw= -3 dB bandwidth, Hz % fs= sample frequency, Hz % [b,a]= vector of filter coefficients % function [b,a]= bp_synth(N,fcenter,bw,fs) f1= fcenter- bw/2; % Hz lower -3 dB frequency f2= fcenter+ bw/2; % Hz upper -3 dB frequency if f2>=fs/2; error('fcenter+ bw/2 must be less than fs/2') end if f1<=0 error('fcenter- bw/2 must be greater than 0') end % find poles of butterworth lpf with Wc = 1 rad/s k= 1:N; theta= (2*k -1)*pi/(2*N); p_lp= -sin(theta) + j*cos(theta); % pre-warp f0, f1, and f2 F1= fs/pi * tan(pi*f1/fs); F2= fs/pi * tan(pi*f2/fs); BW_hz= F2-F1; % Hz prewarped -3 dB bandwidth F0= sqrt(F1*F2); % Hz geometric mean frequency % transform poles for bpf centered at W0 % note: alpha and beta are vectors of length N; pa is a vector of length 2N alpha= BW_hz/F0 * 1/2*p_lp; beta= sqrt(1- (BW_hz/F0*p_lp/2).^2); pa= 2*pi*F0*[alpha+j*beta alpha-j*beta]; % find poles and zeros of digital filter p= (1 + pa/(2*fs))./(1 - pa/(2*fs)); % bilinear transform q= [-ones(1,N) ones(1,N)]; % N zeros at z= -1 (f= fs/2) and z= 1 (f = 0) % convert poles and zeros to numerator and denominator polynomials a= poly(p); a= real(a); b= poly(q); % scale coeffs so that amplitude is 1.0 at f0 f0= sqrt(f1*f2); h= freqz(b,a,[f0 f0],fs); K= 1/abs(h(1)); b= K*b;

## Appendix B. Lowpass to Bandpass Frequency Transformation

In the s-domain, we want to transform a normalized lowpass filter with -3 dB frequency of 1 rad/s to a bandpass filter with a given bandwidth and center frequency [6,7]. Once we have the bandpass H(s), we can use the bilinear transform to obtain H(z). Note it is also possible to perform the lowpass to bandpass transformation in the Z-domain [8].

We define two frequency variables:

Ω = Im(s) imaginary part of the complex frequency

Ω’ = Im(s’) imaginary part of the normalized complex frequency

(As in the previous post, we use Ω for continuous frequency and reserve ω for discrete frequency). To transform H(s’) into H(s), substitute the following for s’:

$$s'=\frac{\Omega_0}{BW}\left(\frac{s}{\Omega_0}+\frac{\Omega_0}{s}\right)\qquad(B-1)$$

where BW = Ω_{2}-Ω_{1} = -3 dB bandwidth and $\Omega_0=\sqrt{\Omega_1\Omega_2}$. If we let s= jΩ, we then have:

$$\Omega'=\frac{\Omega_0}{BW}\left(\frac{\Omega}{\Omega_0}-\frac{\Omega_0}{\Omega}\right)$$

Then, for example, if Ω’= 0, we see that Ω = Ω_{0}. That is, 0 rad/s transforms to Ω_{0} rad/s.

Rearranging B-1, we have:

$$s'=\frac{1}{BW}\left(\frac{s^2+\Omega_0^2}{s}\right)$$

Then, solving for s:

$$s=\Omega_0\left[\frac{BW}{2\Omega_0}s' \pm j\sqrt{1-\left(\frac{BW}{2\Omega_0}s'\right)^2}\ \right]$$

If we replace s’ with a normalized lowpass pole location p_{a}’, we have:

$$p_a=\Omega_0\left[\frac{BW}{2\Omega_0}p'_a \pm j\sqrt{1-\left(\frac{BW}{2\Omega_0}p'_a\right)^2}\ \right]\qquad(B-2)$$

For each lowpass pole p_{a}’ , we get two complex-conjugate bandpass poles. Thus if the lowpass filter has order N, H(s) for the bandpass filter has order 2N. The IIR filter’s poles are computed from p_{a} using the bilinear transform, as discussed in step 4 in the section on filter synthesis.

For the Matlab function, we replace Ω_{0} with 2πF_{0} and replace BW with 2π*BW_{Hz}:

$$p_a=2\pi F_0\left[\frac{BW_{Hz}}{2F_0}p'_a \pm j\sqrt{1-\left(\frac{BW_{Hz}}{2F_0}p'_a\right)^2}\ \right]$$

Figures B.1 and B.2 plot the s-plane poles computed
by equation B-2 for lowpass order N= 2 and N=4, respectively. The zeros at s= 0 will be discussed next.

Figure B.1 s-plane transformation of N= 2 lowpass poles to bandpass poles and zeros.

Left: Normalized lowpass Complex-conjugate poles, N= 2

Right: Bandpass poles and zero for Ω_{0}= 10 rad/s and BW= 2 rad/s

The zero at s= 0 is 2^{nd} order.

Figure B.2 s-plane transformation of N= 4 lowpass poles to bandpass poles and zeros.

Left: Normalized lowpass Complex-conjugate poles, N= 4

Right: Bandpass poles and zeros for Ω_{0}= 10 rad/s and BW= 2 rad/s

The zero at s= 0 is 4^{th} order.

## Zeros of H(s) and H(z)

Consider a single-pole normalized lowpass transfer function in s’:

$$H(s')=\frac{1}{s'-\alpha}$$

It has a zero at s’= infinity. From B-1, we have for the corresponding
bandpass filter:

$$H(s)=\frac{1} {\frac{\Omega_0}{BW}\left(\frac{s}{\Omega_0}+\frac{\Omega_0}{s}\right)-\alpha}$$

H(s) has a zero at s= 0 and a zero at infinity. In general, if H(s’) has N zeros at infinity, then H(s) has N zeros at 0 and N zeros at infinity.

In my last post [1], I showed that the bilinear transform maps a zero at infinity to one at f_{s}/2. Thus z = exp(j2πf/f_{s}) = exp(jπ) = -1. For a zero at s= 0, z= exp(j0)= 1. We can now write H(z) in pole-zero format as:

$$H(z)=K\frac{(z+1)^N(z-1)^N}{(z-p_1)(z-p_2)...(z-p_{2N})}$$

The numerator of H(z) can be written as $(z^2-1)^N$. If you expand this, you get the numerator polynomial coefficients as follows:

N= 1: b= [1 0 -1] N= 2: b= [1 0 -2 0 1] N= 3: b= [1 0 -3 0 3 0 -1]

Thus odd-power numerator coefficients are always zero and all numerator coefficients are integers before scaling by K.

Figure B.3 shows the z-plane pole-zero plot of a bandpass filter based on a 4^{th} order Butterworth lowpass. You can compare it to the pole-zero plot of the continuous-time bandpass filter in Figure B.2.

Figure B.3. Z-plane pole-zero plot of bandpass filter based on 4^{th} order lowpass prototype.

fcenter= 20 Hz, bw= 5 Hz, and fs= 100 Hz. Zeros at z= +/-1 are 4^{th} order.

## Appendix C Notation

In this article, we added to the notation used in the post on lowpass filter synthesis [1].

Notation for normalized lowpass continuous frequencies:

normalized radian frequency Ω’ rad/s

normalized complex frequency s’= σ + jΩ’

Notation for bandpass filter continuous frequencies:

-3 dB frequencies Ω_{1}, Ω_{2} rad/s

-3 dB bandwidth BW rad/s

geometric mean frequency $\Omega_0=\sqrt{\Omega_1\Omega_2}$ rad/s

- 3 dB bandwidth, Hz BW_{Hz} Hz

geometric mean frequency, Hz F_{0} Hz

Notation for bandpass filter discrete frequencies:

-3 dB frequencies, Hz f_{1}, f_{2} Hz

geometric mean frequency, Hz f_{0} Hz

The notation already developed in the previous post:

continuous frequency F Hz

continuous radian frequency Ω radians/s

complex frequency s = σ + jΩ

discrete frequency f Hz

discrete normalized radian frequency ω = 2πf/f_{s}radians, where f_{s} = sample frequency

Finally, the Matlab function represents -3 dB discrete bandwidth (Hz) as bw and the -3 dB continuous bandwidth (Hz) as BW_hz.

## References

1. Robertson, Neil , “Design IIR Butterworth Filters Using 12 Lines of Code”, Dec 2017 https://www.dsprelated.com/showarticle/1119.php

2. Williams, Arthur B. and Taylor, Fred J., __Electronic Filter Design Handbook__, 3^{rd} Ed., McGraw-Hill, 1995, section 2.3

3. Analog Devices Mini Tutorial MT-224, 2012 http://www.analog.com/media/en/training-seminars/tutorials/MT-224.pdf

4. Oppenheim, Alan V. and Shafer, Ronald W., __Discrete-Time Signal Processing__, Prentice Hall, 1989, sections 6.8 and 6.9.

5. Lyons, Richard G., __Understanding Digital Signal Processing__, 2^{nd} Ed., Pearson, 2004, section 6.7

6. Blinchikoff, Herman J., and Zverev,Anatol I., Filtering in the Time and Frequency Domains, Wiley, 1976, section 4.4.

7. Nagendra Krishnapura , “E4215: Analog Filter Synthesis and Design Frequency Transformation”, 4 Mar. 2003 http://www.ee.iitm.ac.in/~nagendra/E4215/2003/handouts/freq_transformation.pdf

8. Oppenheim, Alan V. and Shafer, Ronald W., __Discrete-Time Signal Processing__, Prentice Hall, 1989, section 7.2

Neil Robertson January, 2018

**Previous post by Neil Robertson:**

Design IIR Butterworth Filters Using 12 Lines of Code

**Next post by Neil Robertson:**

Design IIR Band-Reject Filters

## Comments:

- Comments
- Write a Comment Select to add a comment

Thanks Neil,

Another buried secret of Matlab unearthed. I tested N = 4 and got same as matlab butter

[b,a]= butter(4,2*[.1-.1/2,.1+.1/2]);

[b,a]= bp_synth(4,.1,.1,1);

Kaz

Kaz,

Thanks, I had not tested it vs the Matlab function. Nice to know it matches.

Neil

Hi,

Nice to see the 'behind the scenes' equations! Thanks!

This could be interesting to other readers: another way to implement bandpass filter is to use the 'sum-of-allpass' structure, which provides orthogonal control over the 'Q' and 'f0' from each coefficient.

H = 0.5*(1-A(z)), where A(z) is a 2nd order allpass filter. By choosing the lattice structure for A(z), you get orthogonal control over the filter specifications.

For example:

%some settings fs = 48e3; %some bandpass filter specs f0 = 10e3; bw = 1e3; %design using butterworth [b,a] = butter(1,[f0-bw f0+bw]/(fs/2),'bandpass'); Hbutter = filt(b,a); % 0.1163 - 0.1163 z^-2 % Hbutter = ----------------------------- % 1 - 0.4614 z^-1 + 0.7673 z^-2 %design using sum-of-allpass k1 = -cos(2*pi*f0/fs); k2 = (1-tan(2*pi*bw/fs)) / (1+tan(2*pi*bw/fs)); Az = filt([k2 k1*(1+k2) 1],[1 k1*(1+k2) k2]); Hbp = (1/2)*(1-Az); % 0.1163 - 0.1163 z^-2 % Hbp = ---------------------------- % 1 - 0.4574 z^-1 + 0.7673 z^-2 %Show they are the same fvtool(Hbp.num{1}, Hbp.den{1}, b,a,'Fs', fs);

You can see that k1 only depends on f0, and k2 only depends on the bandwidth. This allows for easy tuning. For example, if you keep k1 constant and sweep k2 and vice versa, you get those responses:

Those sum-of-allpass structures are also extremely robust to coefficient and data quantization, less prone to limit cycles and have many other nice properties.

More details, examples and cool things that one can do with allpass filters can be found here:

The Digital All-Pass Filter: A Versatile Signal Processing Building Block

Thanks -- looks interesting.

Neil

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.