Second-order sections for SciPy Python

readme.md

sos functions have been rewritten and moved into SciPy: scipy/scipy#3717

butter_sos_example.py

"""
Translation of example from:
https://ccrma.stanford.edu/~jos/fp/Butterworth_Lowpass_Filter_Example.html
to demonstrate that sosfilt and tf2sos function the same as the Octave version

"""

from __future__ import division
from scipy.signal import butter
from sos import tf2sos, sosfilt
from numpy import shape, zeros, log10
from numpy.fft import fft, fftfreq
from matplotlib.pyplot import plot, figure, axis, grid

fc = 1000 # Cut-off frequency (Hz)
fs = 8192 # Sampling rate (Hz)
order = 5 # Filter order
B, A = butter(order, fc / (fs/2)) # [0:pi] maps to [0:1] here
sos, k = tf2sos(B, A)
print sos
print k

"""
Actual values differ slightly from Octave's, which are:

sos = array([[1.00000,  2.00080,   1.00080,  1.00000,  -0.92223,  0.28087],
             [1.00000,  1.99791,   0.99791,  1.00000,  -1.18573,  0.64684],
             [1.00000,  1.00129,  -0.00000,  1.00000,  -0.42504,  0.00000]])

g = 0.0029714
"""

## Compute and display the amplitude response
#Bs = sos[:, :3] # Section numerator polynomials
#As = sos[:, 3:] # Section denominator polynomials
#nsec = shape(sos)[0]

nsamps = 256 # Number of impulse-response samples

# Note use of input scale-factor k here:
x = zeros(nsamps)
x[0] = k # SCALED impulse signal

x = sosfilt(sos, x)

# Plot impulse response to make sure it has decayed to zero (numerically)
plot(x)

# Plot amplitude of frequency response 
figure(2)
X = fft(x) # sampled frequency response
f = fftfreq(nsamps, 1/fs)
grid(True)

axis([0, fs / 2, -100, 5])

plot(f[:nsamps / 2], 20 * log10(X[:nsamps / 2]))

Butterworth_Lowpass_Filter_Example.png