Functions for frequency, phase, impulse and step response

Let's first define functions to plot filter properties.

from pylab import *
import scipy.signal as signal

#Plot frequency and phase response
def mfreqz(b,a=1):
    w,h = signal.freqz(b,a)
    h_dB = 20 * log10 (abs(h))
    subplot(211)
    plot(w/max(w),h_dB)
    ylim(-150, 5)
    ylabel('Magnitude (db)')
    xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
    title(r'Frequency response')
    subplot(212)
    h_Phase = unwrap(arctan2(imag(h),real(h)))
    plot(w/max(w),h_Phase)
    ylabel('Phase (radians)')
    xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
    title(r'Phase response')
    subplots_adjust(hspace=0.5)

#Plot step and impulse response
def impz(b,a=1):
    l = len(b)
    impulse = repeat(0.,l); impulse[0] =1.
    x = arange(0,l)
    response = signal.lfilter(b,a,impulse)
    subplot(211)
    stem(x, response)
    ylabel('Amplitude')
    xlabel(r'n (samples)')
    title(r'Impulse response')
    subplot(212)
    step = cumsum(response)
    stem(x, step)
    ylabel('Amplitude')
    xlabel(r'n (samples)')
    title(r'Step response')
    subplots_adjust(hspace=0.5)

Lowpass FIR filter

Designing a lowpass FIR filter is very simple to do with SciPy, all you need to do is to define the window length, cut off frequency and the window.

The Hamming window is defined as: \(w(n) = \alpha - \beta\cos\frac{2\pi n}{N-1}\), where \(\alpha=0.54\) and \(\beta=0.46\)

The next code chunk is executed in term mode, see the Python script for syntax. Notice also that Pweave can now catch multiple figures/code chunk.

print("I'm publishing a term chunk")

I'm publishing a term chunk

Let's define a highpass FIR filter, if you compare to original blog post you'll notice that it has become easier since 2009. You don't need to do ' spectral inversion "manually" anymore!

a = 12
b = 10
print(a+b)

\[ \sum_{i=1}^n x_i \]