In [1]:
%pylab inline
import seaborn as sb
In [18]:
# Plot exponential signals
t = np.arange(-5., 5., 0.001)
figure(figsize=(17, 4))
subplot(131)
a = 1.1
plot(t, np.exp(a * t), label="$e^{{{0}t}}$".format(a))
ylim(0, 20.)
xlim(t[0], t[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("Exponentially growing signal", fontsize=20)
legend(loc=0, prop={'size':25})
subplot(132)
a = 0.
plot(t, np.exp(a * t), label="$e^{{{0}t}}$".format(a))
ylim(0, 1.5)
xlim(t[0], t[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("DC signal", fontsize=20)
legend(loc=0, prop={'size':25})
subplot(133)
a = -0.5
plot(t, np.exp(a * t), label="$e^{{{0}t}}$".format(a))
ylim(0, 20.)
xlim(t[0], t[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("Exponentially decaying signal", fontsize=20)
legend(loc=0, prop={'size':25})
savefig("img/exp.svg", format="svg");
In [12]:
# Plot exponential signals
n = np.arange(-10, 10, 1)
figure(figsize=(17, 6))
subplot(221)
a = 1.74
stem(n, np.power(a, n), label="${0}^n$".format(a))
ylim(0, 20.)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("Exponentially growing signal", fontsize=20)
legend(loc=0, prop={'size':25})
subplot(222)
a = 0.76
stem(n, np.power(a, n), label="${0}^n$".format(a))
ylim(0, 20.)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("Exponentially growing signal", fontsize=20)
legend(loc=0, prop={'size':25})
subplot(223)
a = -1.56
stem(n, np.power(a, n), label="${0}^n$".format(a))
ylim(-4.5, 4.5)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
legend(loc=0, prop={'size':25})
subplot(224)
a = -0.56
stem(n, np.power(a, n), label="${0}^n$".format(a))
ylim(-3.0, 3.0)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
legend(loc=0, prop={'size':25})
tight_layout()
savefig("img/disc_exp.svg", format="svg");
In [11]:
# Sinusoids
t = np.arange(-2.0, 2.0, 0.001)
s1 = np.sin(2 * np.pi * 1 * t)
s2 = 2 * np.sin(2 * np.pi * 1 * t)
s3 = np.sin(2 * np.pi * 2 * t)
s4 = np.sin(2 * np.pi * 1 * t + 0.25 * np.pi)
figure(figsize=(17, 6))
subplot(221)
plot(t, s1)
ylim(-2.1, 2.1)
xlim(t[0], t[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$\sin(2\pi t)$", fontsize=25)
subplot(222)
plot(t, s2)
ylim(-2.1, 2.1)
xlim(t[0], t[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$2\sin(2\pi t)$", fontsize=25)
subplot(223)
plot(t, s3)
ylim(-2.1, 2.1)
xlim(t[0], t[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$\sin(4\pi t)$", fontsize=25)
subplot(224)
plot(t, s4)
ylim(-2.1, 2.1)
xlim(t[0], t[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$\sin(2\pi t + \pi/4)$", fontsize=25)
tight_layout()
savefig("img/sinu.svg", format="svg");
In [21]:
# Discrete sinusoid
n = np.arange(-10, 10.0, 1.0)
s1 = np.cos(0.3 * np.pi * n)
s2 = 2 * np.cos(0.3 * np.pi * n)
s3 = np.cos(np.pi * n)
s4 = np.cos(0.3 * np.pi * n + 0.25 * np.pi)
figure(figsize=(17, 6))
subplot(221)
stem(n, s1)
ylim(-2.1, 2.1)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$\cos(0.3\pi n)$", fontsize=25)
subplot(222)
stem(n, s2)
ylim(-2.1, 2.1)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$2\cos(0.3\pi n)$", fontsize=25)
subplot(223)
stem(n, s3)
ylim(-2.1, 2.1)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$\cos(\pi n)$", fontsize=25)
subplot(224)
stem(n, s4)
ylim(-2.1, 2.1)
xlim(n[0], n[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title("$\cos(0.3\pi n + \pi/4)$", fontsize=25)
tight_layout()
savefig("img/disc_sinu.svg", format="svg");
In [42]:
# Exponential sinusoids
t = np.arange(-2.0, 2.0, 0.01)
figure(figsize=(17,8))
# sinuisoid parameters
A = 1.4
f = 2. # omega = 2 * pi * f
p = 6. # phi = pi / p
lbl_txt = "${1}e^{{ {0}t }}\sin({2}\pi t + \pi / {2})$"
# exponential parameter
b = 0.67
subplot(221)
x = np.exp(b * t) * A * np.sin(2 * np.pi * f * t + np.pi / p)
plot(t, x, label=lbl_txt.format(b, A, 2*f, p))
plot(t, np.exp(b * t) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
plot(t, -np.exp(b * t) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
ylim(-6, 9)
xticks(fontsize=18)
yticks(fontsize=18)
legend(prop={'size':25})
# exponential parameter
b = -1.02
subplot(222)
x = np.exp(b * t) * A * np.sin(2 * np.pi * f * t + np.pi / p)
plot(t, x, label=lbl_txt.format(b, A, 2*f, p))
plot(t, np.exp(b * t) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
plot(t, -np.exp(b * t) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
ylim(-10, 11)
xticks(fontsize=18)
yticks(fontsize=18)
legend(prop={'size':25});
# exponential parameter
lbl_txt = "${1} * {0}^n\sin({2}\pi n)$"
n = np.arange(-10, 10, 1)
b = 1.17
f = 0.15
subplot(223)
x = np.power(b, n) * A * np.sin(f * np.pi * n)
stem(n, x, label=lbl_txt.format(b, A, f))
plot(n, np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
plot(n, -np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
ylim(-6, 9)
xlim(n[0], n[-1])
xticks(fontsize=18)
yticks(fontsize=18)
legend(prop={'size':25})
n = np.arange(-10, 10, 1)
b = 0.851
f = 0.15
subplot(224)
x = np.power(b, n) * A * np.sin(f * np.pi * n)
stem(n, x, label=lbl_txt.format(b, A, f))
plot(n, np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
plot(n, -np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
ylim(-6, 9)
xlim(n[0], n[-1])
xticks(fontsize=18)
yticks(fontsize=18)
legend(prop={'size':25})
savefig("img/exp_sin.svg", format="svg");
In [4]:
# Exponential sinusoids
t = np.arange(-2.0, 2.0, 0.01)
figure(figsize=(17,8))
A = 1.4
# exponential parameter
lbl_txt = "${1} * {0}^n\sin({2}\pi n)$"
n = np.arange(-10, 10, 1)
b = 1.17
f = 0.15
subplot(121)
x = np.power(b, n) * A * np.sin(f * np.pi * n)
stem(n, x, label=lbl_txt.format(b, A, f))
plot(n, np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
plot(n, -np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
ylim(-6, 9)
xlim(n[0], n[-1])
xticks(fontsize=18)
yticks(fontsize=18)
legend(prop={'size':25})
n = np.arange(-10, 10, 1)
b = 0.851
f = 0.15
subplot(122)
x = np.power(b, n) * A * np.sin(f * np.pi * n)
stem(n, x, label=lbl_txt.format(b, A, f))
plot(n, np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
plot(n, -np.power(b, n) * A, color=sb.xkcd_rgb["pale red"], linestyle='--')
ylim(-6, 9)
xlim(n[0], n[-1])
xticks(fontsize=18)
yticks(fontsize=18)
legend(prop={'size':25})
savefig("img/disc_exp_sin.svg", format="svg");
In [55]:
# Impusle function demo
# Ordinary function.
def g_value(t):
return 13.0 + 5.25 * np.sin(t - np.pi/2.5)
# Impulse function sequence
def f_value(t, n):
return n if (t <= (1./(2*n)) and t >= -1./(2*n)) else 0.
dt = 0.001
time = np.arange(-1., 4.0, dt)
g = [g_value(t) for t in time]
# f_n for different values of N.
N = [2, 4, 8, 16, 32, 64]
f_times = [0] * len(N)
f_vals = [0] * len(N)
g_vals = [0] * len(N)
for i, n in enumerate(N):
dt = 1./(50*n)
f_times[i] = np.arange(-1.5*(1./(2*n)), 1.5*(1./(2*n)), dt)
f_vals[i] = np.array([f_value(t, n) for t in f_times[i]])
g_vals[i] = np.dot(f_vals[i], np.array([g_value(t) for t in f_times[i]])) * dt
figure(figsize=(17,6))
subplot2grid((1,6), (0,0), rowspan=1, colspan=4)
plot(time, g, label="$g$")
for i, n in enumerate(N):
step(f_times[i], f_vals[i], label="$f_{{{0}}}$".format(n))
xlabel('Time', fontsize=25)
xticks(fontsize=20)
yticks(fontsize=20)
ylim(-0.5, 40)
legend(prop={'size':25});
# this is an inset axes over the main axes
subplot2grid((1,6), (0,4), rowspan=1, colspan=2)
plot(N, g_vals, label="$g_n$")
plot(N, [g_value(0)] * len(N), '--', label="$g(0)$")
xlim(N[0], N[-1])
xticks(fontsize=20)
yticks(fontsize=20)
title('$g_n$ as a function of n', fontsize=25)
xlabel("n", fontsize=25)
legend(prop={'size':25})
tight_layout();
savefig("img/impulse_demo.svg", format="svg");
In [60]:
# Discrete impulse
n = np.arange(-10, 10, 1)
imp = 1.0 * (n == 0)
figure(figsize=(14, 3))
stem(n, imp, label="$\delta[n]$")
ylim(-0.1, 1.1)
xlim(n[0], n[-1])
yticks(fontsize=20)
xticks(fontsize=20)
xlabel("Time", fontsize=20)
legend(prop={'size': 25})
savefig("img/disc_imp.svg", format="svg");
In [9]:
# Unit step function
t = np.arange(-10, 10, 0.001)
n = np.arange(-10, 10, 1.0)
u = 1.0 * (t >= 0)
ud = 1.0 * (n >= 0)
figure(figsize=(14, 3))
stem(n, ud, markerfmt="s", label="$u[n]$")
step(t, u, 'k', lw=1, label="$u(t)$")
ylim(-0.1, 1.1)
xlim(n[0], n[-1])
yticks(fontsize=20)
xticks(fontsize=20)
xlabel("Time", fontsize=20)
legend(loc=2, prop={'size': 25})
savefig("img/step.svg", format="svg");
In [34]:
# Signals as vectors
x = [1., 2.0]
y = [-1., 1.0]
figure(figsize=(7, 7))
plot(x[0], x[1], 'bo')
plot(y[0], y[1], 'bo')
arrow(0, 0, x[0]-0.15, x[1]-0.3, head_width=0.1, head_length=0.3, fc='k', ec='k')
arrow(0, 0, y[0]+0.17, y[1]-0.19, head_width=0.1, head_length=0.3, fc='k', ec='k')
plot([0, 0], [-3, 3], '0.5')
plot([-3, 3], [0, 0], '0.5')
text(x[0] - 0.1, x[1] + 0.2, "$x[n] = \{1, 2\}$", fontsize=25)
text(y[0]-1.5, y[1] + 0.25, "$y[n] = \{-1, 1\}$", fontsize=25)
xlim(-3,3)
ylim(-3,3)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel("$x[0]$", fontsize=25)
ylabel("$x[1]$", fontsize=25)
savefig("img/2dvec.svg", format="svg");
In [13]:
# Signals represented through other signals
x = [1., 2.0]
x0 = [1., 1.0]
x1 = [-1., 1.0]
figure(figsize=(7, 6.5))
# plot(x[0], x[1], 'ro')
arrow(0, 0, x0[0]-0.2, x0[1]-0.2, head_width=0.1, head_length=0.3, fc='k', ec='k')
arrow(0, 0, x1[0]+0.2, x1[1]-0.2, head_width=0.1, head_length=0.3, fc='k', ec='k')
arrow(0, 0, x[0]-0.13, x[1]-0.25, head_width=0.1, head_length=0.3, fc='r', ec='r')
text(x[0] - 0.1, x[1] + 0.2, "$x[n] = \{1, 2\}$", fontsize=25)
text(x0[0]-0.2, x0[1] + 0.15, "$x_0[n] = \{1, 1\}$", fontsize=25)
text(x1[0]-1.75, x1[1] + 0.15, "$x_1[n] = \{-1, 1\}$", fontsize=25)
xlim(-3,3)
ylim(-3,3)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel("$x[0]$", fontsize=25)
ylabel("$x[1]$", fontsize=25)
plot([0, 0], [-3, 3], '0.5')
plot([-3, 3], [0, 0], '0.5')
savefig("img/bases.svg", format="svg");
In [14]:
# Time shifting and time scaling.
def x(t):
return (0. if np.abs(t) > 1 else
1 if t > 0. and t <= 1.0 else
t + 1.0)
figure(figsize=(14, 4))
time = np.arange(-5, 5, 0.01)
plot(time, [x(t) for t in time], label="$x(t)$")
ylim(-0.1, 1.1)
xlim(-5.0, 5.0)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=20)
title('Demonstration of time shifting', fontsize=20)
legend(fontsize=20)
savefig("img/tshift0.svg", format="svg");
figure(figsize=(14, 4))
time = np.arange(-5, 5, 0.01)
plot(time, [x(t) for t in time], label="$x(t)$")
tau = 1.5
plot(time, [x(t-tau) for t in time], label="$x(t-{0})$".format(tau))
ylim(-0.1, 1.1)
xlim(-5.0, 5.0)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=20)
title('Demonstration of time shifting', fontsize=20)
legend(fontsize=20)
savefig("img/tshift1.svg", format="svg");
figure(figsize=(14, 4))
time = np.arange(-5, 5, 0.01)
plot(time, [x(t) for t in time], label="$x(t)$")
tau = 1.5
plot(time, [x(t-tau) for t in time], label="$x(t-{0})$".format(tau))
tau = 2.5
plot(time, [x(t+tau) for t in time], label="$x(t+{0})$".format(tau))
ylim(-0.1, 1.1)
xlim(-5.0, 5.0)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=20)
title('Demonstration of time shifting', fontsize=20)
legend(fontsize=20)
savefig("img/tshift2.svg", format="svg");
In [17]:
# Time scaling.
def x(t):
return (0. if np.abs(t) > 1 else
1 if t > 0. and t <= 1.0 else
t + 1.0)
figure(figsize=(14, 4))
time = np.arange(-5, 5, 0.01)
plot(time, [x(t) for t in time], label="$x(t)$")
ylim(-0.1, 1.1)
xlim(-5.0, 5.0)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=20)
title('Demonstration of time scaling', fontsize=20)
legend(fontsize=20)
savefig("img/tscale0.svg", format="svg");
figure(figsize=(14, 4))
time = np.arange(-5, 5, 0.01)
plot(time, [x(t) for t in time], label="$x(t)$")
a = 0.5
plot(time, [x(a*t) for t in time], label="$x({0}t)$".format(a))
ylim(-0.1, 1.1)
xlim(-5.0, 5.0)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=20)
title('Demonstration of time scaling', fontsize=20)
legend(fontsize=20)
savefig("img/tscale1.svg", format="svg");
figure(figsize=(14, 4))
time = np.arange(-5, 5, 0.01)
plot(time, [x(t) for t in time], label="$x(t)$")
a = 0.5
plot(time, [x(a*t) for t in time], label="$x({0}t)$".format(a))
a = 2.5
plot(time, [x(a*t) for t in time], label="$x({0}t)$".format(a))
ylim(-0.1, 1.1)
xlim(-5.0, 5.0)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=20)
title('Demonstration of time scaling', fontsize=20)
legend(fontsize=20)
savefig("img/tscale2.svg", format="svg");
figure(figsize=(14, 4))
time = np.arange(-5, 5, 0.01)
plot(time, [x(t) for t in time], label="$x(t)$")
a = -1.5
plot(time, [x(a*t) for t in time], label="$x({0}t)$".format(a))
ylim(-0.1, 1.1)
xlim(-5.0, 5.0)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=20)
title('Demonstration of time scaling', fontsize=20)
legend(fontsize=20)
savefig("img/tscale3.svg", format="svg");
In [26]:
# Impulse response of a simple RC circuit.
def response(time, T, RC):
v = np.array([0 if t < 0. else
(1./T) * (1 - np.exp(-t/RC)) if t >= 0. and t <= T else
(1./T) * (1 - np.exp(-T/RC)) * np.exp(-(t-T)/RC)
for t in time])
return v
def vin(time, T):
return (1./T) * (time >= 0) * (time <= T)
T = [0.001, 0.25, 0.5, 0.75, 1.0]
R, C = 1.0, 1.0
time = np.arange(-3, 7, 0.001)
figure(figsize=(14, 5))
for _T in T:
plot(time, response(time, _T, R*C), label="$T={0}$".format(_T))
plot(time, vin(time, _T), '0.4', lw=0.5)
legend(fontsize=20)
xlim(time[0], time[-1])
ylim(-0.1, 1.25)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time ($t$)', fontsize=20)
ylabel('$v_{out}(t)$', fontsize=20)
savefig("img/simple_rc.svg", format="svg");
In [48]:
# Impulse respons mechanics demonstration.
def h(t):
return np.exp(t) if t <= 0 else np.exp(-0.25*t)
def x(t):
return 1 if t >= 0. else 0.
time = np.arange(-10, 10, 0.001)
_x = np.array([x(_t) for _t in time])
t0 = 2.5
_h = np.array([h(t0 - _t) for _t in time])
figure(figsize=(14, 5))
ax = subplot(1,1,1)
plot(time, _x, label=r'$x(\tau)$')
plot(time[time < t0], _h[time < t0], color=[0, 0.7, 0], lw=3, label=r'$h(2.5-\tau), \tau < 2.5$')
plot(time[time > t0], _h[time > t0], color=[0.7, 0, 0], lw=3, label=r'$h(2.5-\tau), \tau > 2.5$')
plot(time[time == t0], _h[time == t0], 'ko', lw=3)
ax.fill_between(time[time < t0], 0, _h[time < t0], color=[0, 0.7, 0], alpha=0.5)
ax.fill_between(time[time > t0], 0, _h[time > t0], color=[0.7, 0, 0], alpha=0.5)
xlim(time[0], time[-1])
ylim(-0.1, 1.1)
xlabel(r'$\tau$', fontsize=30)
yticks(fontsize=20)
xticks(fontsize=20)
legend(loc=2, fontsize=20)
savefig("img/imp_resp_mech.svg", format="svg");
In [24]:
# Rectancular pulse and its Fourier transform
def rect_pulse(t0, t):
return (1/np.sqrt(t0)) * (np.abs(t) <= 0.5 * t0)
def my_sinc(t0, omega):
return np.sqrt(t0) * np.sinc(0.5 * t0 * omega)
t = np.arange(-5, 5, 0.001)
omega = np.arange(-2 * np.pi * 7.5, 2 * np.pi * 7.5, 0.01)
t0 = 1.0
figure(figsize=(16, 3))
subplot(121)
plot(t, rect_pulse(t0, t), label="Rectangular\npulse")
xlim(t[0], t[-1])
ylim(-0.1, 1.1)
xlabel("Time (sec)", fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
legend(loc=2, fontsize=18);
subplot(122)
plot(omega, my_sinc(t0, omega), label="sinc")
xlim(omega[0], omega[-1])
ylim(-0.3, 1.1)
xlabel("Frequency ($\\omega$)", fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
legend(fontsize=18);
savefig("img/ft_rect_sinc.svg", format="svg");
In [35]:
# Triangular pulse
def tri_pulse(t):
return (np.abs(t) <= 1.0) * (1 - np.abs(t))
def my_sinc2(omega):
return np.power(np.sinc(0.5 * omega), 2)
t = np.arange(-5, 5, 0.001)
omega = np.arange(-2 * np.pi * 1.5, 2 * np.pi * 1.5, 0.01)
figure(figsize=(16, 3))
subplot(121)
plot(t, tri_pulse(t), label="Triangular\npulse")
xlim(t[0], t[-1])
ylim(-0.1, 1.1)
xlabel("Time (sec)", fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
legend(loc=2, fontsize=18);
subplot(122)
plot(omega, my_sinc2(omega), label="sinc$^2$")
xlim(omega[0], omega[-1])
ylim(-0.01, 1.01)
xlabel("Frequency ($\\omega$)", fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
legend(fontsize=18);
savefig("img/ft_tri_sinc2.svg", format="svg");
In [31]:
# Fourier series - Gibbs Phenomenon.
def sq_wave(T, t):
return 1.0 * (np.abs(t) <= 0.25 * T)
def fs_partial_sum(N, T, t):
x = 0.5 * np.ones(len(t))
for i in xrange(1, N):
x += np.sinc(0.5 * i) * np.cos(2 * np.pi * i * t / T)
return x
T = 1
t = np.arange(-0.5, 0.5, 0.001)
figure(figsize=(16,6))
plot(t, fs_partial_sum(2, T, t), lw=1, label="$\\hat{x}_{2}(t)$")
plot(t, fs_partial_sum(4, T, t), lw=1, label="$\\hat{x}_{4}(t)$")
plot(t, fs_partial_sum(8, T, t), lw=1, label="$\\hat{x}_{8}(t)$")
plot(t, fs_partial_sum(16, T, t), lw=1, label="$\\hat{x}_{16}(t)$")
plot(t, fs_partial_sum(32, T, t), lw=1, label="$\\hat{x}_{32}(t)$")
plot(t, fs_partial_sum(256, T, t), 'k', lw=1, label="$\\hat{x}_{16}(t)$")
plot(t, sq_wave(T, t), label="$x(t)$")
xlim(t[0], t[-1]+0.2)
ylim(-0.2, 1.2)
xlabel('Time (s)', fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
legend(fontsize=20)
title('Partial reconstruction from Fourier series coefficients', fontsize=20)
savefig("img/gibbs_fs.svg", format="svg");
In [54]:
# Complex exponentials
t = np.arange(-5.0, 5.0, 0.001)
figure(figsize=(17,8))
# s is real.
subplot2grid((2, 2), (0, 0), colspan=1)
sig1 = 1.3
sig2 = -0.4
plot(t, np.exp(sig1 * t), label=r'$e^{{{0}t}}$'.format(sig1))
plot(t, np.exp(sig2 * t), label=r'$e^{{{0}t}}$'.format(sig2), color=sb.xkcd_rgb["pale red"])
xlim(t[0], t[-1])
ylim(-0.5, 10)
xticks(fontsize=20)
yticks(fontsize=20)
title(r'$s\,\mathrm{is\,real}$', fontsize=25)
legend(prop={'size':25})
# s is imaginary.
subplot2grid((2, 2), (0, 1), colspan=1)
sig = 0.0
f = 0.5
plot(t, np.exp(sig * t) * np.cos(2 * np.pi * f * t),
label=r'$\Re \,e^{{j2\pi({0})t}}$'.format(f))
plot(t, np.exp(sig * t) * np.sin(2 * np.pi * f * t),
label=r'$\Im \,e^{{j2\pi({0})t}}$'.format(f),
color=sb.xkcd_rgb["pale red"])
xlim(t[0], t[-1])
ylim(-1.1, 3.1)
xticks(fontsize=20)
yticks(fontsize=20)
title(r'$s\,\mathrm{is\,imaginary}$', fontsize=25)
legend(prop={'size':25})
# s is complex.
subplot2grid((2, 2), (1, 0), colspan=1)
sig = 0.53
f = 1.5
plot(t, np.exp(sig * t) * np.cos(2 * np.pi * f * t),
label=r'$\Re \,e^{{{0} + j2\pi({1})t}}$'.format(sig, f))
plot(t, np.exp(sig * t) * np.sin(2 * np.pi * f * t),
label=r'$\Im \,e^{{{0} + j2\pi({1})t}}$'.format(sig, f),
color=sb.xkcd_rgb["pale red"])
xlim(t[0], t[-1])
ylim(-4.1, 4.1)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=18)
title(r'$s\,\mathrm{is\,complex}$', fontsize=25)
legend(loc=2, prop={'size':25})
# s is imaginary.
subplot2grid((2, 2), (1, 1), colspan=1)
sig = -0.41
f = 3.0
plot(t, np.exp(sig * t) * np.cos(2 * np.pi * f * t),
label=r'$\Re \,e^{{{0} + j2\pi({1})t}}$'.format(sig, f))
plot(t, np.exp(sig * t) * np.sin(2 * np.pi * f * t),
label=r'$\Im \,e^{{{0} + j2\pi({1})t}}$'.format(sig, f),
color=sb.xkcd_rgb["pale red"])
xlim(t[0], t[-1])
ylim(-4.1, 4.1)
xticks(fontsize=20)
yticks(fontsize=20)
xlabel('Time (sec)', fontsize=18)
title(r'$s\,\mathrm{is\,complex}$', fontsize=25)
legend(loc=1, prop={'size':25})
tight_layout()
savefig("img/complex_exp.svg", format="svg", dpi=600)
savefig("img/complex_exp.eps", format="eps", dpi=600);
In [49]:
# S plane, poles and zeros plot
z1 = (1, 0)
z2 = (-2, 0)
p1 = (-0.75, 0)
p2 = (-0.5, np.sqrt(3)/2.)
p3 = (-0.5, -np.sqrt(3)/2.)
fig = figure(figsize=(10,10))
plot([-6, 6], [0, 0], 'k')
plot([0, 0], [-6, 6], 'k')
# zeros plots
plot(z1[0], z1[1], 'ko', ms=20, mfc='none')
plot(z2[0], z2[1], 'ko', ms=20, mfc='none')
# poles plots
plot(p1[0], p1[1], 'kx', ms=20, mfc='none')
plot(p2[0], p2[1], 'kx', ms=20, mfc='none')
plot(p3[0], p3[1], 'kx', ms=20, mfc='none')
title(r'$s$ plane', fontsize=25)
xlabel(r'$\Re \{s\} = \sigma$', fontsize=30)
ylabel(r'$\Im \{s\} = \omega$', fontsize=30)
xticks(np.arange(-2.5, 3.0, 0.5), fontsize=16)
yticks(np.arange(-2.5, 3.0, 0.5), fontsize=16)
grid(True)
axis('equal')
xlim(-2.5, 2.5)
ylim(-2.5, 2.5)
savefig("img/splane.svg", format="svg", dpi=600)
savefig("img/splane.eps", format="eps", dpi=600);
In [42]:
# S plane, poles and zeros plot
z1 = (1, 0)
p1 = (-0.75, 0)
w = 1.5
fig = figure(figsize=(10.5,10.5))
plot([-6, 6], [0, 0], '0.5')
plot([0, 0], [-6, 6], '0.5')
# zeros plots
plot(z1[0], z1[1], 'ko', ms=20, mfc='none')
# poles plots
plot(p1[0], p1[1], 'kx', ms=20, mfc='none')
# Pole and zero arrows to omega.
sz = 0.2
ang = np.arctan2((w - z1[1]), (- z1[0]))
arrow(z1[0], z1[1],
0 - z1[0] - sz * np.cos(ang),
w - z1[1] - sz * np.sin(ang),
head_width=0.1, head_length=sz, fc='k', ec='k')
# Indicate angle.
_t = np.arange(0, ang, 0.01)
_x, _y = 0.25 * np.cos(_t), .25 * np.sin(_t)
plot(_x + z1[0], _y + z1[1], 'b')
ang = np.arctan2((w - p1[1]), (- p1[0]))
arrow(p1[0], p1[1],
0 - p1[0] - sz * np.cos(ang),
w - p1[1] - sz * np.sin(ang),
head_width=0.1, head_length=sz, fc='k', ec='k')
# Indicate angle.
_t = np.arange(0, ang, 0.01)
_x, _y = 0.25 * np.cos(_t), .25 * np.sin(_t)
plot(_x + p1[0], _y + p1[1], 'b')
# Text on the plot.
text(z1[0], -0.3, r'$z_1$', fontsize=25)
text(p1[0]-0.2, -0.3, r'$p_1$', fontsize=25)
text(0-0.3, w, r'$j\omega$', fontsize=25)
text(0.4, 1.0, r'$\left|j\omega - z_1\right|$', fontsize=25)
text(-0.9, 1.0, r'$\left|j\omega - p_1\right|$', fontsize=25)
text(z1[0]+0.15, z1[1]+0.25, r'$\angle (j\omega - z_1)$', fontsize=25)
text(p1[0]+0.25, p1[1]+0.15, r'$\angle (j\omega - p_1)$', fontsize=25)
text(z1[0]-0.2, -0.2, r'$A$', fontsize=25)
text(p1[0], -0.2, r'$B$', fontsize=25)
text(0.05, w, r'$C$', fontsize=25)
text(-1.75, 1.5, r'$H(s) = \frac{s+z_1}{s+p_1}$', fontsize=25)
text(-1.75, -1.0,
r'$\left|H(\omega)\right| = \frac{\left|j\omega+z_1\right|}{\left|j\omega+p_1\right|}=\frac{\left|\overline{AC}\right|}{\left|\overline{BC}\right|}$', fontsize=25)
text(-1.75, -1.5, r'$\angle H(\omega) = \angle (j\omega+z_1) - \angle (j\omega+p_1) = (\pi - \angle BAC) - \angle ABC$', fontsize=25)
title(r'$s$ plane', fontsize=25)
xlabel(r'$\Re \{s\} = \sigma$', fontsize=30)
ylabel(r'$\Im \{s\} = \omega$', fontsize=30)
xticks(np.arange(-2.5, 3.0, 0.5), fontsize=16)
yticks(np.arange(-2.5, 3.0, 0.5), fontsize=16)
grid(True)
axis('equal')
xlim(-2, 2)
ylim(-2, 2)
tight_layout()
savefig("img/geo_freq_resp.svg", format="svg", dpi=600)
savefig("img/geo_freq_resp.eps", format="eps", dpi=600);
In [34]:
# Frequency response of first order system
w = np.arange(0, 2 * np.pi * 100, 0.001)
# Low pass
p = 2 * np.pi * 0.5
mag = np.array([p / np.sqrt(_w**2 + p**2) for _w in w])
ph = np.array([-np.arctan2(_w, p) for _w in w])
figure(figsize=(16, 8))
subplot(2,2,1)
plot(w / (2 * np.pi), mag)
ylabel('Magnitude', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.1 * w[-1] / (2 * np.pi))
xticks(fontsize=16)
yticks(fontsize=16)
text(5.0, 0.7, r'$H(s) = \frac{{{0:0.3f}}}{{s+{0:0.3f}}}$'.format(p), fontsize=25)
subplot(2,2,2)
semilogx(w / (2 * np.pi), 20 * np.log10(mag))
ylabel('Magnitude (dB)', fontsize=18)
xlim(w[0] / (2 * np.pi), w[-1] / (2 * np.pi))
ylim(-50, 1.0)
xticks(fontsize=16)
yticks(fontsize=16)
subplot(2,2,3)
plot(w / (2 * np.pi), ph)
xlabel('Frequency (Hz)', fontsize=18)
ylabel('Phase', fontsize=18)
xticks(fontsize=16)
yticks(fontsize=16)
xlim(w[0] / (2 * np.pi), 0.1 * w[-1] / (2 * np.pi))
subplot(2,2,4)
semilogx(w / (2 * np.pi), ph)
xlabel('Frequency (Hz)', fontsize=18)
ylabel('Phase', fontsize=18)
xticks(fontsize=16)
yticks(fontsize=16)
xlim(w[0] / (2 * np.pi), w[-1] / (2 * np.pi))
tight_layout()
savefig("img/1st_sys.svg", format="svg", dpi=600)
savefig("img/1st_sys.eps", format="eps", dpi=600);
In [53]:
def mag_resp(wn, Q, w):
"""Returns the magnitude response of the 2nd order system."""
return np.array([(wn**2) / np.sqrt((wn**2 - _w**2)**2 + (_w*wn/Q)**2) for _w in w])
def ph_resp(wn, Q, w):
"""Returns the phase response of the 2nd order system."""
return np.array([-np.arctan2(_w*wn/Q, (wn**2 - _w**2)) for _w in w])
# Frequency response of second order system
w = np.arange(0, 2 * np.pi * 100, 0.001)
# Low pass
wn = 2 * np.pi * 0.5
Q = [0.5, 1.0, 2.0]
mag = [mag_resp(wn, _q, w) for _q in Q]
ph = [ph_resp(wn, _q, w) for _q in Q]
figure(figsize=(16, 8))
subplot(2,2,1)
for i in xrange(len(Q)):
plot(w / (2 * np.pi), mag[i], label=r'$Q={0:0.1f}$'.format(Q[i]))
legend(fontsize=18)
ylabel('Magnitude', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.05 * w[-1] / (2 * np.pi))
xticks(fontsize=16)
yticks(fontsize=16)
text(1.5, 0.7,
r'$H(s) = \frac{{\omega_n^2}}{{s^2+(Q/\omega_n)s+\omega_n^2}}$',
fontsize=25)
subplot(2,2,2)
for i in xrange(len(Q)):
semilogx(w / (2 * np.pi), 20 * np.log10(mag[i]), label=r'$Q={0:0.1f}$'.format(Q[i]))
legend(fontsize=18)
ylabel('Magnitude (dB)', fontsize=18)
xlim(w[0] / (2 * np.pi), w[-1] / (2 * np.pi))
xticks(fontsize=16)
yticks(fontsize=16)
subplot(2,2,3)
for i in xrange(len(Q)):
plot(w / (2 * np.pi), ph[i], label=r'$Q={0:0.1f}$'.format(Q[i]))
legend(fontsize=18)
xlabel('Frequency (Hz)', fontsize=18)
ylabel('Phase', fontsize=18)
xticks(fontsize=16)
yticks(fontsize=16)
ylim(-3.3, 0.2)
xlim(w[0] / (2 * np.pi), 0.05 * w[-1] / (2 * np.pi))
subplot(2,2,4)
for i in xrange(len(Q)):
semilogx(w / (2 * np.pi), ph[i], label=r'$Q={0:0.1f}$'.format(Q[i]))
legend(fontsize=18)
xlabel('Frequency (Hz)', fontsize=18)
ylabel('Phase', fontsize=18)
xticks(fontsize=16)
yticks(fontsize=16)
ylim(-3.3, 0.2)
xlim(w[0] / (2 * np.pi), w[-1] / (2 * np.pi))
tight_layout()
savefig("img/2nd_sys.svg", format="svg", dpi=600)
savefig("img/2nd_sys.eps", format="eps", dpi=600);
In [4]:
def pole_mag(p, w):
"""Returns the magnitude response in dB for a simple pole."""
return 20 * np.log(np.array([p / np.sqrt(_w**2 + p**2) for _w in w]))
def pole_ph(p, w):
"""Returns the phase response in degrees for a simple pole."""
return np.array([-np.arctan2(_w, p) * 180 / np.pi for _w in w])
def pole_mag_bode(p, w):
"""Returns the Bode magnitude response in dB for a simple pole."""
return 20 * np.log(np.array([1 if _w < p else p / _w
for _w in w]))
def pole_ph_bode(p, w):
"""Returns the Bode phase response in degrees for a simple pole."""
delta_w = np.log(10. * p) - np.log(0.1 * p)
return np.array(
[np.arctan2(0, p) if _w < 0.1 * p
else - np.arctan2(_w, 0) if _w > 10 * p
else (- np.arctan2(_w, 0) - np.arctan2(0, p)) * (np.log(_w) - np.log(0.1 * p))/(delta_w)
for _w in w]) * 180 / np.pi
def zero_mag(z, w):
"""Returns the magnitude response in dB for a simple zero."""
return 20 * np.log(np.array([np.sqrt(_w**2 + z**2) / z for _w in w]))
def zero_mag_bode(z, w):
"""Returns the magnitude response in dB for a simple zero."""
return 20 * np.log(np.array([1 if _w < z else _w / z
for _w in w]))
def zero_ph(p, w):
"""Returns the phase response in degrees for a simple zero."""
return np.array([np.arctan2(_w, p) * 180 / np.pi for _w in w])
def zero_ph_bode(p, w):
"""Returns the Bode phase response in degrees for a simple zero."""
delta_w = np.log(10. * p) - np.log(0.1 * p)
return np.array(
[np.arctan2(0, p) if _w < 0.1 * p
else np.arctan2(_w, 0) if _w > 10 * p
else (np.arctan2(_w, 0) - np.arctan2(0, p)) * (np.log(_w) - np.log(0.1 * p))/(delta_w)
for _w in w]) * 180 / np.pi
# BODE PLOT OF POLES, ZEROs.
w = np.arange(0, 2 * np.pi * 100, 0.001)
# simple pole.
p = 2 * np.pi * 0.5
mag = pole_mag(p, w)
mag_b = pole_mag_bode(p, w)
ph = pole_ph(p, w)
ph_b = pole_ph_bode(p, w)
figure(figsize=(16, 8))
subplot(2,2,1)
semilogx(w / (2 * np.pi), mag, label='Actual magnitude response')
semilogx(w / (2 * np.pi), mag_b, label='Bode approximation')
semilogx((p / (2*np.pi)) * np.ones(2), [-80, 5], 'k--', lw=1)
ylabel('Magnitude (dB)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-80, 5)
xticks(fontsize=16)
yticks(fontsize=16)
title('Bode approximation of single simple real pole', fontsize=20)
legend(loc=3, fontsize=16)
subplot(2,2,3)
semilogx(w / (2 * np.pi), ph, label='Actual phase response')
semilogx(w / (2 * np.pi), ph_b, label='Bode approximation')
semilogx((p / (2*np.pi)) * np.ones(2), [-100, 10], 'k--', lw=1)
ylabel('Phase (deg)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-100, 10)
xticks(fontsize=16)
yticks(fontsize=16)
legend(loc=3, fontsize=16)
# simple zero.
z = 2 * np.pi * 1
mag = zero_mag(z, w)
mag_b = zero_mag_bode(z, w)
ph = zero_ph(p, w)
ph_b = zero_ph_bode(p, w)
subplot(2,2,2)
semilogx(w / (2 * np.pi), mag, label='Actual magnitude response')
semilogx(w / (2 * np.pi), mag_b, label='Bode approximation')
semilogx((z / (2*np.pi)) * np.ones(2), [-5, 50], 'k--', lw=1)
ylabel('Magnitude (dB)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-5, 50)
xticks(fontsize=16)
yticks(fontsize=16)
title('Bode approximation of single simple real zero', fontsize=20)
legend(loc=2, fontsize=16)
subplot(2,2,4)
semilogx(w / (2 * np.pi), ph, label='Actual phase response')
semilogx(w / (2 * np.pi), ph_b, label='Bode approximation')
semilogx((p / (2*np.pi)) * np.ones(2), [-10, 100], 'k--', lw=1)
ylabel('Phase (deg)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-10, 100)
xticks(fontsize=16)
yticks(fontsize=16)
legend(loc=2, fontsize=16)
tight_layout()
savefig("img/1st_bode.svg", format="svg", dpi=600)
savefig("img/1st_bode.eps", format="eps", dpi=600);
In [4]:
def pole_mag(wn, Q, w):
"""Returns the magnitude response of 2nd order poles."""
return 20 * np.log(np.array([(wn**2) / np.sqrt((wn**2 - _w**2)**2 + (_w*wn/Q)**2)
for _w in w]))
def pole_ph(wn, Q, w):
"""Returns the phase response of 2nd order poles."""
return np.array([-np.arctan2(_w*wn/Q, (wn**2 - _w**2)) for _w in w])* 180./np.pi
def pole_mag_bode(wn, Q, w):
"""Returns the Bode magnitude response in dB for 2nd order poles."""
return - 20 * np.log(np.array([1 if _w < wn else (_w / wn)**2
for _w in w]))
def pole_ph_bode(wn, Q, w):
"""Returns the Bode phase response in degrees for 2nd order poles."""
delta_w = np.log(10. * wn) - np.log(0.1 * wn)
return np.array(
[0 if _w < 0.1 * wn
else -np.pi if _w > 10 * wn
else - np.pi * (np.log(_w) - np.log(0.1 * wn))/(delta_w)
for _w in w]) * 180 / np.pi
def zero_mag(wn, Q, w):
"""Returns the magnitude response of 2nd order poles."""
return 20 * np.log(np.array([(np.sqrt((wn**2 - _w**2)**2 + (_w*wn/Q)**2) / wn**2)
for _w in w]))
def zero_ph(wn, Q, w):
"""Returns the phase response of 2nd order poles."""
return np.array([np.arctan2(_w*wn/Q, (wn**2 - _w**2)) for _w in w])* 180./np.pi
def zero_mag_bode(wn, Q, w):
"""Returns the Bode magnitude response in dB for 2nd order poles."""
return 20 * np.log(np.array([1 if _w < wn else (_w / wn)**2
for _w in w]))
def zero_ph_bode(wn, Q, w):
"""Returns the Bode phase response in degrees for 2nd order poles."""
delta_w = np.log(10. * wn) - np.log(0.1 * wn)
return np.array(
[0 if _w < 0.1 * wn
else np.pi if _w > 10 * wn
else np.pi * (np.log(_w) - np.log(0.1 * wn))/(delta_w)
for _w in w]) * 180 / np.pi
In [ ]:
# BODE PLOT OF POLES, ZEROs.
w = np.arange(0, 2 * np.pi * 100, 0.001)
# 2nd order poles
wn = 2 * np.pi * 0.5
Q = [0.5, 4.0]
mag = [pole_mag(wn, _q, w) for _q in Q]
mag_b = [pole_mag_bode(wn, _q, w) for _q in Q]
ph = [pole_ph(wn, _q, w) for _q in Q]
ph_b = [pole_ph_bode(wn, _q, w) for _q in Q]
figure(figsize=(16, 8))
subplot(2,2,1)
semilogx(w / (2 * np.pi), mag[0], label='Actual magnitude response (Q={0:0.1f})'.format(Q[0]))
semilogx(w / (2 * np.pi), mag[1], label='Actual magnitude response (Q={0:0.1f})'.format(Q[1]))
semilogx(w / (2 * np.pi), mag_b[0], label='Bode approximation')
semilogx((wn / (2*np.pi)) * np.ones(2), [-160, 40], 'k--', lw=1)
ylabel('Magnitude (dB)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-160, 40)
xticks(fontsize=16)
yticks(fontsize=16)
title('Bode approximation of 2$^{nd}$ order pole.', fontsize=20)
legend(loc=3, fontsize=16)
subplot(2,2,3)
semilogx(w / (2 * np.pi), ph[0], label='Actual phase response (Q={0:0.1f})'.format(Q[0]))
semilogx(w / (2 * np.pi), ph[1], label='Actual phase response (Q={0:0.1f})'.format(Q[1]))
semilogx(w / (2 * np.pi), ph_b[0], label='Bode approximation')
semilogx((wn / (2*np.pi)) * np.ones(2), [-190, 10], 'k--', lw=1)
ylabel('Phase (deg)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-190, 10)
xticks(fontsize=16)
yticks(fontsize=16)
legend(loc=3, fontsize=16)
# 2nd order zero.
wn = 2 * np.pi * 1.0
Q = [0.5, 4.0]
mag = [zero_mag(wn, _q, w) for _q in Q]
mag_b = [zero_mag_bode(wn, _q, w) for _q in Q]
ph = [zero_ph(wn, _q, w) for _q in Q]
ph_b = [zero_ph_bode(wn, _q, w) for _q in Q]
subplot(2,2,2)
semilogx(w / (2 * np.pi), mag[0], label='Actual magnitude response (Q={0:0.1f})'.format(Q[0]))
semilogx(w / (2 * np.pi), mag[1], label='Actual magnitude response (Q={0:0.1f})'.format(Q[1]))
semilogx(w / (2 * np.pi), mag_b[0], label='Bode approximation')
semilogx((wn / (2*np.pi)) * np.ones(2), [-40, 160], 'k--', lw=1)
ylabel('Magnitude (dB)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-40, 160)
xticks(fontsize=16)
yticks(fontsize=16)
title('Bode approximation of 2$^{nd}$ order zero.', fontsize=20)
legend(loc=2, fontsize=16)
subplot(2,2,4)
semilogx(w / (2 * np.pi), ph[0], label='Actual phase response (Q={0:0.1f})'.format(Q[0]))
semilogx(w / (2 * np.pi), ph[1], label='Actual phase response (Q={0:0.1f})'.format(Q[1]))
semilogx(w / (2 * np.pi), ph_b[0], label='Bode approximation')
semilogx((wn / (2*np.pi)) * np.ones(2), [-10, 190], 'k--', lw=1)
ylabel('Phase (deg)', fontsize=18)
xlim(w[0] / (2 * np.pi), 0.5 * w[-1] / (2 * np.pi))
ylim(-10, 190)
xticks(fontsize=16)
yticks(fontsize=16)
legend(loc=2, fontsize=16)
tight_layout()
savefig("img/2nd_bode.svg", format="svg", dpi=600)
savefig("img/2nd_bode.eps", format="eps", dpi=600);
In [4]:
# Arbitrary transfer function
def pole_mag(p, w):
"""Returns the magnitude response in dB for a simple pole."""
return 20 * np.log(np.array([p / np.sqrt(_w**2 + p**2) for _w in w]))
def pole_ph(p, w):
"""Returns the phase response in degrees for a simple pole."""
return np.array([-np.arctan2(_w, p) * 180 / np.pi for _w in w])
def pole_mag_bode(p, w):
"""Returns the Bode magnitude response in dB for a simple pole."""
return 20 * np.log(np.array([1 if _w < p else p / _w
for _w in w]))
def pole_ph_bode(p, w):
"""Returns the Bode phase response in degrees for a simple pole."""
delta_w = np.log(10. * p) - np.log(0.1 * p)
return np.array(
[np.arctan2(0, p) if _w < 0.1 * p
else - np.arctan2(_w, 0) if _w > 10 * p
else (- np.arctan2(_w, 0) - np.arctan2(0, p)) * (np.log(_w) - np.log(0.1 * p))/(delta_w)
for _w in w]) * 180 / np.pi
def zero_mag(z, w):
"""Returns the magnitude response in dB for a simple zero."""
return 20 * np.log(np.array([np.sqrt(_w**2 + z**2) / z for _w in w]))
def zero_mag_bode(z, w):
"""Returns the magnitude response in dB for a simple zero."""
return 20 * np.log(np.array([1 if _w < z else _w / z
for _w in w]))
def zero_ph(p, w):
"""Returns the phase response in degrees for a simple zero."""
return np.array([np.arctan2(_w, p) * 180 / np.pi for _w in w])
def zero_ph_bode(p, w):
"""Returns the Bode phase response in degrees for a simple zero."""
delta_w = np.log(10. * p) - np.log(0.1 * p)
return np.array(
[np.arctan2(0, p) if _w < 0.1 * p
else np.arctan2(_w, 0) if _w > 10 * p
else (np.arctan2(_w, 0) - np.arctan2(0, p)) * (np.log(_w) - np.log(0.1 * p))/(delta_w)
for _w in w]) * 180 / np.pi
# Poles and zeros.
p = 2*np.pi*np.array([10., 500.0, 1000., 7500.])
z = 2*np.pi*np.array([50., 1000., 5000.])
w = np.arange(0, 2*np.pi*100000, 1)
# magnitudes
mag_p = [pole_mag(_p, w) for _p in p]
mag_z = [zero_mag(_z, w) for _z in z]
mag_p_b = [pole_mag_bode(_p, w) for _p in p]
mag_z_b = [zero_mag_bode(_z, w) for _z in z]
# phases
ph_p = [pole_ph(_p, w) for _p in p]
ph_z = [zero_ph(_z, w) for _z in z]
ph_p_b = [pole_ph_bode(_p, w) for _p in p]
ph_z_b = [zero_ph_bode(_z, w) for _z in z]
mag = np.zeros(len(w))
mag_b = np.zeros(len(w))
ph = np.zeros(len(w))
ph_b = np.zeros(len(w))
for i in xrange(len(mag_z)):
mag += mag_z[i]
ph += ph_z[i]
mag_b += mag_z_b[i]
ph_b += ph_z_b[i]
for i in xrange(len(mag_p)):
mag += mag_p[i]
ph += ph_p[i]
mag_b += mag_p_b[i]
ph_b += ph_p_b[i]
In [12]:
figure(figsize=(10, 10))
subplot(211)
for i in xrange(len(mag_z)):
semilogx(w / (2 * np.pi), mag_z[i], '0.5', lw=1)
semilogx(w / (2 * np.pi), mag_z_b[i], 'k', lw=0.5)
for i in xrange(len(mag_p)):
semilogx(w / (2 * np.pi), mag_p[i], '0.5', lw=1)
semilogx(w / (2 * np.pi), mag_p_b[i], 'k', lw=0.5)
semilogx(w / (2 * np.pi), mag, lw=2, label='Actual magnitude response')
semilogx(w / (2 * np.pi), mag_b, lw=2, color=sb.xkcd_rgb["pale red"], label='Bode approximation')
ylabel('Magnitude (dB)', fontsize=20)
xticks(fontsize=16)
yticks(fontsize=16)
legend(loc=2, fontsize=18)
title('Magnitude plot of the transfer function.', fontsize=20)
subplot(212)
for i in xrange(len(ph_z)):
semilogx(w / (2 * np.pi), ph_z[i], '0.5', lw=1)
semilogx(w / (2 * np.pi), ph_z_b[i], 'k', lw=0.5)
for i in xrange(len(mag_p)):
semilogx(w / (2 * np.pi), ph_p[i], '0.5', lw=1)
semilogx(w / (2 * np.pi), ph_p_b[i], 'k', lw=0.5)
semilogx(w / (2 * np.pi), ph, lw=2, label='Actual magnitude response')
semilogx(w / (2 * np.pi), ph_b, lw=2, color=sb.xkcd_rgb["pale red"], label='Bode approximation')
ylabel('Phase (deg)', fontsize=20)
xticks(fontsize=16)
yticks(fontsize=16)
legend(loc=2, fontsize=18)
title('Phase plot of the transfer function.', fontsize=20)
tight_layout()
savefig("img/bode.svg", format="svg", dpi=600)
savefig("img/bode.eps", format="eps", dpi=600);
In [74]:
# Signal and noise spectra.
f1 = np.arange(0.1, 1000, 0.1)
sig = 0.01 * np.cumsum(np.exp(-0.01*(f1-1000)**2))[::-1]
noise = 0.001 * np.cumsum(np.exp(-0.01*(f1-100)**2))
figure(figsize=(8,3))
semilogx(f1, sig, label=r'Signal $\left|X_s(\omega)\right|^2$')
semilogx(f1, noise, label=r'Noise $\left|X_n(\omega)\right|^2$',
color=sb.xkcd_rgb["pale red"])
legend(fontsize=16)
xticks(fontsize=16)
yticks(fontsize=16)
xlabel('Frequency (Hz)', fontsize=16)
ylim(-0.1, 1.1)
tight_layout()
savefig("img/sig_noise_sep.svg", format="svg", dpi=600)
savefig("img/sig_noise_sep.eps", format="eps", dpi=600);
In [75]:
# Signal and noise spectra, with idea filter.
f1 = np.arange(0.1, 1000, 0.1)
sig = 0.01 * np.cumsum(np.exp(-0.01*(f1-1000)**2))[::-1]
noise = 0.001 * np.cumsum(np.exp(-0.01*(f1-100)**2))
filt = 1.0 * (f1 < 40.)
figure(figsize=(8,3))
semilogx(f1, sig, label=r'Signal $\left|X_s(\omega)\right|^2$')
semilogx(f1, noise, label=r'Noise $\left|X_n(\omega)\right|^2$',
color=sb.xkcd_rgb["pale red"])
semilogx(f1, filt, label=r'Ideal filter $H_{ideal}(\omega)$',
color='k', ls='--')
legend(loc=3, fontsize=16)
xticks(fontsize=16)
yticks(fontsize=16)
xlabel('Frequency (Hz)', fontsize=16)
ylim(-0.1, 1.1)
savefig("img/sig_noise_ideal_filt.svg", format="svg", dpi=600)
savefig("img/sig_noise_ideal_filt.eps", format="eps", dpi=600);
In [69]:
# Signal and noise spectra, with idea filter.
f1 = np.arange(0.0, 1000, 0.1)
real_filt_ub = [1.0 if _f < 450.
else 0.1 if _f > 550. else
np.NaN for _f in f1]
real_filt_lb = [0.8 if _f < 450.
else np.NaN for _f in f1]
figure(figsize=(8,5))
plot(f1, real_filt_ub, color='k', lw=2)
plot(f1, real_filt_lb, color='k', lw=2)
plot([-100, 1000], [0., 0.], color='0.5', lw=1)
plot([0.1, 0.1], [-0.2, 1.3], color='0.5', lw=1)
plot([450, 450], [0., 1.], color='0', lw=1, ls='--')
plot([550, 550], [0.1, 0.], color='0', lw=1, ls='--')
plot([0, 1000], [0.1, 0.1], color='0', lw=1, ls='--')
legend(loc=3, fontsize=16)
xticks(fontsize=16)
yticks(fontsize=16)
# xlabel('Frequency (Hz)', fontsize=16)
text(-35, 0.97, r'$1$', fontsize=20)
text(-100, 0.77, r'$1-\epsilon$', fontsize=20)
text(-35, 0.07, r'$\delta$', fontsize=20)
text(435, -0.1, r'$\omega_p$', fontsize=20)
text(535, -0.1, r'$\omega_s$', fontsize=20)
# text(-35, 0.07, r'$\delta$', fontsize=20)
xlim(-100, 1000)
ylim(-0.15, 1.1)
yticks([])
xticks([])
tight_layout()
savefig("img/real_filt.svg", format="svg", dpi=600)
savefig("img/real_filt.eps", format="eps", dpi=600);
In [16]:
# Butterworth filter.
def get_butterworth_response(w, w0, N):
"""Returns the Butterworth filter magnitude response.
"""
return np.array([1. / np.sqrt(1 + np.power((_w/w0), 2*N)) for _w in w])
w = np.arange(0.01, 100, 0.01)
N_orders = [1, 2, 4, 8, 16]
mag = [get_butterworth_response(w, 1, _n) for _n in N_orders]
figure(figsize=(10, 4))
for i in xrange(len(N_orders)):
semilogx(w, mag[i], label=r'$N={0}$'.format(N_orders[i]))
legend(fontsize=20)
title('Butterworth filter magnitude response', fontsize=20)
xlabel('Frequency (Hz)', fontsize=20)
ylabel(r'$|H(\omega)|$', fontsize=20)
xticks(fontsize=18)
yticks(fontsize=18)
xlim(0.01, 100)
ylim(-0.1, 1.1)
tight_layout();
savefig("img/butter.svg", format="svg", dpi=600)
savefig("img/butter.eps", format="eps", dpi=600);
In [13]:
# Butterworth response of different orders.
def get_butterworth_response(w, w0, N):
"""Returns the Butterworth filter magnitude response.
"""
return np.array([1. / np.sqrt(1 + np.power((_w/w0), 2*N)) for _w in w])
def get_butterworth_poles(w0, N):
"""Returns the coordinates of all the poles of the
Butterworth transfer function H(s)H(-s)."""
if N % 2 == 1:
# Odd order
return [[w0 * np.cos(n * np.pi / N),
w0 * np.sin(n * np.pi / N)] for n in xrange(2*N)]
else:
# Even order
return [[w0 * np.cos((2*n + 1) * np.pi / (2 * N)),
w0 * np.sin((2*n + 1) * np.pi / (2 * N))] for n in xrange(2*N)]
In [15]:
#### Butterworth filter magnitude response and its poles.
n_ord = [1, 2, 5, 10]
w = np.arange(0.01, 100., 0.01)
poles = [get_butterworth_poles(1., n) for n in n_ord]
# unit circle.
teta = np.arange(0., 2 * np.pi, 0.01)
_x, _y = np.cos(teta), np.sin(teta)
figure(figsize=(16, 4))
for c in xrange(len(n_ord)):
subplot(1, 4, c+1)
for p in poles[c]:
plot(p[0], p[1], 'x', color='k', ms=10)
plot(_x, _y, '0.7', ls='--')
plot([-1.5, 1.5], [0, 0], '0.8')
plot([0., 0.], [-1.5, 1.5], '0.8')
xlim(-1.5, 1.5)
ylim(-1.5, 1.5)
xticks(fontsize=14)
yticks(fontsize=14)
title('Order: {0}'.format(n_ord[c]), fontsize=18)
tight_layout()
savefig("img/butter_poles.svg", format="svg", dpi=600)
savefig("img/butter_poles.eps", format="eps", dpi=600);
In [160]:
# Design of Butterworth filter.
# Butterworth response of different orders.
def get_butterworth_response(w, w0, N):
"""Returns the Butterworth filter magnitude response.
"""
return np.array([1. / np.sqrt(1 + np.power((_w/w0), 2*N)) for _w in w])
def get_butterworth_poles(w0, N):
"""Returns the coordinates of all the poles of the
Butterworth transfer function H(s)H(-s)."""
if N % 2 == 1:
# Odd order
return [[w0 * np.cos(n * np.pi / N),
w0 * np.sin(n * np.pi / N)] for n in xrange(2*N)]
else:
# Even order
return [[w0 * np.cos((2*n + 1) * np.pi / (2 * N)),
w0 * np.sin((2*n + 1) * np.pi / (2 * N))] for n in xrange(2*N)]
def get_butterworth_param(epsilon, delta, wp, ws):
"""Returns the Butterworth lowpass filter parameters.
"""
_n1 = (1/np.power(epsilon, 2)) - 1
_n2 = (1/np.power(delta, 2)) - 1
_d = wp/ws
N = int(np.ceil(0.5 * np.log(_n1/_n2)/np.log(_d)))
w0 = wp / np.power(_n1, 1./(2*N))
return int(N), w0
# Butterworth filter specifications.
epsilon, delta, wp, ws = 0.9, 0.1, 50*np.pi, 100*np.pi
# Filter parameters
N, w0 = get_butterworth_param(epsilon, delta, wp, ws)
# magnitude response and poles.
w = 2 * np.pi * np.arange(1., 1000., 0.1)
mag = get_butterworth_response(w, w0, N)
poles = get_butterworth_poles(w0, N)
figure(figsize=(12.5, 4.75))
subplot2grid((1,3), (0,0), rowspan=1, colspan=2)
semilogx([w[0], ws], np.ones(2), '0.5', ls='--')
semilogx(w, mag, lw=2)
semilogx([w[0], wp], epsilon*np.ones(2), '0.5', ls='--')
semilogx([ws, w[-1]], delta*np.ones(2), '0.5', ls='--')
semilogx(wp*np.ones(2), [-0.01, epsilon], '0.5', ls='--')
semilogx(ws*np.ones(2), [delta, 1.], '0.5', ls='--')
title(r'Butterworth Filter ($N={0}, \omega_0={1:0.3f}$)'.format(N, w0), fontsize=20)
ylim(-0.01, 1.05)
xlim(w[0], w[-1])
xticks(fontsize=18)
yticks(fontsize=18)
text(45, 0.5, r'$\omega_p={0}\pi$'.format(wp/np.pi), fontsize=20)
text(350, 0.5, r'$\omega_s={0}\pi$'.format(ws/np.pi), fontsize=20)
text(10, 0.8, r'$\epsilon={0}$'.format(epsilon), fontsize=20)
text(1000, 0.15, r'$\delta={0}$'.format(delta), fontsize=20)
xlabel(r'Frequency $(\omega)$', fontsize=18)
ylabel(r'$|H(j\omega)|$', fontsize=18)
subplot2grid((1,3), (0,2), rowspan=1, colspan=1)
teta = np.arange(0., 2 * np.pi, 0.01)
_x, _y = w0 * np.cos(teta), w0 * np.sin(teta)
for p in poles:
if p[0] < 0:
plot(p[0], p[1], 'x', color='k', ms=10)
plot(_x, _y, '0.7', ls='--')
plot([-1.1 * w0, 1.1 * w0], [0, 0], '0.8')
plot([0., 0.], [-1.1 * w0, 1.1 * w0], '0.8')
xlim(-1.1 * w0, 1.1 * w0)
ylim(-1.1 * w0, 1.1 * w0)
xticks(fontsize=10)
yticks(fontsize=10)
title('Filter poles'.format(N, w0), fontsize=18)
tight_layout()
savefig("img/butter_filt_example.svg", format="svg", dpi=600)
savefig("img/butter_filt_example.eps", format="eps", dpi=600);
In [14]:
# Chebyshev polynomial
def get_chebyshev_poly(x, N):
"""Returns the value of the Chebyshev polynomial.
"""
return np.array([np.cos(N * np.arccos(_x)) if np.abs(_x) < 1
else np.cosh(N * np.arccosh(_x))
for _x in x])
# Chebyshev filter.
def get_chebyshev_response(w, w0, epsilon, N):
"""Returns the Chebyshev filter magnitude response.
"""
cheby_pol = get_chebyshev_poly(w/w0, N)
return 1. / np.sqrt(1 + np.power(epsilon * cheby_pol, 2))
w = np.arange(0.01, 100, 0.01)
epsilon = 0.707
N_orders = [2, 3, 4, 5]
mag = [get_chebyshev_response(w, 1, epsilon, _n) for _n in N_orders]
figure(figsize=(10, 6))
for i in xrange(len(N_orders)):
semilogx(w, mag[i], label=r'$N={0}$'.format(N_orders[i]))
legend(fontsize=20)
title('Chebyshev filter magnitude response', fontsize=20)
xlabel('Frequency (Hz)', fontsize=20)
ylabel(r'$|H(j\omega)|$', fontsize=20)
xticks(fontsize=18)
yticks(fontsize=18)
xlim(0.01, 100)
ylim(-0.1, 1.1)
tight_layout();
savefig("img/cheby.svg", format="svg", dpi=600)
savefig("img/cheby.eps", format="eps", dpi=600);
In [4]:
# Design of Butterworth filter.
# Chebyshev polynomial
def get_chebyshev_poly(x, N):
"""Returns the value of the Chebyshev polynomial.
"""
return np.array([np.cos(N * np.arccos(_x)) if np.abs(_x) < 1
else np.cosh(N * np.arccosh(_x))
for _x in x])
# Chebyshev filter.
def get_chebyshev_response(w, w0, gamma, N):
"""Returns the Chebyshev filter magnitude response.
"""
cheby_pol = get_chebyshev_poly(w/w0, N)
return 1. / np.sqrt(1 + np.power(gamma * cheby_pol, 2))
def get_chebyshev_param(epsilon, delta, wp, ws):
"""Returns the Chebyshev lowpass filter parameters.
"""
w0 = wp
gamma = np.sqrt(1/np.power(epsilon, 2) - 1)
_n = np.arccosh(np.sqrt((1/np.power(delta, 2) - 1)/np.power(gamma, 2)))
_d = np.arccosh(ws/wp)
N = int(np.ceil(_n/_d))
return int(N), w0, gamma
# Butterworth filter specifications.
epsilon, delta, wp, ws = 0.9, 0.1, 50*np.pi, 100*np.pi
# Filter parameters
N, w0, gamma = get_chebyshev_param(epsilon, delta, wp, ws)
# magnitude response and poles.
w = 2 * np.pi * np.arange(1., 1000., 0.1)
mag = get_chebyshev_response(w, w0, gamma, N)
# poles = get_butterworth_poles(w0, N)
figure(figsize=(10, 5.0))
# subplot2grid((1,3), (0,0), rowspan=1, colspan=2)
semilogx([w[0], ws], np.ones(2), '0.5', ls='--')
semilogx(w, mag, lw=2)
semilogx([w[0], wp], epsilon*np.ones(2), '0.5', ls='--')
semilogx([ws, w[-1]], delta*np.ones(2), '0.5', ls='--')
semilogx(wp*np.ones(2), [-0.01, epsilon], '0.5', ls='--')
semilogx(ws*np.ones(2), [delta, 1.], '0.5', ls='--')
title(r'Chebyshev Filter ($N={0}, \omega_0={1:0.3f}, \gamma={2:0.3f}$)'.format(N, w0, gamma),
fontsize=20)
ylim(-0.01, 1.05)
xlim(w[0], w[-1])
xticks(fontsize=18)
yticks(fontsize=18)
text(45, 0.5, r'$\omega_p={0}\pi$'.format(wp/np.pi), fontsize=20)
text(350, 0.5, r'$\omega_s={0}\pi$'.format(ws/np.pi), fontsize=20)
text(10, 0.8, r'$\epsilon={0}$'.format(epsilon), fontsize=20)
text(1000, 0.15, r'$\delta={0}$'.format(delta), fontsize=20)
xlabel(r'Frequency $(\omega)$', fontsize=18)
ylabel(r'$|H(j\omega)|$', fontsize=18)
tight_layout()
savefig("img/cheby_filt_example.svg", format="svg", dpi=600)
savefig("img/cheby_filt_example.eps", format="eps", dpi=600);
In [5]:
from scipy import signal
In [14]:
# Lowpass filter.
filtp = (0.1, 0.1, 10 * np.pi, 100 * np.pi)
g_max = 20 * np.log10(1 - filtp[0])
g_min = 20 * np.log10(filtp[1])
wp = filtp[2]
ws = filtp[3]
N, wn = signal.buttord(wp, ws, g_max, g_min, analog=True)
print "Lowpass Filter"
print "--------------"
print N, wn
# Highpass filter.
filtp = (0.1, 0.1, 100 * np.pi, 10 * np.pi)
g_max = 20 * np.log10(1 - filtp[0])
g_min = 20 * np.log10(filtp[1])
wp = filtp[2]
ws = filtp[3]
N, wn = signal.cheb1ord(wp, ws, g_max, g_min, analog=True)
print "Highpass Filter"
print "---------------"
print N, wn
# Bandpass filter.
filtp = (0.1, 0.1, [50 * np.pi, 1000 * np.pi], [10 * np.pi, 5000 * np.pi])
g_max = 20 * np.log10(1 - filtp[0])
g_min = 20 * np.log10(filtp[1])
wp = filtp[2]
ws = filtp[3]
N, wn = signal.cheb1ord(wp, ws, g_max, g_min, analog=True)
print "Bandpass Filter"
print "---------------"
print N, wn
# Bandstop filter.
filtp = (0.1, 0.1, [2 * np.pi, 10000 * np.pi], [20 * np.pi, 1000 * np.pi])
N, wn = signal.buttord(wp, ws, g_max, g_min, analog=True)
print "Bandstop Filter"
print "---------------"
print N, wn
In [ ]: