# How do I manually plot the frequency response of a bandpass Butterworth filter in MATLAB without freqz function?

I have code like below that applies a bandpass filter onto a signal. I am quite a noob at DSP and I want to understand what is going on behind the scenes before I proceed.

To do this, I want to know how to plot the frequency response of the filter without using freqz.

[b, a] = butter(order, [flo fhi]);
filtered_signal = filter(b, a, unfiltered_signal)


Given the outputs [b, a] how would I do this? This seems like it would be a simple task, but I'm having a hard time finding what I need in the documentation or online.

I'd also like to understand how to do this as quickly as possible, e.g. using an fft or other fast algorithm.

We know that in general transfer function of a filter is given by:

$$H(z)=\dfrac{\sum_{k=0}^{M}b_kz^{-k}}{\sum_{k=0}^{N}a_kz^{-k}}$$

Now substitute $$z=e^{j\omega}$$ to evaluate the transfer function on the unit circle:

$$H(e^{j\omega})=\dfrac{\sum_{k=0}^{M}b_ke^{-j\omega k}}{\sum_{k=0}^{N}a_ke^{-j\omega k}}$$

Thus this becomes only a problem of polynomial evaluation at a given $$\omega$$. Here are the steps:

• Create a vector of angular frequencies $$\omega = [0, \ldots,\pi]$$ for the first half of spectrum (no need to go up to $$2\pi$$) and save it in w.
• Pre-compute exponent $$e^{-j\omega}$$ at all of them and store it in variable ze.
• Use the polyval function to calculate the values of numerator and denominator by calling: polyval(b, ze), divide them and store in H. Because we are interested in amplitude, then take the absolute value of the result.
• Convert to dB scale by using: $$H_{dB}=20\log_{10} H$$ - in this case $$1$$ is the reference value.

Putting all of that in code:

%% Filter definition
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
b = [1 3 2];

%% My freqz calculation
N = 1024; % Number of points to evaluate at
upp = pi; % Evaluate only up to fs/2
% Create the vector of angular frequencies at one more point.
% After that remove the last element (Nyquist frequency)
w = linspace(0, pi, N+1);
w(end) = [];
ze = exp(-1j*w); % Pre-compute exponent
H = polyval(b, ze)./polyval(a, ze); % Evaluate transfer function and take the amplitude
Ha = abs(H);
Hdb  = 20*log10(Ha); % Convert to dB scale
wn   = w/pi;
% Plot and set axis limits
xlim = ([0 1]);
plot(wn, Hdb)
grid on

%% MATLAB freqz
figure
freqz(b,a)


Original output of freqz: And the output of my script: And quick comparison in linear scale - looks great!

[h_f, w_f] = freqz(b,a);
figure
xlim = ([0 1]);
plot(w, Ha) % mine
grid on
hold on
plot(w_f, abs(h_f), '--r') % MATLAB
legend({'my freqz','MATLAB freqz'}) Now you can rewrite it into some function and add few conditions to make it more useful.

Another way (previously proposed is more reliable) would be to use the fundamental property, that frequency response of a filter is a Fourier Transform of its impulse response:

$$H(\omega)=\mathcal{F}\{h(t)\}$$

Therefore you must feed into your system $$\delta(t)$$ signal, calculate the response of your filter and take the FFT of it:

d = [zeros(1,length(w_f)) 1 zeros(1,length(w_f)-1)];
h = filter(b, a, d);
HH = abs(fft(h));
HH = HH(1:length(w_f));


