MATLAB and SciPy give different results for 'buttord' function

Question

I'm trying to design an analog Butterworth filter using the buttord function (actually, I'm porting a program where this function is called from MATLAB to Python).

My parameters are:

Passband frequency (Fp) = 10 Hz, giving Wp = 2*pi*10 Hz

Stopband frequency (Fs) = 100 Hz, giving Ws = 2*pi*100 Hz

The passband and stopband losses/attenuations (Rp, Rs) are 3 and 80 dB respectively.

In MATLAB I use this code:

Wp = 2 * pi * 10
Ws = 2 * pi * 100
Rp = 3
Rs = 80
[N, Wn] = buttord(Wp, Ws, Rp, Rs, 's')

that gives me N = 5, Wn = 99.581776302.

In SciPy I tried to do the same:

from numpy import pi
from scipy import signal
Wp = 2 * pi * 10
Ws = 2 * pi * 100
Rp = 3
Rs = 80
(N, Wn) = signal.buttord(Wp, Ws, Rp, Rs, analog=True)

and I get N = 5 and Wn = 62.861698649592753. Wn is different than the value that MATLAB gives, and is strangely close to Wp. What is wrong here?

Digging into SciPy's sources and issues, I found this pull request which might explain things: turns out MATLAB and SciPy have different design goals (MATLAB tries to optimize for matching the stopband frequency and SciPy tries to optimize for matching the passband frequency).

I'm using MATLAB R2013a, Python 3.4.2 and SciPy 0.15.0 if that matters.

Matlab signal processing tools typically require frequency related parameters be expressed as normalized frequency with respect to the Nyquist. That is, Wp = Fp / Fn and Ws = Fs / Fn, where Fn is your Nyquist frequency. — Juderb
@Juderb I'm working with analog filters where the Nyquist frequency doesn't matter. From SciPy's doc, For analog filters, wp and ws are angular frequencies (e.g. rad/s). — Renan

Warren Weckesser Warren Weckesser · Accepted Answer · 2015-01-21T17:15:13

(I also posted the following on the scipy mailing list.)

When you design a Butterworth filter with buttord, there aren't enough degrees of freedom to meet all the design constraints exactly. So there is a choice of which end of the transition region hits the constraints and which end is "over-designed". A change made in scipy 0.14.0 switched that choice from the stop-band edge to the pass-band edge.

A picture will make it clear. The script below generates the following plot. (I changed Rp from 3 to 1.5. -3 dB coincides with the gain at Wn, that's why your Wn was the same as Wp.) The filters generated using either the old or new convention both satisfy the design constraints. With the new convention, the response just bumps against the constraint at the end of the pass-band.

frequency response

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


# Design results for:
Wp = 2*np.pi*10
Ws = 2*np.pi*100
Rp = 1.5      # instead of 3
Rs = 80

n_old = 5
wn_old = 99.581776302787929

n_new, wn_new = buttord(Wp, Ws, Rp, Rs, analog=True)

b_old, a_old = butter(n_old, wn_old, analog=True)
w_old, h_old = freqs(b_old, a_old)

b_new, a_new = butter(n_new, wn_new, analog=True)
w_new, h_new = freqs(b_new, a_new)


db_old = 20*np.log10(np.abs(h_old))
db_new = 20*np.log10(np.abs(h_new))

plt.semilogx(w_old, db_old, 'b--', label='old')
plt.axvline(wn_old, color='b', alpha=0.25)
plt.semilogx(w_new, db_new, 'g', label='new')
plt.axvline(wn_new, color='g', alpha=0.25)

plt.axhline(-3, color='k', ls=':', alpha=0.5, label='-3 dB')

plt.xlim(40, 1000)
plt.ylim(-100, 5)

xbounds = plt.xlim()
ybounds = plt.ylim()
rect = plt.Rectangle((Wp, ybounds[0]), Ws - Wp, ybounds[1] - ybounds[0],
                     facecolor="#000000", edgecolor='none', alpha=0.1, hatch='//')
plt.gca().add_patch(rect)
rect = plt.Rectangle((xbounds[0], -Rp), Wp - xbounds[0], 2*Rp,
                     facecolor="#FF0000", edgecolor='none', alpha=0.25)
plt.gca().add_patch(rect)
rect = plt.Rectangle((Ws, ybounds[0]), xbounds[1] - Ws, -Rs - ybounds[0],
                     facecolor="#FF0000", edgecolor='none', alpha=0.25)
plt.gca().add_patch(rect)

plt.annotate("Pass", (0.5*(xbounds[0] + Wp), Rp+0.5), ha='center')
plt.annotate("Stop", (0.5*(Ws + xbounds[1]), -Rs+0.5), ha='center')
plt.annotate("Don't Care", (0.1*(8*Wp + 2*Ws), -Rs+10), ha='center')

plt.legend(loc='best')
plt.xlabel('Frequency [rad/s]')
plt.ylabel('Gain [dB]')
plt.show()

MATLAB and SciPy give different results for 'buttord' function

1 Answers