Design IIR BandReject Filters
In this post, I show how to design IIR Butterworth bandreject filters, and provide two Matlab functions for bandreject filter synthesis. Earlier posts covered IIR Butterworth lowpass [1] and bandpass [2] filters. Here, the function br_synth1.m designs bandreject filters based on null frequency and upper 3 dB frequency, while br_synth2.m designs them based on lower and upper 3 dB frequencies. I’ll discuss the differences between the two approaches later in this article. Here is an example function call to br_synth1.m:
N= 2; % order of prototype LPF f0= 15; % Hz Null frequency fU= 16; % Hz fU= upper 3 dB freq fs= 100; % Hz sample frequency [b,a]= br_synth1(N,f0,fU,fs) b = 0.9167 2.1554 3.1004 2.1554 0.9167 a = 1.0000 2.2492 3.0935 2.0616 0.8404
There are five “b” (numerator) and five “a” (denominator) coefficients, so H(z) is 4th order. The frequency response is computed as follows, and is plotted in Figure 1.
[h,f]= freqz(b,a,4096,fs); H= 20*log10(abs(h));
This article is available in PDF format for easy printing.
Figure 1. Magnitude response of bandreject filter based on N= 2 lowpass prototype.
f0= 15 Hz, fU= 16 Hz, and fs= 100 Hz.
Filter Synthesis
Here is a summary of the steps for computing the bandreject filter coefficients. Note 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 the null frequency f_{0} and upper 3 dB frequency f_{U} of the digital bandreject filter, find the corresponding frequencies of the analog bandreject filter (prewarping).
 Transform the analog lowpass poles to analog bandreject poles.
 Transform the poles from the splane to the zplane, using the bilinear transform.
 Add N zeros at z= exp(jω_{0}) and N zeros at z= z= exp(jω_{0}), 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 bandpass design procedure, with steps 2, 3, 5, and 6 modified for the bandreject case. Now let’s look at the design procedure in detail. The Matlab function br_synth1.m 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 [3,4]:
