%pylab inline
rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division
Populating the interactive namespace from numpy and matplotlib
sig = random.random(2048)*2 -1
Pxx, freqs, times, im = specgram(sig);
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x057A9AA8>
Pxx.shape
(129, 15)
2048/128
16.0
plot(Pxx[:,10])
[<matplotlib.lines.Line2D at 0x599d950>]
filtered = (sig + r_[0,sig[:-1]])
Pxx2, freqs, times, im = specgram(filtered);
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x05BAC1C0>
plot(Pxx2[:,10])
[<matplotlib.lines.Line2D at 0x71738b0>]
plot(Pxx[:,10])
plot(Pxx2[:,10])
[<matplotlib.lines.Line2D at 0x7154f30>]
plot(Pxx[:,10])
plot(Pxx2[:,10]/4.0)
[<matplotlib.lines.Line2D at 0x7f72708bad90>]
Difference equation:
$$y(n) = b_0x(n) + b_1x(n-1) + ... + b_Mx(n-M)$$for this simple filter:
$$y(n) = x(n) + x(n-1)$$You can estimate the effect of frequency from the difference equation!
from scipy.signal import freqz
frq, resp = freqz([1,1])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()
plot(frq, angle(resp))
title('Phase Response')
grid()
filtered2 = (sig + 0.5*r_[0,sig[:-1]])
Pxx3, freqs, times, im = specgram(filtered2);
plot(Pxx3[:,10])
[<matplotlib.lines.Line2D at 0x7f0feb0>]
plot(Pxx[:,10])
plot(Pxx2[:,10]/4.0)
plot(Pxx3[:,10]/(1.5**2))
[<matplotlib.lines.Line2D at 0x7f6d4b0>]
from scipy.signal import freqz
frq, resp = freqz([1,1])
semilogy(frq, abs(resp))
frq, resp = freqz([1,0.5])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()
frq, resp = freqz([1,1])
semilogy(frq, abs(resp))
frq, resp = freqz([1,2])
semilogy(frq, abs(resp))
frq, resp = freqz([1,0.5])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()
frq, resp = freqz([1,-1])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()
z-transform on the difference equation:
$$Y(z) = b_0z^{0}X(z) + b_1z^{-1}X(z) + ... + b_Mz^{-M}X(z) - a_1z^{-1}Y(z) - a_2z^{-2}Y(z) - ... - a_Mz^{-M}Y(z)$$$$ [1 + a_1z^{-1} + a_2z^{-2} - ... + a_Mz^{-M}] \cdot Y(z) = [b_0z^{0} + b_1z^{-1} + ... + b_Mz^{-M}]\cdot X(z) $$$$H(z) = \frac{Y(z)}{X(z)} = \frac{b_0z^{0} + b_1z^{-1} + ... + b_Mz^{-M}}{ 1 + a_1z^{-1} - a_2z^{-2} - ... - a_Mz^{-M}}$$from scipy.signal import lfilter
filtered4 = lfilter([1],[1, 1], sig)
specgram(filtered4);
frq, Y = freqz([1], [1,1])
semilogy(frq, abs(Y))
[<matplotlib.lines.Line2D at 0xdde2310>]
frq, Y = freqz([1], [1,0.5])
semilogy(frq, abs(Y))
[<matplotlib.lines.Line2D at 0xe150810>]
frq, Y = freqz([1], [1,-0.5])
semilogy(frq, abs(Y))
[<matplotlib.lines.Line2D at 0xe00f330>]
frq, Y = freqz([1], [1,-0.5])
plot(frq, angle(Y))
[<matplotlib.lines.Line2D at 0x796e7d0>]
from scipy.signal import tf2zpk
tf2zpk([1,2,1],[1])
(array([-1., -1.]), array([], dtype=float64), 1.0)
tf2zpk([1],[1, 1])
(array([], dtype=float64), array([-1.]), 1.0)
tf2zpk([2, 2],[1])
(array([-1.]), array([], dtype=float64), 2.0)
def PoleZeroPlot(b, a):
(zeros,poles,gain) = tf2zpk(b, a)
angle = np.linspace(-np.pi,np.pi,50)
cirx = np.sin(angle)
ciry = np.cos(angle)
figure()
plot(poles.real, poles.imag, 'x', zeros.real, zeros.imag, 'o', cirx,ciry, 'k-')
grid()
xlim((-2, 2))
xlabel('Real')
ylim((-1.5, 1.5))
ylabel('Imag')
gcf().set_figwidth(5)
return (zeros,poles,gain)
PoleZeroPlot([1,-2,1],[1])
(array([ 1., 1.]), array([], dtype=float64), 1.0)
PoleZeroPlot([1],[1,1])
(array([], dtype=float64), array([-1.]), 1.0)
PoleZeroPlot([1],[1,1.1])
(array([], dtype=float64), array([-1.1]), 1.0)
When a pole is outside the unit circle in the z-plane the filter is unstable!
PoleZeroPlot([1],[1.1,1.1])
(array([], dtype=float64), array([-1.]), 0.90909090909090906)
b = [1,0,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726f20f350>]
b = [1,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726f1b1d90>]
b = [1,0,0,0,0,0,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
semilogy(frq, abs(resp))
twinx()
plot(frq, angle(resp), 'r')
[<matplotlib.lines.Line2D at 0x7f726f0df9d0>]
b = [1,0,0,0,0,0,0,-1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
b = [1,0,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726efbc290>]
b = [1,0,0,0,0,0,0,-1]
a = [1,0,0,0,0,0,0,1]
frq, resp = freqz(b,a)
semilogy(frq, abs(resp));
b = [1]
a = [1,0,0,0,0,0,0,1]
frq, resp = freqz(b,a, worN=8192)
semilogy(frq, abs(resp));
b = [1, 1, 1]
a = [1]
frq, resp = freqz(b,a)
semilogy(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726efd4690>]
R = 0.5
theta_c = 1.0
b = [1, -2*R*cos(theta_c), R**2]
a = [1]
frq, resp = freqz(b,a)
semilogy(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726f011c90>]
R = 1
theta_c = 2.0
b = [1, -2*R*cos(theta_c), R**2]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726eaccc50>]
R = 1
theta_c = 1.5
b = [1, -2*R*cos(theta_c), R**2]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726ec3d450>]
PoleZeroPlot(b,a)
(array([ 0.0707372+0.99749499j, 0.0707372-0.99749499j]), array([], dtype=float64), 1.0)
R = 1
theta_c = 2.0
b = [1]
a = [1, -2*R*cos(2.0), R**2]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726e90fd90>]
poles, zeros, k = PoleZeroPlot(b,a)
abs(zeros)
array([ 1., 1.])
R = 0.9
theta_c = 1.5
b = [1]
a = [1, -2*R*cos(theta_c), R**2]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726e9aa2d0>]
poles, zeros, k = PoleZeroPlot(b,a)
Great "cookbook": http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
# low shelf-filter
Fs = 44100
f0 = 10000.0
dBgain = 30.0
S = 1.0 # shelf slope
# -----------------------
A = 10**(dBgain/40)
w0 = 2*pi*f0/Fs
alpha = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 )
b0 = A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha )
b1 = 2*A*( (A-1) - (A+1)*cos(w0) )
b2 = A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha )
a0 = (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha
a1 = -2*( (A-1) + (A+1)*cos(w0) )
a2 = (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha
w, h = freqz([b0, b1, b2], [a0, a1, a2])
semilogy(w,abs(h))
[<matplotlib.lines.Line2D at 0x7f726ed60810>]
Fs = 44100
f0 = 10000.0
dBgain = -20.0
S = 2.0 # shelf slope
A = 10**(dBgain/40)
w0 = 2*pi*f0/Fs
alpha = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 )
b0 = A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha )
b1 = 2*A*( (A-1) - (A+1)*cos(w0) )
b2 = A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha )
a0 = (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha
a1 = -2*( (A-1) + (A+1)*cos(w0) )
a2 = (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha
w, h = freqz([b0, b1, b2], [a0, a1, a2])
plot(w,abs(h))
[<matplotlib.lines.Line2D at 0x7f726e4c6e50>]
shelved = lfilter([b0, b1, b2], [a0, a1, a2], sig)
specgram(shelved);
PoleZeroPlot([b0, b1, b2], [a0, a1, a2]);
http://en.wikipedia.org/wiki/Chebyshev_filter
from scipy.signal import cheby1
b, a = cheby1(2, 0.4, 0.5)
b,a
(array([ 0.38096036, 0.76192073, 0.38096036]), array([ 1. , 0.32450947, 0.27114842]))
w, h = freqz(b,a)
plot(w, np.abs(h), 'b')
[<matplotlib.lines.Line2D at 0x7f726e6cb750>]
PoleZeroPlot(b,a)
(array([-1., -1.]), array([-0.16225474+0.49479472j, -0.16225474-0.49479472j]), 0.3809603640233169)
ripples = [2, 0.5, 0.1, 0.01]
for ripple in ripples:
b, a = cheby1(4, ripple, 0.6)
frq, resp = freqz(b,a)
plot(frq, abs(resp))
legend(ripples)
title('Different ripple values for a Chebyshev I filter')
grid()
orders = [2,3,4,6]
for order in orders:
b, a = cheby1(order, 2, 0.6)
frq, resp = freqz(b,a)
plot(frq, abs(resp))
legend(orders)
title('Different orders for a Chebyshev I filter')
grid()
from scipy.signal import cheby2
ripples = [12, 15, 20, 40]
for ripple in ripples:
b, a = cheby2(4, ripple, 0.6)
frq, resp = freqz(b,a)
plot(frq, abs(resp))
legend(ripples)
title('Different ripple values for a Chebyshev II filter')
grid()
from scipy.signal import iirdesign
Wp = 0.5 # Cutoff frequency
Ws = 0.6 # Stop frequency
Rp = 0.1 # passband maximum loss (gpass)
As = 60 # stoppand min attenuation (gstop)
b,a = iirdesign(Wp, Ws, Rp, As, ftype='butter')
frq, resp = freqz(b,a)
plot(frq, abs(resp))
twinx()
plot(frq, angle(resp), 'r')
title('Butterworth filter')
grid()
Wp = 0.5 # Cutoff frequency
Ws = 0.6 # Stop frequency
Rp = 1 # passband maximum loss (gpass)
As = 20 # stoppand min attenuation (gstop)
types = ['butter', 'ellip', 'cheby1', 'cheby2']
for t in types:
b,a = iirdesign(Wp, Ws, Rp, As, ftype=t)
frq, resp = freqz(b,a)
plot(frq, abs(resp))
legend(types)
title('Filters')
grid()
You can use the iirfilter function from scipy.signal to design and filter in a single step.
http://docs.scipy.org/doc/scipy/reference/signal.html#filter-design
There are several techniques that allow you to arbitrarily define target filter responses.
from scipy.signal import firwin2
freqs = [0.0, 0.5, 1.0]
gains = [1.0, 1.0, 0.0]
order = 150
taps = firwin2(order, freqs, gains)
frq, resp = freqz(taps)
plot(frq, abs(resp))
twinx()
plot(frq, angle(resp), 'r')
[<matplotlib.lines.Line2D at 0x7f726e554b90>]
freqs = [0.0, 0.3, 0.5, 0.8, 1.0]
gains = [1.0, 0.2, 0.5, 0.1, 0.0]
order = 200
taps = firwin2(order, freqs, gains)
frq, resp = freqz(taps)
plot(frq, abs(resp))
twinx()
plot(frq, angle(resp), 'r')
[<matplotlib.lines.Line2D at 0x7f726dd1d2d0>]
shaped = lfilter(taps, [1], sig)
specgram(shaped);
from scipy.signal import remez
freqs = [0, 0.1, 0.2, 0.4, 0.45, 0.5]
gains = [0, 1, 0]
taps = remez(72, freqs, gains, type='bandpass')
frq, resp = freqz(taps)
plot(frq/(2*pi), abs(resp))
twinx()
plot(frq/(2*pi), angle(resp), 'r')
[<matplotlib.lines.Line2D at 0x7f726dbee4d0>]
from scipy.signal import remez
freqs = [0, 0.1, 0.3, 0.4, 0.45, 0.5]
gains = [0, 1, 0]
taps = remez(100, freqs, gains, type='bandpass')
frq, resp = freqz(taps)
plot(frq/(2*pi), abs(resp))
twinx()
plot(frq/(2*pi), angle(resp), 'r')
[<matplotlib.lines.Line2D at 0x7f726daa4150>]
Many ways of doing them. One is:
$$y(n) = -gx(n) + x(n-D) + gy(y-D)$$from scipy.signal import freqz
g = 0.9
b = [-g, 0,0, 1]
a = [1, 0,0,-g]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f726d9515d0>]
plot(frq, angle(resp))
xlabel('freq')
ylabel('angle')
<matplotlib.text.Text at 0x7f726da249d0>
sig = sin(linspace(0, 2*pi*5, 1000, endpoint=False)) + sin(linspace(0, 2*pi*30, 1000, endpoint=False))
plot(sig)
[<matplotlib.lines.Line2D at 0x7f726d83fb50>]
from scipy.signal import lfilter
f = lfilter(b, a, sig)
plot(abs(rfft(f)))
plot(abs(rfft(sig)))
xlim((0,80))
(0, 80)
plot(angle(rfft(f)))
plot(angle(rfft(sig)))
legend(['original', 'all-pass'])
<matplotlib.legend.Legend at 0x7f726d97e750>
stem(angle(rfft(f)),linefmt='b:')
stem(angle(rfft(sig)), 'g.-.')
legend(['original', 'all-pass'])
xlim(0,50)
(0, 50)
angle(resp)[5], angle(resp)[30]
(-1.4376450639128875, -2.7742892239674024)
angle(rfft(f))[5] - angle(rfft(sig))[5]
-1.4165245589590689
angle(rfft(f))[30] - angle(rfft(sig))[30]
3.5012099242089638
hmmm....
angle(rfft(f))[30] - angle(rfft(sig))[30] - (2 * pi)
-2.7819753829706224
plot(f)
plot(sig)
xlim((0, 200))
grid()
By: Andrés Cabrera mantaraya36@gmail.com For MAT course MAT 201A at UCSB
This ipython notebook is licensed under the CC-BY-NC-SA license: http://creativecommons.org/licenses/by-nc-sa/4.0/