Transfer functions for Scipy and Matlab's Butterworth filter don't seem to match theory

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')
ax.legend() 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)')
ylabel('|H(f)|^2')
• 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)$ Jun 18 '21 at 17:31
• 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. Jun 18 '21 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:

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
figure;
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)');
ylabel('|H(f)|^2');
axis tight;

% Theory vs analog Butterworth
figure;
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)');
ylabel('|H(f)|^2');
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: Here are some methods I've bumped to:

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