$$p'_{ak}= sin(\theta)+jcos(\theta)$$
where $$\theta=\frac{(2k1)\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 bandreject poles.
2. Given null frequency f_{0} and upper 3 dB frequency f_{U} of the digital bandreject filter, find the corresponding frequencies of the analog bandreject filter. As before, we’ll adjust (prewarp) the analog frequencies to take the nonlinearity of the bilinear transform into account:
$$F_0=\frac{f_s}{\pi}tan\left(\frac{\pi f_0}{f_s}\right)$$
$$F_U=\frac{f_s}{\pi}tan\left(\frac{\pi f_U}{f_s}\right)$$
We also need to find the lower 3 dB frequency F_{L} and the bandwidth BW_{Hz}:
$F_L=F_0^2/F_U$
$BW_{Hz}=F_UF_L$_{}
3. Transform the analog lowpass poles to analog bandreject poles. See Appendix B for a derivation of the transformation. For each lowpass pole p_{a}' , we get two complexconjugate bandreject poles:
$$p_a=2\pi F_0\left[\frac{BW_{Hz}}{2F_0p'_a} \pm
j\sqrt{1\left(\frac{BW_{Hz}}{2F_0p'_a}\right)^2}\ \right]$$
4. Transform the poles from the splane to the zplane, using the bilinear transform (Note there are 2N poles). This is the same as for the IIR bandpass:
$$p_k=\frac{1+p_{ak}/(2f_s)}{1p_{ak}/(2f_s)},\quad k=1:2N$$
5. Add N zeros on the unit circle at z= exp(jω_{0}) and N zeros at z= exp(jω_{0}), where ω_{0}= 2πf_{0}/f_{s}. See Appendix B for details. We can now write H(z) as:
$$H(z)=K\frac{(z e^{j\omega_0})^N(ze^{j\omega_0})^N}{(zp_1)(zp_2)...(zp_{2N})}\qquad(1)$$
In br_synth1.m, we represent the N zeros at exp(jω_{0}) and N zeros at exp(jω_{0}) as a vector:
q= [exp(j*w0)*ones(1,N) exp(j*w0)*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}}\qquad(2)$$
The Matlab code to perform the expansion is:
a= poly(p) a= real(a) b= poly(q) b= real(b)
We want H(z) to have a gain of 1 at ω= 0. Letting z= e^{jω}, we have z= 1. Then, referring to equation 2, we have gain at ω= 0 of:
$$H(z=1)=K\frac{\sum b}{\sum a}$$
So, for gain of 1 at ω= 0, we make $K=\sum a/\sum b$.
Example 1
Let’s repeat the example of Figure 1, adding some more detail. We’ll plot the poles and zeros computed by br_synth1.m, and we’ll also look at the group delay. (Note the pole and zero values are not printed by the code listed in Appendix A). Repeating the function call from above:
N= 2; % order of prototype LPF f0= 15; % Hz Null frequency fU= 16; % Hz fU= upper 3 dB freq fs= 100; % Hz sample frequency [b,a]= br_synth1(N,f0,fU,fs)
Here are the poles and zeros:
p = 0.5968 + 0.7504i 0.5278 + 0.7973i 0.5278  0.7973i 0.5968  0.7504i q = 0.5878 + 0.8090i 0.5878 + 0.8090i 0.5878  0.8090i 0.5878  0.8090i
The zeros are on the unit circle, with N= 2 zeros at exp(jω0) and 2 zeros at exp(jω0). The poles and zeros are plotted in Figure 2. Now let’s look at the filter coefficients.
Numerator coefficients before scaling (note symmetry with respect to center):
b = 1.0000 2.3511 3.3820 2.3511 1.0000
Numerator scale factor:
K = 0.9167
Filter coefficients b and a:
b = 0.9167 2.1554 3.1004 2.1554 0.9167 a = 1.0000 2.2492 3.0935 2.0616 0.8404
From K, b, and a, we can write the filter’s transfer function:
$$H(z)=.9167*\frac{12.3511z^{1}+3.382z^{2}2.3511z^{3}+z^{4}}{12.2492z^{1}+3.0935z^{2}2.0616z^{3}+.8404z^{4}}$$
The magnitude response was plotted in Figure
1. Figure 3 shows the magnitude and group
delay responses over a 10 Hz bandwidth centered at the null frequency. The group delay response is not very flat,
even outside the 1 dB frequency of the amplitude response.
Figure 2. Poles and zeros of bandreject filter based on 2nd order lowpass prototype.
f0= 15 Hz, fU= 16 Hz, and fs= 100 Hz. Zeros at z= exp(+jω_{0}) are 2^{nd} order.
Figure 3. Response of bandreject filter based on 2nd order lowpass prototype.
f0= 15 Hz, fU= 16 Hz, and fs= 100 Hz.
Top: Magnitude Response Bottom: Group Delay [gd,f]= grpdelay(b,a,4096,fs);
Example 2
Let’s look at the magnitude response for different filter orders, for a null frequency f_{0}= 30 MHz and 3 dB bandwidth bw of about 4 Hz. For a narrowband filter, the 3 dB frequencies f_{U} and f_{L} are approximately f_{0} + bw/2 and f_{0 } bw/2, respectively. So we have f_{U}= 30 + 2 = 32 Hz. The function call is:
[b,a]= br_synth1(N,f0,fU,fs)
Letting N = 1, 2, and 3, we compute the magnitude response and obtain the plots shown in Figure 4. All of the plots have f_{U} of exactly 32 Hz and f_{L} of slightly less than 28 Hz.
Figure 4. Magnitude response of bandreject filters vs. order of lowpass prototype.
f0= 30 Hz, fU= 32 Hz, and fs= 100 Hz.
Comparing the two Matlab functions for Bandreject Filter Synthesis
From Appendix B, the null frequency of the analog bandreject filter is:
$$\Omega_0=\sqrt{\Omega_L\Omega_U}\qquad(3)$$
where Ω_{L} and Ω_{H} are the lower and upper 3 dB frequencies. In designing a filter, we can choose two of the three parameters in Equation 3, and the other is then determined. For br_synth1.m, we chose Ω_{0} and Ω_{H}, while for br_synth2.m, we chose Ω_{L} and Ω_{H}. The disadvantage of the latter choice is that the null frequency falls at the geometric mean of Ω_{L} and Ω_{H}, which is not convenient if we want a deep null at a particular frequency.
Figure 5 compares the responses for the two functions. The top plot is for br_synth1.m, with f0= 30 and fU= 34. The bottom plot is for br_synth2.m, with fL= 26 and fU= 34, which results in a null frequency of 30.168 Hz. Table 1 shows how each Matlab function computes frequencies.
Table 1. Frequency Computations of br_synth1.m and br_synth2.m
br_synth1(N,f0,fU,fs) 
br_synth2(N,fL,fU,fs) 

discrete input frequencies 
f0, fU 
fL, fU 
Analog prewarped frequencies 
F0= fs/pi * tan(pi*f0/fs) FU= fs/pi * tan(pi*fU/fs) FL= F0^2/FU

FL= fs/pi * tan(pi*fL/fs) FU= fs/pi * tan(pi*fU/fs) F0= sqrt(FL*FU); 
resulting discrete frequencies 
fL= fs/pi*atan(pi*FL/fs) 
f0= fs/pi*atan(pi*F0/fs) 
Note that [b,a]= br_synth2(N,fL,fU,fs) gives the same results as the Matlab function
[b,a]= butter(N,[fL fU]*2/fs,’stop’).
Figure 5. Magnitude Responses of two Matlab functions, fs= 100 Hz
Top: [b,a]= br_synth1(2,30,34,fs)
Bottom: [b,a]= br_synth2(2,26,34,fs)
Coefficient Quantization and Null Depth
The frequency response of a digital filter is H(e^{jω}), with 0 < ω < 2π. In other words, we evaluate H(z) on the unit circle. As shown in Figure 2, a bandreject filter has N zeros at exp(jω_{0}) and N zeros at exp(jω_{0}), which produce a null in the response at ω_{0}. For a realworld filter, the numerator coefficients b are quantized, which means the zeros don’t fall exactly at the desired spots on the unit circle. As a result, a practical bandreject filter does not have a perfect null at ω_{0}.
Let’s look at how quantization of the numerator coefficients affects the response of the filter in Example 1, which is based on an N= 2 lowpass prototype. Repeating the function call once again:
N= 2; % order of prototype LPF f0= 15; % Hz Null frequency fU= 16; % Hz fU= upper 3 dB freq fs= 100; % Hz sample frequency [b,a]= br_synth1(N,f0,fU,fs);
The numerator coefficients are: b = 0.9167 2.1554 3.1004 2.1554 0.9167
Now let’s quantize the numerator coefficients to 2^{13} = 8192 steps per unit of coefficient amplitude:
bq= round(b*8192)/8192;
The poles and zeros are then:
q= roots(bq); p= roots(a);
The amplitude response is computed using bq and a:
[h,f]= freqz(bq,a,Nfft,fs); H= 20*log10(abs(h));
The resulting zeros and amplitude response are shown in the top row of Figure 6. In the polezero plot, which shows just quadrant 1 of the zplane, you can just notice that the two zeros do not fall at exactly the same place. The depth of the null is about 45 dB. If we quantize to 4096 steps, we get the plots in the bottom row of Figure 6. Null depth is about 33 dB.
Quantizing to 2048 steps gives a null depth of only 25 dB, as shown in the top row of Figure 7. One way to get a deeper null without adding coefficient bits is to experiment with slight offsets of f_{0}, which may result in quantized zero locations that are closer to a point on the unit circle. For example, the bottom row of Figure 7 has 2048 steps and f_{0} = 15.03 Hz, resulting in a null depth of about 43 dB.
For lowpass prototype N >2 (bandreject order > 4), the effects of quantization of both numerator and denominator coefficients are more pronounced. Quantization of denominator coefficients can cause a pole to move outside the unit circle, making the filter unstable. In many cases, it may be advisable to break the filter up into a cascade of second order sections [5,6].
Figure 6. Quantization of numerator coefficients of BRF based on N= 2 lowpass.
Top row: quantization with 8192 steps per unit of coefficient amplitude
Bottom row: quantization with 4096 steps per unit of coefficient amplitude
Figure 7. Quantization of numerator coefficients of BRF based on N= 2 lowpass.
Top row: quantization with 2048 steps per unit of coefficient amplitude
Bottom row: quantization with 2048 steps with f_{0 }offset to 15.03 Hz
Appendix A. Matlab Functions br_synth1.m and br_synth2.m
These programs are provided asis 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 programs.
%function [b,a]= br_synth1(N,f0,fU,fs) 1/15/18 Neil Robertson % synthesize band reject IIR Butterworth filter with null frequency f0 and % upper 3 dB frequency fU. % % N= order of prototype LPF % f0= null frequency, Hz % fU= upper 3 dB frequency, Hz % fs= sample frequency, Hz function [b,a]= br_synth1(N,f0,fU,fs) if fU>=fs/2 error('fU must be less than fs/2') end % prewarp f0 and fu F0= fs/pi * tan(pi*f0/fs); FU= fs/pi * tan(pi*fU/fs); FL= F0^2/FU; % Hz analog lower 3 dB freq BW_hz= FUFL; % Hz analog 3 dB bandwidth %fL= fs/pi*atan(pi*FL/fs) % lower 3 dB freq % 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); % transform poles for brf centered at 2*pi*F0 % 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./(2*p_lp)).^2); pa= 2*pi*F0*[alpha+j*beta alphaj*beta]; % find poles and zeros of digital filter p= (1 + pa/(2*fs))./(1  pa/(2*fs)); % poles from bilinear transform w0= 2*pi*f0/fs; q= [exp(j*w0)*ones(1,N) exp(j*w0)*ones(1,N)]; % zeros at w0 on unit circle % convert poles and zeros to numerator and denominator polynomials a= poly(p); a= real(a); b= poly(q); b= real(b); % scale coeffs so amplitude is 1.0 at f= 0 K= sum(a)/sum(b); b= K*b;
%function [b,a]= br_synth2(N,fL,fU,fs) 1/15/18 Neil Robertson % synthesize band reject IIR Butterworth filter with 3 dB frequencies fL and %fU. Note null frequency f0 is offset from center frequency. % % N= order of prototype LPF % fL= lower 3 dB freq, Hz % fU= upper 3 dB freq, Hz % fs= sample frequency, Hz function [b,a]= br_synth2(N,fL,fU,fs) if fU>=fs/2 error('fU must be less than fs/2') end % prewarp fL and fU FL= fs/pi * tan(pi*fL/fs); FU= fs/pi * tan(pi*fU/fs); BW_hz= FUFL; % Hz analog 3 dB bandwidth F0= sqrt(FL*FU); % Hz analog null frequency f0= fs/pi*atan(pi*F0/fs); % Hz discrete null frequency % 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); % transform poles for brf centered at 2*pi*F0 % 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./(2*p_lp)).^2); pa= 2*pi*F0*[alpha+j*beta alphaj*beta]; % find poles and zeros of digital filter p= (1 + pa/(2*fs))./(1  pa/(2*fs)); % poles from bilinear transform w0= 2*pi*f0/fs; q= [exp(j*w0)*ones(1,N) exp(j*w0)*ones(1,N)]; % zeros at w0 on unit circle % convert poles and zeros to numerator and denominator polynomials a= poly(p); a= real(a); b= poly(q); b= real(b); % scale coeffs for gain of 1 at f= 0 K= sum(a)/sum(b); b= K*b;
Appendix B. Lowpass to Bandreject Frequency Transformation
The lowpass to bandreject transformation is the inverse of the lowpass to bandpass transformation we performed in the last post. In the sdomain, we want to transform a normalized lowpass filter with 3 dB frequency of 1 rad/s to a bandreject filter with a given bandwidth and null frequency [7,8]. Once we have the bandreject H(s), we can use the bilinear transform to obtain H(z). 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 two posts, we use Ω for continuous frequency and reserve ω for discrete frequency). To transform H(s’) into H(s), substitute the following for s’:
$$s'=\frac{BW}{\Omega_0} \frac{1}{\left(\frac{s}{\Omega_0}+\frac{\Omega_0}{s}\right)}\qquad B1$$
where BW = Ω_{U}Ω_{L} = 3 dB bandwidth and $\Omega_0=\sqrt{\Omega_L\Omega_U}$. As we will show, Ω_{0} is the null frequency of the filter. Rearranging B1, we have:
$$s'=BW\frac{s}{s_2+\Omega_0^2}\qquad B2$$
Then, solving for s:
$$s=\Omega_0\left[\frac{BW}{2\Omega_0s'} \pm j\sqrt{1\left(\frac{BW}{2\Omega_0s'}\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_0p'_a} \pm j\sqrt{1\left(\frac{BW}{2\Omega_0p'_a}\right)^2}\ \right]\qquad B3$$
For each lowpass pole p_{a}’ , we get two complexconjugate bandreject poles. Thus if the lowpass filter has order N, H(s) for the bandreject 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_0p'_a} \pm j\sqrt{1\left(\frac{BW_{Hz}}{2F_0p'_a}\right)^2}\ \right]\qquad B4$$
Oddly enough, the poles of Butterworth bandreject and bandpass filters with the same order, center frequency and bandwidth are identical. (This results from the lowpass prototype complexconjugate poles falling on a unit circle in the splane). Thus the feedback (denominator) coefficients of the bandpass and bandreject filters are the same.
Zeros of H(s) and H(z)
Consider a singlepole normalized lowpass transfer function in s’:
$$H(s')=\frac{1}{s'\alpha}$$
It has a zero at s’= infinity. From B2, we have for the corresponding
bandreject filter:
$$H(s)=\frac{1}{BW\frac{s}{s_2+\Omega_0^2}\alpha}$$
$$H(s)=\frac{s^2+\Omega_0^2}{BW*s\alpha(s^2+\Omega_0^2)}$$
$$H(s)=\frac{(s+j\Omega_0)(sj\Omega_0)}{BW*s\alpha(s^2+\Omega_0^2)}$$
H(s) has zeros at +jΩ_{0} and  jΩ_{0}. In general, if H(s’) has N zeros at infinity, then H(s) has N zeros at +jΩ_{0} and N zeros at  jΩ_{0}. Thus Ω_{0} is the null frequency of the analog bandreject filter.
In the zdomain, H(z) has zeros on the unit circle at exp(jω_{0}) and exp(jω_{0}), where ω_{0} is related to Ω_{0} by the frequency mapping of the bilinear transform [1]:
$$\Omega_0=2f_stan(\omega_0/2)$$
or
$$\omega_0=2tan^{1}\left(\frac{\Omega_0}{2f_s}\right)$$
We can now write H(z) in polezero format
as:
$$H(z)=K\frac{(z e^{j\omega_0})^N(ze^{j\omega_0})^N}{(zp_1)(zp_2)...(zp_{2N})}\qquad B5$$
A zplane polezero plot illustrating the zeros on the unit circle is shown in Example 1 in the main text (Figure 2).
In the Matlab functions br_synth1.m and br_synth2.m, we expand the numerator and denominator of H(z) to get the filter coefficients. Note you can begin to expand the numerator by hand as follows:
$$numerator=[z^2z(e^{j\omega_0} +e^{j\omega_0})+1]^N$$
$$=[z^22cos(\omega_0)z+1]^N$$
It then turns out that the numerator
coefficients are symmetric with respect to the center coefficient, and the
first and last coefficients are unity (before scaling by K).
Appendix C Notation
Notation for normalized lowpass continuous frequencies:
normalized radian frequency Ω’ rad/s
normalized complex frequency s’= σ + jΩ’
Notation for bandreject filter continuous frequencies:
3 dB frequencies Ω_{L }, Ω_{U} rad/s
3 dB bandwidth BW rad/s
null frequency $\Omega_0=\sqrt{\Omega_L\Omega_U}$ rad/s
 3 dB bandwidth, Hz BW_{Hz} Hz