In comparison this will produce the following: • Detailed explanation – partida Dec 13 '14 at 14:49
• I am thinking this line a = [1 -0.5 -0.25]; % Some filter with lot's of static gain. Can you explain me the selection of these parameters here, please. I am reading my Matlab's manual and it says [h,w] = freqz(hfilt,n) in one part of synapse. You are giving two filters (a,b) into freqz. Are both filters in hfilt? Or one in n? – Léo Léopold Hertz 준영 Apr 21 '15 at 21:57
• Just to clarify for others: "No need to go up to 2 pi" is when the coefficients are real. There are applications for filters with complex coefficients and in that case the spectrum will no longer be symmetric and would in that case want to extend the frequency to 2 pi. – Dan Boschen May 14 '17 at 14:10
• shouldn't be that ze=exp(+jw) since the the formula of H(exp(+jw))=( b0*exp(2jw) + b1*exp(1jw) + b2 )/( exp(2*jw) +a1*exp(jw) + a2 ). then if we multiply it by exp(-2*jw), we get back to the original formula? also this fits exactly the polyval that has the highest order from the left as if we use exp(+j*w) – Yaakov May 8 '20 at 12:53

this is just an addendum to jojek's answer which is more general and perfectly good when double-precision math is used. when there is less precision, there is a "cosine problem" that crops up when either the frequency in the frequency response is very low (much lower than Nyquist) and also when the resonant frequencies of the filter are very low.

when you compute the magnitude (squared) frequency response $|H(e^{j\omega})|^2$ these complex exponentials will be converted to sines and cosines, but when the math is cranked out, only the cosines will survive. this should be clear because the magnitude is an even function $|H(e^{-j\omega})| = |H(e^{j\omega})|$ w.r.t. frequency and only the cosines are even functions.

consider this trig identity:

$$\cos(\omega) \ = \ 1 - 2 \sin^2 \left( \frac{\omega}{2} \right)$$

now, looking at the right-hand side of the equation, all of the information regarding frequency is in the $\sin^2 \left( \frac{\omega}{2} \right)$ term which is being subtracted from 1. and this term gets exceedingly small as $\omega \to 0$. so small that the term and the frequency information in that term are getting lost when the term (or it's negative) are added to 1. this is the case even with floating point, but is less a problem with double-precision floats. but some of us putting this frequency response function into gear might not have the resource of double-precision floating point, or any floating point.

so, what i have done is use the trig identity above and eliminate all of the cosine terms, essentially replacing them with terms looking like $\sin^2 \left( \frac{\omega}{2} \right)$ and some constants that get combined with other constants. i'll show you the answer for the case of a 2nd-order IIR filter (a.k.a. a "biquad"):

$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{a_0 + a_1 z^{-1} + a_2 z^{-2}}$$

which has complex frequency response

$$H(e^{j\omega}) = \frac{b_0 + b_1 e^{-j\omega} + b_2 e^{-j2\omega}}{a_0 + a_1 e^{-j\omega} + a_2 e^{-j2\omega}}$$

which has magnitude squared:

\begin{align} |H(e^{j\omega})|^2 &= \frac{|b_0 + b_1 e^{-j\omega} + b_2 e^{-j2\omega}|^2}{|a_0 + a_1 e^{-j\omega} + a_2 e^{-j2\omega}|^2} \\ \\ &= \frac{\big(b_0 + b_1\cos(\omega) + b_2\cos(2\omega)\big)^2 + \big(b_1\sin(\omega) + b_2\sin(2\omega)\big)^2}{\big(a_0 + a_1\cos(\omega) + a_2\cos(2\omega)\big)^2 + \big(a_1\sin(\omega) + a_2\sin(2\omega)\big)^2} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\cos(\omega) + 2b_0b_2\cos(2\omega)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\cos(\omega) + 2a_0a_2\cos(2\omega)} \\ \end{align}

so, one can see that the magnitude frequency response $|H(e^{j\omega})|$ is an even symmetry function and depends only on the cosines $\cos(\omega)$ and $\cos(2\omega)$. for very low $\omega$, the values of those cosines are so close to $1$ that, with single-precision fixed or floating point, there are few bits remaining that differentiate those values from $1$. that is the "cosine problem".

using the trig identity above, you get for magnitude squared:

