Instantly share code, notes, and snippets.

# endolith/Butterworth_Lowpass_Filter_Example.png Last active Dec 11, 2015

Second-order sections for SciPy Python
 """ 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]))