143 lines
3.7 KiB
Python
143 lines
3.7 KiB
Python
# ==== Lab #3: Convolution & Sinusoids ====
|
|
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from scipy.integrate import solve_ivp
|
|
|
|
# Components (same as Lab #2)
|
|
R1, C1 = 10e3, 1e-6
|
|
R2, C2 = 4.7e3, 1e-6
|
|
tau1, tau2 = R1*C1, R2*C2
|
|
R4, R5 = 1e3, 4.7e3
|
|
G = 1 + R5/R4
|
|
|
|
# Time grid
|
|
tau_min = min(tau1, tau2)
|
|
tau_max = max(tau1, tau2)
|
|
dt = 0.0001
|
|
t = np.arange(0, 5 + dt/2, dt)
|
|
N = t.size
|
|
|
|
# Impulse responses (causal)
|
|
h1 = (1/tau1) * np.exp(-t/tau1) # u(t) implicit (0 for t<0)
|
|
h2 = (1/tau2) * np.exp(-t/tau2)
|
|
|
|
# Overall impulse response for cascaded stages (discrete conv scaled by dt)
|
|
h = G * dt * np.convolve(h1, h2, mode='full')
|
|
|
|
# Helper: numeric convolution with step dt (Riemann sums)
|
|
def conv_num(a, b, dt):
|
|
"""Full-length discrete convolution approximating integrals via Riemann sums with step dt."""
|
|
return np.convolve(a, b, mode='full') * dt
|
|
|
|
# ======================
|
|
# 1) Single sinusoids
|
|
# ======================
|
|
A = 1.0
|
|
freqs = [1, 5, 10, 50]
|
|
xmax = 2.0
|
|
|
|
plt.figure(1)
|
|
plt.clf()
|
|
for k, f in enumerate(freqs, start=1):
|
|
w = 2*np.pi*f
|
|
x = A * np.sin(w*t) * (t >= 0)
|
|
y_full = conv_num(x, h, dt)
|
|
y = y_full[:N] # keep same length as x
|
|
|
|
# ODE: tau1*tau2*y'' + (tau1+tau2)*y' + y = G*x(t)
|
|
def rhs(ti, z):
|
|
x_t = A*np.sin(w*ti) * (ti >= 0.0)
|
|
return [
|
|
z[1],
|
|
(G*x_t - (tau1 + tau2)*z[1] - z[0]) / (tau1*tau2)
|
|
]
|
|
z0 = [0.0, 0.0] # zero initial conditions
|
|
sol = solve_ivp(rhs, (t[0], t[-1]), z0, t_eval=t, rtol=1e-9, atol=1e-12)
|
|
y_ode = sol.y[0]
|
|
|
|
ax = plt.subplot(len(freqs), 1, k)
|
|
ax.plot(t, x, label='input')
|
|
ax.plot(t, y, label='output (h*x)')
|
|
ax.plot(t, y_ode, '--', label='output (ODE)', color='r')
|
|
ax.grid(True)
|
|
ax.set_xlabel('t (s)')
|
|
ax.set_ylabel('V')
|
|
ax.set_xlim([0, xmax])
|
|
ax.set_title(f'Single sinusoid: f = {f:g} Hz')
|
|
ax.legend()
|
|
xmax = xmax/2
|
|
plt.tight_layout()
|
|
|
|
plt.show()
|
|
exit()
|
|
|
|
# =========================================
|
|
# 2) Sum of sinusoids
|
|
# =========================================
|
|
A1, f1 = 1.0, 5.0
|
|
A2, f2 = 0.5, 20.0
|
|
x1 = A1 * np.sin(2*np.pi*f1*t) * (t >= 0)
|
|
x2 = A2 * np.sin(2*np.pi*f2*t) * (t >= 0)
|
|
alpha, beta = 1.5, 2.0
|
|
x = alpha*x1 + beta*x2
|
|
Nx = x.size
|
|
|
|
y = conv_num(x, h, dt)[:Nx]
|
|
y1 = conv_num(x1, h, dt)[:Nx]
|
|
y2 = conv_num(x2, h, dt)[:Nx]
|
|
|
|
plt.figure(2)
|
|
plt.clf()
|
|
plt.plot(t, y, label='y = h*(x1+x2)')
|
|
plt.plot(t, alpha*y1 + beta*y2, '--', label='alpha*y1 + beta*y2')
|
|
plt.grid(True)
|
|
plt.xlabel('t (s)')
|
|
plt.ylabel('V')
|
|
plt.xlim([0, 1])
|
|
plt.legend()
|
|
plt.title('Linearity check: y[u1+u2] vs y[u1]+y[u2]')
|
|
|
|
# ==================================
|
|
# 3) Time invariance (delay \Delta)
|
|
# ==================================
|
|
Delta = 0.1
|
|
nD = int(round(Delta/dt))
|
|
x_del = np.concatenate([np.zeros(nD), x[:-nD]]) if nD < Nx else np.zeros_like(x)
|
|
y_del = conv_num(x_del, h, dt)[:x_del.size]
|
|
|
|
plt.figure(3)
|
|
plt.clf()
|
|
plt.plot(t, y, label='y(t)')
|
|
plt.plot(t, np.concatenate([np.zeros(nD), y[:-nD]]), label=r'$y(t-\Delta)$')
|
|
plt.plot(t, y_del, '--', label=r'$h*u(t-\Delta)$')
|
|
plt.grid(True)
|
|
plt.xlabel('t (s)')
|
|
plt.ylabel('V')
|
|
plt.xlim([0, 1])
|
|
plt.legend()
|
|
plt.title('Time invariance check')
|
|
plt.tight_layout()
|
|
|
|
|
|
# =======================================
|
|
# 4) Commutativity of stages (h1*h2)
|
|
# =======================================
|
|
h12 = G * dt * np.convolve(h1, h2, mode='full')
|
|
h21 = G * dt * np.convolve(h2, h1, mode='full')
|
|
rmse = np.sqrt(np.mean((h12 - h21)**2))
|
|
|
|
plt.figure(4)
|
|
plt.clf()
|
|
plt.plot(t, h12[:N], label=r'$h_1 \ast h_2$')
|
|
plt.plot(t, h21[:N], '--', label=r'$h_2 \ast h_1$')
|
|
plt.grid(True)
|
|
plt.xlim([0, 0.1])
|
|
plt.xlabel('t (s)')
|
|
plt.ylabel('V')
|
|
plt.legend()
|
|
plt.title(rf'Commutativity - RMSE of $ e = h_{{12}} - h_{{21}}: \sqrt{{\mathrm{{mean}}(e^2)}} = {rmse:.3e} $')
|
|
|
|
plt.tight_layout()
|
|
plt.show()
|