null frequency, Hz F_{0} Hz
Notation for bandreject filter discrete frequencies:
3 dB frequencies, Hz f_{L }, f_{U} Hz
null frequency, Hz f_{0} Hz
General notation:
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 functions represent 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. Robertson, Neil , “Design IIR Bandpass Filters”, Jan 2017 https://www.dsprelated.com/showarticle/1128.php
3. Williams, Arthur B. and Taylor, Fred J., Electronic Filter Design Handbook, 3^{rd} Ed., McGrawHill, 1995, section 2.3
4. Analog Devices Mini Tutorial MT224, 2012 http://www.analog.com/media/en/trainingseminars/tutorials/MT224.pdf
5. Oppenheim, Alan V. and Shafer, Ronald W., DiscreteTime Signal Processing, Prentice Hall, 1989, sections 6.8 and 6.9.
6. Lyons, Richard G., Understanding Digital Signal Processing, 2^{nd} Ed., Pearson, 2004, sections 6.7 and 6.8
7. Blinchikoff, Herman J., and Zverev,Anatol I., Filtering in the Time and Frequency Domains, Wiley, 1976, section 4.7.
8. 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
Neil Robertson January, 2018
Previous post by Neil Robertson:
Design IIR Bandpass Filters
Next post by Neil Robertson:
Design IIR Highpass Filters
Comments:
 Comments
 Write a Comment Select to add a comment
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.