Filtering and Scipy.Signal

Demonstration of FIR filtering using scipy.signal. Part of a lecture given in ELEG 306 Digital Signal Processing, Oct 6, 2014.

Charles Boncelet, boncelet@udel.edu

In [1]: 
%pylab inline
import scipy.signal as sig

Populating the interactive namespace from numpy and matplotlib
In [2]: 
x = randn(1200)
In [56]: 
plot(x)
Out[56]:
[<matplotlib.lines.Line2D at 0x1196a6390>]

Low Pass filter

In [55]: 
#sig.firwin?
In [6]: 
h = sig.firwin(31,0.1)
plot(h)
Out[6]:
[<matplotlib.lines.Line2D at 0x10ce9ec10>]

The connected dots plot above isn't as informative as the stem plot below:

In [8]: 
stem(range(len(h)),h)
Out[8]:
<Container object of 3 artists>

Below we plot the magnitude response of the filter versus frequency (only versus positive frequencies--the negative half is the mirror image).

Notice, this is a linear phase filter whose sidelobes are about 120 dB below the passband. Try building an analog filter like this!

In [11]: 
w, v = sig.freqz(h)
plot(w,20*log(abs(v)))
Out[11]:
[<matplotlib.lines.Line2D at 0x10dfc8450>]
In [12]: 
y=convolve(x,h)
plot(y)
Out[12]:
[<matplotlib.lines.Line2D at 0x10e415490>]

Below, we tighten the filter to have a passband of 0.01\(\pi\) wide (ok, as well as we can do with 51 taps--more taps results in a response closer to ideal at the expense of more delay and more computation).

In [13]: 
h2 = sig.firwin(51,0.01)
stem(range(len(h2)),h2)
Out[13]:
<Container object of 3 artists>

Notice how the signal is much smoother now. This is clearer in the scatter plot below.

One easily missed point: The y-axis autoscales, resulting in plots of about the same "visual magnitude". But look carefully: the original plot of x has a scale from -4 to 4 while this one goes from about -0.4 to 0.4. The signal has gotten much smaller. Why? Because the filter has removed about 99% of the energy!

In [14]: 
y2 = convolve(x,h2)
plot(y2)
Out[14]:
[<matplotlib.lines.Line2D at 0x10e73f210>]

Show the dots. Notice how close consecutive dots are to each other. Also, notice how the dots have "momentum". This is a characteristic of a low pass signal like this.

In [17]: 
scatter(range(len(y2)),y2,s=2)
Out[17]:
<matplotlib.collections.PathCollection at 0x10f1e0450>

Bandpass Filter

Notice how the bandpass filter's impulse response matches the intended signal. It oscillates slowly between positive and negative values, much like a sinusoid of frequency 0.15\(\pi\) does.

In [20]: 
hbp = sig.firwin(51, [0.1,0.2], pass_zero=False)
stem(range(len(hbp)),hbp)
Out[20]:
<Container object of 3 artists>
In [21]: 
w, v = sig.freqz(hbp)
plot(w, 20*log(abs(v)))
Out[21]:
[<matplotlib.lines.Line2D at 0x10f7aa3d0>]
In [22]: 
ybp = convolve(x,hbp)
plot(ybp)
Out[22]:
[<matplotlib.lines.Line2D at 0x10f7dd2d0>]
In [25]: 
stem(range(len(ybp)),ybp,s=1)
Out[25]:
<Container object of 3 artists>

This plot is a total mess. It has too many points and the dots are too large. Let's fix both:

In [51]: 
segment = ybp[100:200]
stem(range(len(segment)),segment,s=2)
Out[51]:
<Container object of 3 artists>

Now we can see what happened. The signal looks like a sinusoid, oscillating slowly between positive and negative values. It isn't a perfect sinusoid for two reasons: 1) the signal is random (we starting with white Gaussian noise) and 2) the output has components at all frequencies within (and near) the passband.

High Pass Filter

In [31]: 
hhp = sig.firwin(51, 0.9, pass_zero=False)
w, v = sig.freqz(hhp)
plot(w, 20*log(abs(v)))
Out[31]:
[<matplotlib.lines.Line2D at 0x1184e2110>]
In [32]: 
yhp = convolve(x,hhp)
plot(yhp)
Out[32]:
[<matplotlib.lines.Line2D at 0x1189d1850>]
In [33]: 
plot(yhp[500:600])
Out[33]:
[<matplotlib.lines.Line2D at 0x1189fddd0>]
In [34]: 
stem(yhp[500:600])
Out[34]:
<Container object of 3 artists>

Pole-Zero Plots

Let's do some pole-zero plots. Since the filters have no poles, we only have to worry about the zeros.

In [61]: 
z, p, k = sig.tf2zpk(h,1)
z
Out[61]:
array([ 1.53330017+0.j        , -0.99486419+0.10121875j,
       -0.99486419-0.10121875j, -0.95409239+0.29951247j,
       -0.95409239-0.29951247j, -0.87421228+0.48554391j,
       -0.87421228-0.48554391j, -0.75848152+0.65169455j,
       -0.75848152-0.65169455j, -0.61161532+0.7911553j ,
       -0.61161532-0.7911553j , -0.43958666+0.89820018j,
       -0.43958666-0.89820018j, -0.24936678+0.96840911j,
       -0.24936678-0.96840911j, -0.04860325+0.99881816j,
       -0.04860325-0.99881816j,  0.15479768+0.98794619j,
        0.15479768-0.98794619j,  0.35338136+0.93547935j,
        0.35338136-0.93547935j,  0.65218802+0.j        ,
        0.54275400+0.83989172j,  0.54275400-0.83989172j,
        0.81062009+0.77712424j,  0.81062009-0.77712424j,
        0.76018662+0.64970478j,  0.76018662-0.64970478j,
        0.64282531+0.61626296j,  0.64282531-0.61626296j])
In [49]: 
theta = linspace(-pi,pi,201)
plot(cos(theta),sin(theta))
fig = scatter(real(z),imag(z))
axes().set_aspect('equal')
In [50]: 
z, p, k = sig.tf2zpk(hbp,[1])
theta = linspace(-pi,pi,201)
plot(cos(theta),sin(theta))
fig = scatter(real(z),imag(z))
axes().set_aspect('equal')

Remember, zeros can be outside the unit circle without affecting stability.

In []: