The squared transfer function for a Butterworth filter of order $n$ should be

$$ |H(f)|^2 = \frac{1}{1+\left(\frac{f}{f_c}\right)^{2n}} $$

where $f_c$ is the cut-off frequency. (Here's one of many possible references.)

I'm unable to correctly reproduce this with either Python or Matlab. For example, let's say I want a third-order filter ($n=3$) with a $f_\mathrm{c}$ = 20 Hz cut-off frequency for a signal sampled at $f_\mathrm{s}$ = 100 Hz.

In Python, I would implement this as

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfreqz

fc = 20  # Hz, Cut-off frequency
fs = 100  # Hz, sampling frequency
order = 3  # Filter order

nyq = 0.5*fs  # Nyquist frequency
sos = butter(order, fc/nyq, analog=False, btype='low', output='sos')

w, h = sosfreqz(sos, worN=2000)
f = nyq*w/np.pi  # Convert to Hz

Hf2_scipy = np.abs(h)**2
Hf2_theory = 1/(1+(f/fc)**(2*order))

Plotting the result gives curves that don't match

fig, ax = plt.subplots()
ax.plot(f, Hf2_scipy, label='Scipy')
ax.plot(f, Hf2_theory, label='Theory')
ax.set(xlabel='f (Hz)', ylabel='|H(f)|^2', xlim=[0, nyq], ylim=[0, 1])
ax.axhline(0.5, color='0.5')
ax.axvline(fc, color='0.5')

Butterworth transfer function plot

What am I missing?

Here's the equivalent attempt in Matlab, which produces the same result

fc = 20;  % Cut-off frequency
fs = 100;  % Sample frequency
order = 3;  % Filter order

[z, p, k] = butter(order, fc/(fs/2), 'low');
sos = zp2sos(z, p, k);
[h, w] = freqz(sos, 2000);
f = w*fs/(2*pi); % Convert to Hz

Hf2_matlab = abs(h).^2;
Hf2_theory = 1./(1+(f/fc).^(2*order));
plot(f, Hf2_matlab, 'k')
plot(f, Hf2_theory)

plot([fc, fc], [0, 1], 'color', 0.5*[1, 1, 1])
plot([0, fs/2], [0.5, 0.5], 'color', 0.5*[1, 1, 1])
legend('Matlab', 'Theory')
xlabel('f (Hz)')
  • 2
    $\begingroup$ Bilinear warping of the frequency axis. The bilinear transform maps the imaginary axis of the s-plane onto the unit circle of z-plane. That means the digital frequency $\omega = \pi$ maps to the analog frequency of infinity, not $fs/2$. That's why the digital lowpass filter has $H_z(\pi) = 0$ and not $H_z(\pi) = H_s(fs/2)$ $\endgroup$ – Hilmar Jun 18 at 17:31
  • 1
    $\begingroup$ A simpler corollary of @Hilmar's comment is that my discrepancy arose because the cut-off frequency was comparable to the Nyquist frequency (20 Hz vs 50 Hz). If I use a sampling frequency, say, 10 times larger, then the theory matches Matlab/Python result. $\endgroup$ – hugke729 Jun 18 at 18:45

It has to do with the Bilinear Transform, as Hilmar already stated.

The theoretical function is defined for the analog domain, so inherently there will be differences when you convert the response to the digital domain. However, you can still generate the "analog" frequency response and yield the expected result.

I took your code and added/modified it a bit to produce two plots: each comparing theory against the digital or analog version of the filter. This is done in MATLAB:

enter image description here enter image description here

And here's the code:

%% Butterworth filter definitions

fc = 20;  % Cut-off frequency
fs = 100;  % Sample frequency
order = 3;  % Filter order

% Zero-pole gain definition for digital Butterworth filter
[z, p, k] = butter(order, fc/(fs/2), 'low');
sos = zp2sos(z, p, k);
[h, f] = freqz(sos, 2000, fs);

% Polynomial coefficients (alternate definition)
% [b, a] = butter(order, fc/(fs/2), 'low');
% [h, f] = freqz(b, a, 2000, fs);

% Analog prototype
[za, pa, ka] = buttap(order);       % Butterworth prototype
[b, a] = zp2tf(za, pa, ka);
[bt, at] = lp2lp(b, a, 2*pi*fc);    % Change the cut-off frequency
ha = freqs(bt, at, 2*pi.*f);        % Analog frequency response

% Define the squared-magnitude transfer functions
Hf2_matlabDigital = abs(h).^2;
Hf2_matlabAnalog = abs(ha).^2;
Hf2_theory = 1./(1 + (f./fc).^(2*order));

%% Plot comparisons

% Theory vs digital Butterworth
hold on;
plot(f, Hf2_theory);
plot(f, Hf2_matlabDigital, '--k');
plot([fc, fc], [0, 1], 'color', 0.5*[1, 1, 1]);
plot([0, fs/2], [0.5, 0.5], 'color', 0.5*[1, 1, 1]);
legend('Theory', 'Matlab (Digital)');
xlabel('f (Hz)');
axis tight;

% Theory vs analog Butterworth
hold on;
plot(f, Hf2_theory);
plot(f, Hf2_matlabAnalog, '--k');
plot([fc, fc], [0, 1], 'color', 0.5*[1, 1, 1]);
plot([0, fs/2], [0.5, 0.5], 'color', 0.5*[1, 1, 1]);
legend('Theory', 'Matlab (Analog)');
xlabel('f (Hz)');
axis tight;

FYI, if needed, it is possible to match magnitude or phase of digital filter to follow quite well the theory/analog model by using certain methods.

I did this example plot using MIM (Magnitude Invariance Method) method:

enter image description here

Here are some methods I've bumped to:

MIM/PIM https://ieeexplore.ieee.org/document/4629708 , https://ieeexplore.ieee.org/document/7357808 , https://soar.wichita.edu/handle/10057/1564 (implementation)

Massberg https://books.google.fi/books?id=QddcxHLavrMC&pg=PA201&lpg=PA204#v=onepage&q&f=false

Orfandis https://www.semanticscholar.org/paper/Digital-Parametric-Equalizer-Design-with-Prescribed-Orfanidis/b98f2ac93a6ef9b60e109960c6f2f00edde9bf64 , (check särkkä .../pub/ folder (https://users.aalto.fi/~ssarkka/pub/eq-design-demo.zip)

Vicanek https://vicanek.de/articles/BiquadFits.pdf

Särkkä https://users.aalto.fi/~ssarkka/pub/eq-article.pdf

MZTi http://www.khabdha.org/wp-content/uploads/2008/03/optimizing-the-magnitude-response-of-mzt-filters-mzti-2007.pdf , https://www.kvraudio.com/forum/viewtopic.php?t=360690

Jackson https://ieeexplore.ieee.org/document/870677?tp=&arnumber=870677

Mecklenbräuker https://www.sciencedirect.com/science/article/abs/pii/S0165168400001134

Nelatury https://www.sciencedirect.com/science/article/abs/pii/S1051200406001424

Flynn/Reiss http://www.eecs.qmul.ac.uk/~josh/documents/2018/19412.pdf etc.


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.