%pylab inline
rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division
sig = random.random(2048)*2 -1
Pxx, freqs, times, im = specgram(sig, interpolation='nearest');
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f4e9e71cbd8>
Pxx.shape
(129, 15)
2048/256
8.0
plot(Pxx[:,10])
[<matplotlib.lines.Line2D at 0x7f4e96f60790>]
filtered = (sig + r_[0,sig[:-1]])
Pxx2, freqs, times, im = specgram(filtered);
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f4e9e614a70>
plot(Pxx2[:,10])
[<matplotlib.lines.Line2D at 0x7f4e96dac550>]
plot(Pxx[:,10])
plot(Pxx2[:,10])
[<matplotlib.lines.Line2D at 0x7f4e96d3fc50>]
plot(Pxx[:,10])
plot(Pxx2[:,10]/4.0)
[<matplotlib.lines.Line2D at 0x7f4e96c90dd0>]
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 0x7f4e8b41d690>]
plot(Pxx[:,10])
plot(Pxx2[:,10]/4.0)
plot(Pxx3[:,10]/(1.5**2))
[<matplotlib.lines.Line2D at 0x7f4e8b350710>]
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 0x7f4e8b5510d0>]
frq, Y = freqz([1], [1,0.5])
semilogy(frq, abs(Y))
[<matplotlib.lines.Line2D at 0x7f4e8b60fe90>]
frq, Y = freqz([1], [1,-0.5])
semilogy(frq, abs(Y))
[<matplotlib.lines.Line2D at 0x7f4e8b785e90>]
frq, Y = freqz([1], [1,-0.5])
plot(frq, angle(Y))
[<matplotlib.lines.Line2D at 0x7f4e8b031110>]
from scipy.signal import tf2zpk
tf2zpk([1,1],[1])
(array([-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,1],[1])
(array([-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 0x7f4e8b068d50>]
b = [1,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))
[<matplotlib.lines.Line2D at 0x7f4e8ae4d3d0>]
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 0x7f4e8ad76990>]
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 0x7f4e8aeeb750>]
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 0x7f4e8a77e110>]
R = 0.8
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 0x7f4e8aad5410>]
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 0x7f4e8aa91410>]
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 0x7f4e8a685210>]
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 0x7f4e8a553e10>]
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 0x7f4e8a4cecd0>]
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 0x7f4e8a3ebed0>]
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 0x7f4e8a205510>]
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(4, 0.5, 0.1)
b,a
(array([ 0.00018204, 0.00072816, 0.00109224, 0.00072816, 0.00018204]), array([ 1. , -3.53281335, 4.78185631, -2.93275311, 0.68679537]))
w, h = freqz(b,a)
plot(w, np.abs(h), 'b')
[<matplotlib.lines.Line2D at 0x7f4e89cb83d0>]
PoleZeroPlot(b,a)
(array([-1.00021915+0.j , -0.99999998+0.00021913j, -0.99999998-0.00021913j, -0.99978088+0.j ]), array([ 0.89936989+0.29745863j, 0.89936989-0.29745863j, 0.86703679+0.11665586j, 0.86703679-0.11665586j]), 0.0001820402860029431)
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 0x7f4e89907810>]
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 0x7f4e899b72d0>]
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(50, 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 0x7f4e88cc7ad0>]
from scipy.signal import remez
freqs = [0, 0.1, 0.3, 0.4, 0.45, 0.5]
gains = [0, 1, 0]
taps = remez(16, 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 0x7f4e88775c10>]
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 0x7f4e88616c90>]
plot(frq, angle(resp))
xlabel('freq')
ylabel('angle')
<matplotlib.text.Text at 0x7f4e886c8610>
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 0x7f4e885098d0>]
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 0x7f4e8836c950>
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(sig)
plot(f)
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/