\begin{align} |H(e^{j\omega})|^2 &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\cos(\omega) + 2b_0b_2\cos(2\omega)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\cos(\omega) + 2a_0a_2\cos(2\omega)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2b_0b_2\left(1 - 2 \sin^2(\omega)\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2a_0a_2\left(1 - 2 \sin^2(\omega)\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2b_0b_2\left(2\cos^2(\omega) - 1\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2a_0a_2\left(2\cos^2(\omega) - 1\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2b_0b_2\left(2\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right)^2 - 1\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2a_0a_2\left(2\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right)^2 - 1\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)(1 - 2\phi) + 2b_0b_2\left(2(1 - 2\phi )^2 - 1\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)(1 - 2\phi) + 2a_0a_2\left(2(1 - 2\phi)^2 - 1\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)(1 - 2\phi) + 2b_0b_2(1 - 8\phi + 8\phi^2)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)(1 - 2\phi) + 2a_0a_2(1 - 8\phi + 8\phi^2)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1b_0+2b_1b_2 - 4(b_1b_0+b_1b_2)\phi + 2b_0b_2 - 16b_0b_2\phi + 16b_0b_2\phi^2}{a_0^2+a_1^2+a_2^2 + 2a_1a_0+2a_1a_2 - 4(a_1a_0+a_1a_2)\phi + 2a_0a_2 - 16a_0a_2\phi + 16a_0a_2\phi^2} \\ \\ &= \frac{\big(b_0^2+b_1^2+b_2^2 + 2b_1b_0+2b_1b_2+2b_0b_2\big) - 4(b_1b_0+b_1b_2-4b_0b_2)\phi + 16b_0b_2\phi^2}{\big(a_0^2+a_1^2+a_2^2 + 2a_1a_0+2a_1a_2+2a_0a_2\big) - 4(a_1a_0+a_1a_2-4a_0a_2)\phi + 16a_0a_2\phi^2} \\ \\ &= \frac{\tfrac14\big(b_0^2+b_1^2+b_2^2 + 2b_1b_0+2b_1b_2+2b_0b_2\big) - (b_1b_0+b_1b_2-4b_0b_2)\phi + 4b_0b_2\phi^2}{\tfrac14\big(a_0^2+a_1^2+a_2^2 + 2a_1a_0+2a_1a_2+2a_0a_2\big) - (a_1a_0+a_1a_2-4a_0a_2)\phi + 4a_0a_2\phi^2} \\ \\ &= \frac{\left(\frac{b_0+b_1+b_2}{2}\right)^2 - \phi \big(4b_0b_2(1-\phi) + b_1(b_0+b_2)\big)}{\left(\frac{a_0+a_1+a_2}{2}\right)^2 - \phi \big(4a_0a_2(1-\phi) + a_1(a_0+a_2)\big)} \\ \end{align}

where $\phi \triangleq \sin^2 \left( \frac{\omega}{2} \right)$

if your gear is intending to plot this as dB, it comes out as

$$20 \log_{10}|H(e^{j\omega})| \ = \ 10 \log_{10}\left( \left(\tfrac{b_0+b_1+b_2}{2}\right)^2 - \phi \big(4b_0b_2(1-\phi) + b_1(b_0+b_2)\big) \right) \\ - 10 \log_{10}\left(\left(\tfrac{a_0+a_1+a_2}{2}\right)^2 - \phi \big(4a_0a_2(1-\phi) + a_1(a_0+a_2)\big) \right)$$

so your division turns into subtraction, but you have to be able to compute logarithms to some base or another. numerically, you will have much less trouble with this for low frequencies than doing it the apparent way.

• That's really cool, thank you Robert! +1 – jojek Jul 11 '14 at 11:29
• @Robert I "believe" similar to my comment for Jojek above that this only applies as well when the coefficients are real (and therefore the spectrum is symmetric and thus the magnitude converts to cosines as you show)... Am I correct? – Dan Boschen May 14 '17 at 14:12
• yes. that commitment is made when you go from the first line of $|H(e^{j\omega})|^2 = ...$ to the second line. no going back after that. – robert bristow-johnson May 14 '17 at 18:38