The Laplace transform is a powerful tool for solving differential equations by converting them to algebraic equations.
The Laplace transform of a function f(t) is defined as:
$\mathcal{L}\{f(t)\} = F(s) = \int_0^{\infty} e^{-st} f(t) dt$
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import sympy as sp
# Symbolic Laplace transforms using SymPy
t, s = sp.symbols('t s', positive=True)
# Common Laplace transforms
print("Common Laplace Transform Pairs:")
print("================================")
# 1. Constant
f1 = 1
F1 = sp.laplace_transform(f1, t, s)
print(f"L{{1}} = {F1[0]}")
# 2. Exponential
a = sp.Symbol('a', real=True)
f2 = sp.exp(a*t)
F2 = sp.laplace_transform(f2, t, s)
print(f"L{{e^(at)}} = {F2[0]}")
# 3. Sine and Cosine
omega = sp.Symbol('omega', positive=True)
f3 = sp.sin(omega*t)
F3 = sp.laplace_transform(f3, t, s)
print(f"L{{sin(ωt)}} = {F3[0]}")
f4 = sp.cos(omega*t)
F4 = sp.laplace_transform(f4, t, s)
print(f"L{{cos(ωt)}} = {F4[0]}")
# 4. Power functions
n = sp.Symbol('n', positive=True, integer=True)
f5 = t**n
F5 = sp.laplace_transform(f5, t, s)
print(f"L{{t^n}} = {F5[0]}")
# Numerical Laplace transform visualization
def numerical_laplace(f, s_val, t_max=10):
"""Compute Laplace transform numerically at s = s_val"""
def integrand(t):
return np.exp(-s_val * t) * f(t)
result, _ = integrate.quad(integrand, 0, t_max)
return result
# Example: Laplace transform of unit step function
t_vals = np.linspace(0, 5, 1000)
s_vals = np.linspace(0.1, 5, 100)
# Unit step function
def unit_step(t):
return np.ones_like(t)
# Compute transform for different s values
F_vals = [numerical_laplace(unit_step, s) for s in s_vals]
plt.figure(figsize=(12, 5))
# Plot original function
plt.subplot(1, 2, 1)
plt.plot(t_vals, unit_step(t_vals), 'b-', linewidth=2)
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Unit Step Function')
plt.grid(True, alpha=0.3)
plt.ylim(-0.1, 1.5)
# Plot Laplace transform
plt.subplot(1, 2, 2)
plt.plot(s_vals, F_vals, 'r-', linewidth=2, label='Numerical')
plt.plot(s_vals, 1/s_vals, 'k--', linewidth=2, label='Exact: 1/s')
plt.xlabel('s')
plt.ylabel('F(s)')
plt.title('Laplace Transform F(s)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(0, 10)
plt.tight_layout()
plt.show()# Demonstrate key properties
print("\nLaplace Transform Properties:")
print("=============================")
# 1. Linearity
f = 3*sp.sin(2*t) + 5*sp.cos(3*t)
F = sp.laplace_transform(f, t, s)
print(f"1. Linearity: L{{3sin(2t) + 5cos(3t)}} = {F[0]}")
# 2. First derivative
f = sp.Function('f')
derivative_property = sp.laplace_transform(sp.diff(f(t), t), t, s)
print(f"\n2. First derivative: L{{f'(t)}} = sF(s) - f(0)")
# 3. Second derivative
print(f" Second derivative: L{{f''(t)}} = s²F(s) - sf(0) - f'(0)")
# 4. Shifting theorem
a_val = 3
f = sp.exp(-a_val*t) * sp.sin(2*t)
F = sp.laplace_transform(f, t, s)
print(f"\n3. Shifting: L{{e^(-3t)sin(2t)}} = {F[0]}")
# 5. Convolution theorem
print(f"\n4. Convolution: L{{f*g}} = F(s)G(s)")
# Visualization of shifting property
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# Original function and its transform
t_vals = np.linspace(0, 5, 1000)
f_original = np.sin(2*np.pi*t_vals)
axes[0, 0].plot(t_vals, f_original, 'b-', linewidth=2)
axes[0, 0].set_xlabel('t')
axes[0, 0].set_ylabel('f(t)')
axes[0, 0].set_title('Original: sin(2πt)')
axes[0, 0].grid(True, alpha=0.3)
# Shifted function
for a in [0.5, 1, 2]:
f_shifted = np.exp(-a*t_vals) * f_original
axes[0, 1].plot(t_vals, f_shifted, linewidth=2, label=f'a={a}')
axes[0, 1].set_xlabel('t')
axes[0, 1].set_ylabel('f(t)')
axes[0, 1].set_title('Shifted: e^(-at)sin(2πt)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# Time scaling property
scales = [0.5, 1, 2]
for idx, scale in enumerate(scales):
f_scaled = np.sin(2*np.pi*scale*t_vals)
axes[1, 0].plot(t_vals, f_scaled, linewidth=2, label=f'k={scale}')
axes[1, 0].set_xlabel('t')
axes[1, 0].set_ylabel('f(t)')
axes[1, 0].set_title('Time Scaling: sin(2πkt)')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# Derivative property illustration
f = t_vals**2 * np.exp(-t_vals)
f_prime = 2*t_vals*np.exp(-t_vals) - t_vals**2*np.exp(-t_vals)
axes[1, 1].plot(t_vals, f, 'b-', linewidth=2, label='f(t) = t²e^(-t)')
axes[1, 1].plot(t_vals, f_prime, 'r--', linewidth=2, label="f'(t)")
axes[1, 1].set_xlabel('t')
axes[1, 1].set_ylabel('f(t)')
axes[1, 1].set_title('Function and its Derivative')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()The Heaviside (unit step) function is fundamental for modeling switched systems and discontinuous inputs.
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
# Define Heaviside function
def heaviside(t, t0=0):
"""Heaviside step function H(t - t0)"""
return np.where(t >= t0, 1, 0)
# Symbolic Heaviside
t, s, a = sp.symbols('t s a', real=True)
H = sp.Heaviside(t)
print("Heaviside Function Properties:")
print("=============================")
print(f"H(t) = 0 for t < 0")
print(f"H(t) = 1 for t > 0")
print(f"H(0) = 1/2 (by convention)")
print(f"\nLaplace transform: L{{H(t)}} = 1/s")
# Shifted Heaviside
H_shifted = sp.Heaviside(t - a)
L_H_shifted = sp.laplace_transform(H_shifted, t, s)
print(f"L{{H(t-a)}} = e^(-as)/s")
# Visualization
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
# Basic Heaviside function
t_vals = np.linspace(-2, 5, 1000)
y_vals = heaviside(t_vals)
axes[0, 0].plot(t_vals, y_vals, 'b-', linewidth=2)
axes[0, 0].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axes[0, 0].axvline(x=0, color='k', linestyle='-', alpha=0.3)
axes[0, 0].set_xlabel('t')
axes[0, 0].set_ylabel('H(t)')
axes[0, 0].set_title('Unit Step Function H(t)')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_ylim(-0.2, 1.2)
# Shifted Heaviside functions
shifts = [1, 2, 3]
for shift in shifts:
y_shifted = heaviside(t_vals, shift)
axes[0, 1].plot(t_vals, y_shifted, linewidth=2, label=f'H(t-{shift})')
axes[0, 1].set_xlabel('t')
axes[0, 1].set_ylabel('H(t-a)')
axes[0, 1].set_title('Shifted Step Functions')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].set_ylim(-0.2, 1.2)
# Rectangular pulse
def rect_pulse(t, a, b):
"""Rectangular pulse from a to b"""
return heaviside(t, a) - heaviside(t, b)
y_rect = rect_pulse(t_vals, 1, 3)
axes[0, 2].plot(t_vals, y_rect, 'g-', linewidth=2)
axes[0, 2].fill_between(t_vals, 0, y_rect, alpha=0.3, color='green')
axes[0, 2].set_xlabel('t')
axes[0, 2].set_ylabel('f(t)')
axes[0, 2].set_title('Rectangular Pulse: H(t-1) - H(t-3)')
axes[0, 2].grid(True, alpha=0.3)
axes[0, 2].set_ylim(-0.2, 1.2)
# Piecewise functions using Heaviside
# f(t) = t for 02
f_piecewise = t_vals * heaviside(t_vals) - (t_vals - 2) * heaviside(t_vals, 2)
axes[1, 0].plot(t_vals, f_piecewise, 'r-', linewidth=2)
axes[1, 0].set_xlabel('t')
axes[1, 0].set_ylabel('f(t)')
axes[1, 0].set_title('Piecewise: f(t) = t·H(t) - (t-2)·H(t-2)')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_ylim(-0.5, 3)
# Staircase function
staircase = np.zeros_like(t_vals)
for i in range(5):
staircase += heaviside(t_vals, i)
axes[1, 1].plot(t_vals, staircase, 'b-', linewidth=2)
axes[1, 1].set_xlabel('t')
axes[1, 1].set_ylabel('f(t)')
axes[1, 1].set_title('Staircase Function')
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].set_ylim(-0.5, 5.5)
# Switching function
# Turn on at t=1, off at t=3, on again at t=4
switching = heaviside(t_vals, 1) - heaviside(t_vals, 3) + heaviside(t_vals, 4)
axes[1, 2].plot(t_vals, switching, 'm-', linewidth=2)
axes[1, 2].fill_between(t_vals, 0, switching, alpha=0.3, color='magenta')
axes[1, 2].set_xlabel('t')
axes[1, 2].set_ylabel('f(t)')
axes[1, 2].set_title('Switching Function')
axes[1, 2].grid(True, alpha=0.3)
axes[1, 2].set_ylim(-0.2, 1.2)
plt.tight_layout()
plt.show()
# Applications with differential equations
print("\nApplication: Solving ODEs with discontinuous forcing")
print("=================================================")
print("Consider: y' + 2y = H(t-1), y(0) = 0")
print("This models a system with input turned on at t=1")
# Numerical solution
from scipy.integrate import odeint
def system_ode(y, t):
"""y' + 2y = H(t-1)"""
forcing = heaviside(t, 1)
return -2*y + forcing
t_span = np.linspace(0, 5, 1000)
y0 = [0]
solution = odeint(system_ode, y0, t_span)
plt.figure(figsize=(10, 6))
plt.plot(t_span, solution[:, 0], 'b-', linewidth=2, label='System response')
plt.plot(t_span, heaviside(t_span, 1), 'r--', linewidth=2, label='Input H(t-1)')
plt.axvline(x=1, color='gray', linestyle=':', alpha=0.7, label='Switch on')
plt.xlabel('Time t')
plt.ylabel('y(t)')
plt.title('Response to Step Input')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show() # Demonstrate the second shifting theorem
print("\nSecond Shifting Theorem:")
print("======================")
print("If L{f(t)} = F(s), then:")
print("L{H(t-a)·f(t-a)} = e^(-as)·F(s)")
# Example: Delayed exponential decay
t_vals = np.linspace(0, 10, 1000)
a = 2 # delay
# Original function: e^(-t)
f_original = np.exp(-t_vals)
# Delayed function: H(t-2)·e^(-(t-2))
f_delayed = heaviside(t_vals, a) * np.exp(-(t_vals - a))
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t_vals, f_original, 'b-', linewidth=2, label='e^(-t)')
plt.plot(t_vals, f_delayed, 'r--', linewidth=2, label='H(t-2)·e^(-(t-2))')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Original and Delayed Functions')
plt.legend()
plt.grid(True, alpha=0.3)
# Multiple delays
plt.subplot(1, 2, 2)
delays = [0, 1, 2, 3]
for delay in delays:
f_multi = heaviside(t_vals, delay) * np.exp(-(t_vals - delay))
plt.plot(t_vals, f_multi, linewidth=2, label=f'Delay = {delay}')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Multiple Delayed Exponentials')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()Real-world systems often involve discontinuous inputs like switches, pulses, and periodic signals.
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import sympy as sp
# Example 1: Square wave
def square_wave(t, period=2*np.pi, duty=0.5):
"""Generate square wave with given period and duty cycle"""
phase = (t % period) / period
return np.where(phase < duty, 1, -1)
# Example 2: Sawtooth wave
def sawtooth_wave(t, period=2*np.pi):
"""Generate sawtooth wave"""
return 2 * ((t % period) / period) - 1
# Example 3: Triangular wave
def triangular_wave(t, period=2*np.pi):
"""Generate triangular wave"""
phase = (t % period) / period
return np.where(phase < 0.5, 4*phase - 1, 3 - 4*phase)
# Visualization
t = np.linspace(0, 6*np.pi, 1000)
fig, axes = plt.subplots(3, 2, figsize=(12, 10))
# Square wave
y_square = square_wave(t)
axes[0, 0].plot(t, y_square, 'b-', linewidth=2)
axes[0, 0].set_xlabel('t')
axes[0, 0].set_ylabel('f(t)')
axes[0, 0].set_title('Square Wave')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_ylim(-1.5, 1.5)
# Fourier series approximation of square wave
fourier_approx = np.zeros_like(t)
for n in range(1, 20, 2):
fourier_approx += (4/np.pi) * np.sin(n*t) / n
axes[0, 1].plot(t, y_square, 'b-', linewidth=2, alpha=0.5, label='Square wave')
axes[0, 1].plot(t, fourier_approx, 'r--', linewidth=2, label='Fourier approx')
axes[0, 1].set_xlabel('t')
axes[0, 1].set_ylabel('f(t)')
axes[0, 1].set_title('Fourier Series Approximation')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# Sawtooth wave
y_sawtooth = sawtooth_wave(t)
axes[1, 0].plot(t, y_sawtooth, 'g-', linewidth=2)
axes[1, 0].set_xlabel('t')
axes[1, 0].set_ylabel('f(t)')
axes[1, 0].set_title('Sawtooth Wave')
axes[1, 0].grid(True, alpha=0.3)
# Triangular wave
y_triangular = triangular_wave(t)
axes[1, 1].plot(t, y_triangular, 'm-', linewidth=2)
axes[1, 1].set_xlabel('t')
axes[1, 1].set_ylabel('f(t)')
axes[1, 1].set_title('Triangular Wave')
axes[1, 1].grid(True, alpha=0.3)
# Pulse train
def pulse_train(t, width=0.1, period=1):
"""Generate pulse train"""
result = np.zeros_like(t)
for n in range(int(t[-1] / period) + 1):
result += (heaviside(t, n*period) - heaviside(t, n*period + width))
return result
def heaviside(t, t0):
return np.where(t >= t0, 1, 0)
y_pulse = pulse_train(t/np.pi, width=0.2, period=2)
axes[2, 0].plot(t, y_pulse, 'r-', linewidth=2)
axes[2, 0].set_xlabel('t')
axes[2, 0].set_ylabel('f(t)')
axes[2, 0].set_title('Pulse Train')
axes[2, 0].grid(True, alpha=0.3)
axes[2, 0].set_ylim(-0.2, 1.2)
# Rectified sine wave
y_rectified = np.abs(np.sin(t))
axes[2, 1].plot(t, y_rectified, 'c-', linewidth=2)
axes[2, 1].set_xlabel('t')
axes[2, 1].set_ylabel('f(t)')
axes[2, 1].set_title('Rectified Sine Wave')
axes[2, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Laplace transform of periodic functions
print("Laplace Transform of Periodic Functions")
print("======================================")
print("For a periodic function f(t) with period T:")
print("L{f(t)} = (1/(1-e^(-sT))) * ∫[0 to T] e^(-st)f(t)dt")
# Example: Laplace transform of square wave
t_sym, s_sym = sp.symbols('t s', positive=True)
T = sp.Symbol('T', positive=True)
# Square wave from 0 to T/2 is 1, from T/2 to T is 0
print("\nSquare wave (0 to T/2: 1, T/2 to T: 0):")
print("L{f(t)} = (1-e^(-sT/2))/(s(1-e^(-sT)))")
# Numerical example
def laplace_square_wave(s_val, T=2*np.pi):
"""Laplace transform of square wave"""
if s_val == 0:
return np.inf
return (1 - np.exp(-s_val*T/2)) / (s_val * (1 - np.exp(-s_val*T)))
s_vals = np.linspace(0.1, 5, 100)
F_vals = [laplace_square_wave(s) for s in s_vals]
plt.figure(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.plot(s_vals, np.real(F_vals), 'b-', linewidth=2)
plt.xlabel('s')
plt.ylabel('Re[F(s)]')
plt.title('Laplace Transform of Square Wave (Real Part)')
plt.grid(True, alpha=0.3)
# Phase plot
plt.subplot(1, 2, 2)
poles = [2*np.pi*n*1j for n in range(-5, 6) if n != 0]
plt.scatter([0], [0], color='red', s=100, marker='x', label='Pole at s=0')
for pole in poles:
plt.scatter([0], [pole.imag], color='red', s=50, marker='x')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
plt.xlabel('Real(s)')
plt.ylabel('Imag(s)')
plt.title('Poles of F(s)')
plt.grid(True, alpha=0.3)
plt.xlim(-2, 2)
plt.ylim(-20, 20)
plt.legend()
plt.tight_layout()
plt.show()# Common window functions
t = np.linspace(-1, 5, 1000)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# Rectangular window
rect = heaviside(t, 0) - heaviside(t, 2)
axes[0, 0].plot(t, rect, 'b-', linewidth=2)
axes[0, 0].fill_between(t, 0, rect, alpha=0.3)
axes[0, 0].set_title('Rectangular Window')
axes[0, 0].set_xlabel('t')
axes[0, 0].set_ylabel('f(t)')
axes[0, 0].grid(True, alpha=0.3)
# Exponential window
exp_window = heaviside(t, 0) * np.exp(-t) * (1 - heaviside(t, 3))
axes[0, 1].plot(t, exp_window, 'g-', linewidth=2)
axes[0, 1].fill_between(t, 0, exp_window, alpha=0.3, color='green')
axes[0, 1].set_title('Exponential Window')
axes[0, 1].set_xlabel('t')
axes[0, 1].set_ylabel('f(t)')
axes[0, 1].grid(True, alpha=0.3)
# Gaussian-like pulse
gaussian_pulse = heaviside(t, 0) * np.exp(-2*(t-1)**2) * (1 - heaviside(t, 3))
axes[1, 0].plot(t, gaussian_pulse, 'r-', linewidth=2)
axes[1, 0].fill_between(t, 0, gaussian_pulse, alpha=0.3, color='red')
axes[1, 0].set_title('Gaussian-like Pulse')
axes[1, 0].set_xlabel('t')
axes[1, 0].set_ylabel('f(t)')
axes[1, 0].grid(True, alpha=0.3)
# Ramp window
ramp = t * heaviside(t, 0) - t * heaviside(t, 1) + heaviside(t, 1) - heaviside(t, 2)
axes[1, 1].plot(t, ramp, 'm-', linewidth=2)
axes[1, 1].fill_between(t, 0, ramp, alpha=0.3, color='magenta')
axes[1, 1].set_title('Ramp Window')
axes[1, 1].set_xlabel('t')
axes[1, 1].set_ylabel('f(t)')
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Application: RC circuit with pulse input
print("\nApplication: RC Circuit with Pulse Input")
print("======================================")
from scipy.integrate import odeint
# RC circuit: RC*dv/dt + v = V_in(t)
R = 1.0 # Resistance
C = 0.5 # Capacitance
RC = R * C
def rc_circuit(v, t):
# Pulse input from t=1 to t=3
V_in = heaviside(t, 1) - heaviside(t, 3)
return (V_in - v) / RC
t_span = np.linspace(0, 6, 1000)
v0 = [0] # Initial voltage
solution = odeint(rc_circuit, v0, t_span)
V_in = heaviside(t_span, 1) - heaviside(t_span, 3)
plt.figure(figsize=(10, 6))
plt.plot(t_span, V_in, 'r--', linewidth=2, label='Input voltage')
plt.plot(t_span, solution[:, 0], 'b-', linewidth=2, label='Capacitor voltage')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.title('RC Circuit Response to Pulse Input')
plt.legend()
plt.grid(True, alpha=0.3)
plt.axvline(x=1, color='gray', linestyle=':', alpha=0.5)
plt.axvline(x=3, color='gray', linestyle=':', alpha=0.5)
plt.show()Converting from the s-domain back to the time domain using partial fractions and standard transforms.
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from sympy import symbols, apart, inverse_laplace_transform, exp, sin, cos
# Define symbols
s, t = symbols('s t', real=True)
print("Inverse Laplace Transform Examples")
print("=================================\n")
# Example 1: Simple poles
F1 = 1 / ((s + 1) * (s + 2))
print("Example 1: F(s) = 1/((s+1)(s+2))")
print(f"Partial fractions: {apart(F1, s)}")
f1 = inverse_laplace_transform(F1, s, t)
print(f"f(t) = {f1}\n")
# Example 2: Repeated poles
F2 = 1 / (s * (s + 1)**2)
print("Example 2: F(s) = 1/(s(s+1)²)")
print(f"Partial fractions: {apart(F2, s)}")
f2 = inverse_laplace_transform(F2, s, t)
print(f"f(t) = {f2}\n")
# Example 3: Complex poles
F3 = (s + 1) / (s**2 + 2*s + 5)
print("Example 3: F(s) = (s+1)/(s²+2s+5)")
f3 = inverse_laplace_transform(F3, s, t)
print(f"f(t) = {f3}\n")
# Numerical evaluation and plotting
t_vals = np.linspace(0, 5, 1000)
# Convert symbolic expressions to numerical functions
f1_num = sp.lambdify(t, f1 * sp.Heaviside(t), 'numpy')
f2_num = sp.lambdify(t, f2 * sp.Heaviside(t), 'numpy')
f3_num = sp.lambdify(t, f3 * sp.Heaviside(t), 'numpy')
plt.figure(figsize=(12, 8))
# Plot Example 1
plt.subplot(2, 2, 1)
y1 = f1_num(t_vals)
plt.plot(t_vals, y1, 'b-', linewidth=2)
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Example 1: Simple Poles')
plt.grid(True, alpha=0.3)
# Plot Example 2
plt.subplot(2, 2, 2)
y2 = f2_num(t_vals)
plt.plot(t_vals, y2, 'g-', linewidth=2)
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Example 2: Repeated Poles')
plt.grid(True, alpha=0.3)
# Plot Example 3
plt.subplot(2, 2, 3)
y3 = f3_num(t_vals)
plt.plot(t_vals, y3, 'r-', linewidth=2)
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Example 3: Complex Poles (Damped Oscillation)')
plt.grid(True, alpha=0.3)
# Example 4: Step response of second-order system
F4 = 1 / (s * (s**2 + 2*s + 2))
f4 = inverse_laplace_transform(F4, s, t)
f4_num = sp.lambdify(t, f4 * sp.Heaviside(t), 'numpy')
plt.subplot(2, 2, 4)
y4 = f4_num(t_vals)
plt.plot(t_vals, y4, 'm-', linewidth=2)
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Example 4: Second-Order Step Response')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()# Demonstrate various inverse Laplace techniques
print("\nInverse Laplace Transform Techniques")
print("===================================")
# 1. Completing the square
print("\n1. Completing the Square:")
F = 1 / (s**2 + 4*s + 13)
print(f" F(s) = {F}")
print(f" Rewrite as: 1/((s+2)² + 9)")
f = inverse_laplace_transform(F, s, t)
print(f" f(t) = {f}")
# 2. First Shifting Theorem
print("\n2. First Shifting Theorem:")
print(" If L{f(t)} = F(s), then L{e^(at)f(t)} = F(s-a)")
F_shifted = 1 / ((s - 2)**2 + 4)
f_shifted = inverse_laplace_transform(F_shifted, s, t)
print(f" F(s) = 1/((s-2)² + 4)")
print(f" f(t) = {f_shifted}")
# 3. Differentiation in s-domain
print("\n3. Differentiation in s-domain:")
print(" If L{f(t)} = F(s), then L{t·f(t)} = -dF/ds")
F_base = 1 / (s + 3)
F_diff = sp.diff(F_base, s)
print(f" F(s) = 1/(s+3), -dF/ds = {-F_diff}")
f_mult = inverse_laplace_transform(-F_diff, s, t)
print(f" f(t) = {f_mult}")
# Visualization of transform pairs
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
# Common transform pairs
transform_pairs = [
(1/s, "1", "Unit step"),
(1/s**2, "t", "Ramp"),
(1/(s**2 + 4), "(1/2)sin(2t)", "Sine"),
(s/(s**2 + 4), "cos(2t)", "Cosine"),
(1/(s + 2), "e^(-2t)", "Exponential"),
(1/(s + 1)**2, "t·e^(-t)", "t×exponential")
]
t_vals = np.linspace(0, 5, 500)
for idx, (F_expr, f_str, title) in enumerate(transform_pairs):
row = idx // 3
col = idx % 3
# Get inverse transform
f_expr = inverse_laplace_transform(F_expr, s, t)
f_num = sp.lambdify(t, f_expr * sp.Heaviside(t), 'numpy')
y_vals = f_num(t_vals)
axes[row, col].plot(t_vals, y_vals, linewidth=2)
axes[row, col].set_xlabel('t')
axes[row, col].set_ylabel('f(t)')
axes[row, col].set_title(f'{title}\nF(s) = {F_expr}')
axes[row, col].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Table of common transforms
print("\nCommon Laplace Transform Pairs:")
print("==============================")
print("f(t) | F(s)")
print("--------------------- | ---------------------")
print("1 | 1/s")
print("t^n | n!/s^(n+1)")
print("e^(at) | 1/(s-a)")
print("sin(ωt) | ω/(s²+ω²)")
print("cos(ωt) | s/(s²+ω²)")
print("e^(at)sin(ωt) | ω/((s-a)²+ω²)")
print("e^(at)cos(ωt) | (s-a)/((s-a)²+ω²)")
print("t^n·e^(at) | n!/(s-a)^(n+1)")
print("H(t-a) | e^(-as)/s")
print("δ(t-a) | e^(-as)")Enter coefficients for F(s) = (as + b) / (s² + cs + d):
Applying Laplace transforms to solve real-world engineering problems in mechanical and electrical systems.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import sympy as sp
# Mass-Spring-Damper System
print("Mass-Spring-Damper System")
print("========================")
print("Equation: m·x'' + c·x' + k·x = F(t)")
print("where m = mass, c = damping, k = spring constant\n")
# System parameters
m = 1.0 # mass (kg)
c = 0.5 # damping coefficient (N·s/m)
k = 4.0 # spring constant (N/m)
# Natural frequency and damping ratio
omega_n = np.sqrt(k/m)
zeta = c / (2 * np.sqrt(k * m))
print(f"Natural frequency ωn = {omega_n:.2f} rad/s")
print(f"Damping ratio ζ = {zeta:.2f}")
if zeta < 1:
print("System is underdamped")
elif zeta == 1:
print("System is critically damped")
else:
print("System is overdamped")
# Symbolic Laplace analysis
s, t = sp.symbols('s t', real=True)
# Transfer function H(s) = X(s)/F(s) = 1/(ms² + cs + k)
H_s = 1 / (m*s**2 + c*s + k)
print(f"\nTransfer function: H(s) = {H_s}")
# Step response (F(t) = 1 for t > 0)
X_s = H_s / s # Step input has Laplace transform 1/s
x_step = sp.inverse_laplace_transform(X_s, s, t)
print(f"Step response: x(t) = {x_step}")
# Numerical simulation
def mass_spring_damper(state, t, F_func):
x, v = state
F = F_func(t)
x_dot = v
v_dot = (F - c*v - k*x) / m
return [x_dot, v_dot]
t_vals = np.linspace(0, 10, 1000)
# Different forcing functions
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
# 1. Step input
F_step = lambda t: 1.0
y0 = [0, 0] # Initial position and velocity
sol_step = odeint(mass_spring_damper, y0, t_vals, args=(F_step,))
axes[0, 0].plot(t_vals, sol_step[:, 0], 'b-', linewidth=2)
axes[0, 0].axhline(y=1/k, color='r', linestyle='--', label=f'Steady state = {1/k:.3f}')
axes[0, 0].set_xlabel('Time (s)')
axes[0, 0].set_ylabel('Position (m)')
axes[0, 0].set_title('Step Response')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# 2. Impulse input
F_impulse = lambda t: 10.0 if t < 0.1 else 0
sol_impulse = odeint(mass_spring_damper, y0, t_vals, args=(F_impulse,))
axes[0, 1].plot(t_vals, sol_impulse[:, 0], 'g-', linewidth=2)
axes[0, 1].set_xlabel('Time (s)')
axes[0, 1].set_ylabel('Position (m)')
axes[0, 1].set_title('Impulse Response')
axes[0, 1].grid(True, alpha=0.3)
# 3. Sinusoidal input
omega_force = 2.0 # Forcing frequency
F_sin = lambda t: np.sin(omega_force * t)
sol_sin = odeint(mass_spring_damper, y0, t_vals, args=(F_sin,))
axes[0, 2].plot(t_vals, sol_sin[:, 0], 'r-', linewidth=2, label='Response')
axes[0, 2].plot(t_vals, F_sin(t_vals)/k, 'k--', alpha=0.5, label='Input/k')
axes[0, 2].set_xlabel('Time (s)')
axes[0, 2].set_ylabel('Position (m)')
axes[0, 2].set_title('Sinusoidal Response')
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)
# Frequency response (Bode plot)
omega_range = np.logspace(-1, 1, 200)
H_mag = np.zeros_like(omega_range)
H_phase = np.zeros_like(omega_range)
for i, omega in enumerate(omega_range):
H_jw = 1 / (k - m*omega**2 + 1j*c*omega)
H_mag[i] = np.abs(H_jw)
H_phase[i] = np.angle(H_jw)
axes[1, 0].loglog(omega_range, H_mag, 'b-', linewidth=2)
axes[1, 0].axvline(x=omega_n, color='r', linestyle='--', label='ωn')
axes[1, 0].set_xlabel('Frequency (rad/s)')
axes[1, 0].set_ylabel('|H(jω)|')
axes[1, 0].set_title('Magnitude Response')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 1].semilogx(omega_range, np.degrees(H_phase), 'g-', linewidth=2)
axes[1, 1].axvline(x=omega_n, color='r', linestyle='--', label='ωn')
axes[1, 1].set_xlabel('Frequency (rad/s)')
axes[1, 1].set_ylabel('Phase (degrees)')
axes[1, 1].set_title('Phase Response')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
# Different damping ratios
zeta_values = [0.1, 0.5, 1.0, 2.0]
for zeta_val in zeta_values:
c_val = 2 * zeta_val * np.sqrt(k * m)
sol = odeint(mass_spring_damper, y0, t_vals, args=(F_step,))
# Recompute with new damping
def temp_system(state, t):
x, v = state
return [v, (1 - c_val*v - k*x) / m]
sol = odeint(temp_system, y0, t_vals)
axes[1, 2].plot(t_vals, sol[:, 0], linewidth=2, label=f'ζ={zeta_val}')
axes[1, 2].set_xlabel('Time (s)')
axes[1, 2].set_ylabel('Position (m)')
axes[1, 2].set_title('Effect of Damping Ratio')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()# RLC Circuit Analysis
print("\nRLC Circuit Analysis")
print("===================")
print("Series RLC: L·i'' + R·i' + (1/C)·i = v'(t)")
print("Parallel RLC: C·v'' + (1/R)·v' + (1/L)·v = i(t)\n")
# Series RLC parameters
R = 10.0 # Resistance (Ω)
L = 0.1 # Inductance (H)
C = 1e-3 # Capacitance (F)
# Characteristic parameters
omega_0 = 1 / np.sqrt(L * C) # Resonant frequency
alpha = R / (2 * L) # Damping factor
Q = omega_0 / (2 * alpha) # Quality factor
print(f"Resonant frequency f0 = {omega_0/(2*np.pi):.1f} Hz")
print(f"Damping factor α = {alpha:.1f} s^-1")
print(f"Quality factor Q = {Q:.2f}")
# Transfer function for voltage across capacitor
# H(s) = V_C(s)/V_in(s) = 1/(LCs² + RCs + 1)
H_rlc = 1 / (L*C*s**2 + R*C*s + 1)
print(f"\nTransfer function: H(s) = {H_rlc}")
# Different circuit configurations
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
# 1. Step response of RLC circuit
def rlc_circuit(state, t, V_func):
i, v_c = state # Current and capacitor voltage
v_in = V_func(t)
di_dt = (v_in - R*i - v_c) / L
dv_c_dt = i / C
return [di_dt, dv_c_dt]
t_vals = np.linspace(0, 0.01, 1000) # 10 ms
V_step = lambda t: 1.0 # 1V step
y0 = [0, 0] # Initial current and capacitor voltage
sol_rlc = odeint(rlc_circuit, y0, t_vals, args=(V_step,))
axes[0, 0].plot(t_vals*1000, sol_rlc[:, 1], 'b-', linewidth=2)
axes[0, 0].set_xlabel('Time (ms)')
axes[0, 0].set_ylabel('V_C (V)')
axes[0, 0].set_title('RLC Step Response')
axes[0, 0].grid(True, alpha=0.3)
# 2. Current through inductor
axes[0, 1].plot(t_vals*1000, sol_rlc[:, 0]*1000, 'g-', linewidth=2)
axes[0, 1].set_xlabel('Time (ms)')
axes[0, 1].set_ylabel('Current (mA)')
axes[0, 1].set_title('Inductor Current')
axes[0, 1].grid(True, alpha=0.3)
# 3. AC steady-state response
freq_range = np.logspace(1, 5, 200) # 10 Hz to 100 kHz
H_mag = np.zeros_like(freq_range)
H_phase = np.zeros_like(freq_range)
for i, f in enumerate(freq_range):
omega = 2 * np.pi * f
H_jw = 1 / (1 - L*C*omega**2 + 1j*R*C*omega)
H_mag[i] = 20 * np.log10(np.abs(H_jw)) # Convert to dB
H_phase[i] = np.degrees(np.angle(H_jw))
axes[0, 2].semilogx(freq_range, H_mag, 'r-', linewidth=2)
axes[0, 2].axvline(x=omega_0/(2*np.pi), color='k', linestyle='--', label='f0')
axes[0, 2].set_xlabel('Frequency (Hz)')
axes[0, 2].set_ylabel('Magnitude (dB)')
axes[0, 2].set_title('Frequency Response')
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)
# 4. Op-amp active filter
print("\nActive Filter (Sallen-Key Low-pass)")
print("===================================")
# Sallen-Key parameters
R1 = R2 = 10e3 # 10 kΩ
C1 = C2 = 10e-9 # 10 nF
fc = 1 / (2 * np.pi * np.sqrt(R1*R2*C1*C2))
print(f"Cutoff frequency: {fc:.1f} Hz")
# Butterworth response
H_butter = 1 / (1 + np.sqrt(2)*s/(2*np.pi*fc) + (s/(2*np.pi*fc))**2)
# Plot filter responses
filter_types = [
("Butterworth", 1, np.sqrt(2)),
("Chebyshev", 1, 1.0),
("Bessel", 1, 2.0)
]
for idx, (name, gain, damping) in enumerate(filter_types):
H_filter = gain / (1 + damping*s/(2*np.pi*fc) + (s/(2*np.pi*fc))**2)
# Evaluate frequency response
H_vals = []
for f in freq_range:
s_val = 2j * np.pi * f
H_val = gain / (1 + damping*s_val/(2*np.pi*fc) + (s_val/(2*np.pi*fc))**2)
H_vals.append(20 * np.log10(np.abs(H_val)))
axes[1, 0].semilogx(freq_range, H_vals, linewidth=2, label=name)
axes[1, 0].axvline(x=fc, color='k', linestyle='--', alpha=0.5)
axes[1, 0].axhline(y=-3, color='k', linestyle=':', alpha=0.5)
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('Magnitude (dB)')
axes[1, 0].set_title('Active Filter Comparison')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_ylim(-40, 5)
# 5. Transmission line (telegraphers equation)
print("\nTransmission Line Model")
print("======================")
# Distributed parameters (per unit length)
R_line = 0.1 # Ω/m
L_line = 1e-6 # H/m
C_line = 100e-12 # F/m
G_line = 1e-6 # S/m
# Characteristic impedance and propagation constant
Z0 = np.sqrt((R_line + 1j*2*np.pi*1e6*L_line) / (G_line + 1j*2*np.pi*1e6*C_line))
gamma = np.sqrt((R_line + 1j*2*np.pi*1e6*L_line) * (G_line + 1j*2*np.pi*1e6*C_line))
print(f"Characteristic impedance Z0 = {np.abs(Z0):.1f} Ω")
print(f"Attenuation constant α = {np.real(gamma):.3f} Np/m")
print(f"Phase constant β = {np.imag(gamma):.3f} rad/m")
# Plot transmission line response
length = 100 # meters
position = np.linspace(0, length, 200)
# Voltage along line for matched load
V_in = 1.0
V_line = V_in * np.exp(-gamma * position)
axes[1, 1].plot(position, np.abs(V_line), 'b-', linewidth=2, label='Magnitude')
axes[1, 1].set_xlabel('Position (m)')
axes[1, 1].set_ylabel('|V(x)| / V_in')
axes[1, 1].set_title('Voltage Along Transmission Line')
axes[1, 1].grid(True, alpha=0.3)
# 6. Switch-mode power supply model
t_switch = np.linspace(0, 0.1, 1000)
duty_cycle = 0.3
V_in_switch = 12 # Input voltage
# PWM signal
pwm = np.where((t_switch % 0.001) < duty_cycle * 0.001, V_in_switch, 0)
# Low-pass filter output (simplified)
V_out = duty_cycle * V_in_switch * np.ones_like(t_switch)
ripple = 0.5 * np.sin(2*np.pi*1000*t_switch) * np.exp(-100*t_switch)
axes[1, 2].plot(t_switch*1000, pwm, 'r-', alpha=0.5, label='PWM input')
axes[1, 2].plot(t_switch*1000, V_out + ripple, 'b-', linewidth=2, label='Filtered output')
axes[1, 2].set_xlabel('Time (ms)')
axes[1, 2].set_ylabel('Voltage (V)')
axes[1, 2].set_title('Switch-Mode Power Supply')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)
axes[1, 2].set_xlim(0, 10)
plt.tight_layout()
plt.show()Exploring how Laplace transforms help analyze and design audio systems and musical instruments.
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import sympy as sp
# Audio equalizer design using Laplace transforms
print("Audio Equalizer Design")
print("=====================\n")
# Standard audio frequency bands
freq_bands = {
'Sub-bass': (20, 60),
'Bass': (60, 250),
'Low-mid': (250, 500),
'Mid': (500, 2000),
'High-mid': (2000, 4000),
'Presence': (4000, 6000),
'Brilliance': (6000, 20000)
}
# Design filters for each band
fs = 44100 # Sampling frequency (Hz)
fig, axes = plt.subplots(3, 2, figsize=(12, 10))
# 1. Parametric EQ filter
print("Parametric EQ Design:")
fc = 1000 # Center frequency
Q = 2 # Quality factor
gain_db = 6 # Gain in dB
# Convert to angular frequency
omega_c = 2 * np.pi * fc
# Transfer function for peaking filter
s = sp.Symbol('s')
A = 10**(gain_db/40)
H_peak = (s**2 + (A/Q)*omega_c*s + omega_c**2) / (s**2 + (omega_c/Q)*s + omega_c**2)
print(f"Center frequency: {fc} Hz")
print(f"Q factor: {Q}")
print(f"Gain: {gain_db} dB")
# Frequency response
w = np.logspace(1, 4.3, 1000) * 2 * np.pi
H_vals = []
for omega in w:
H_val = complex(H_peak.subs(s, 1j*omega))
H_vals.append(H_val)
H_mag = 20 * np.log10(np.abs(H_vals))
H_phase = np.degrees(np.angle(H_vals))
axes[0, 0].semilogx(w/(2*np.pi), H_mag, 'b-', linewidth=2)
axes[0, 0].axvline(x=fc, color='r', linestyle='--', alpha=0.5, label='fc')
axes[0, 0].axhline(y=gain_db, color='g', linestyle=':', alpha=0.5, label=f'{gain_db} dB')
axes[0, 0].set_xlabel('Frequency (Hz)')
axes[0, 0].set_ylabel('Magnitude (dB)')
axes[0, 0].set_title('Parametric EQ Response')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_ylim(-3, gain_db+3)
# 2. Crossover network design
print("\nCrossover Network Design:")
fc_low = 800 # Crossover frequency for woofer/midrange
fc_high = 5000 # Crossover frequency for midrange/tweeter
# Linkwitz-Riley 4th order filters
omega_low = 2 * np.pi * fc_low
omega_high = 2 * np.pi * fc_high
# Low-pass for woofer
H_lp = omega_low**4 / ((s**2 + np.sqrt(2)*omega_low*s + omega_low**2)**2)
# Band-pass for midrange
H_bp_num = (omega_high**2 - omega_low**2)**2 * s**4
H_bp_den = ((s**2 + np.sqrt(2)*omega_low*s + omega_low**2) *
(s**2 + np.sqrt(2)*omega_high*s + omega_high**2))**2
H_bp = H_bp_num / H_bp_den
# High-pass for tweeter
H_hp = s**4 / ((s**2 + np.sqrt(2)*omega_high*s + omega_high**2)**2)
# Plot crossover responses
for H, label, color in [(H_lp, 'Woofer', 'r'),
(H_bp, 'Midrange', 'g'),
(H_hp, 'Tweeter', 'b')]:
H_cross = []
for omega in w:
try:
H_val = complex(H.subs(s, 1j*omega))
H_cross.append(H_val)
except:
H_cross.append(0)
H_mag = 20 * np.log10(np.abs(H_cross) + 1e-10)
axes[0, 1].semilogx(w/(2*np.pi), H_mag, color=color, linewidth=2, label=label)
axes[0, 1].axvline(x=fc_low, color='k', linestyle='--', alpha=0.3)
axes[0, 1].axvline(x=fc_high, color='k', linestyle='--', alpha=0.3)
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('Magnitude (dB)')
axes[0, 1].set_title('3-Way Crossover Network')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].set_ylim(-60, 5)
axes[0, 1].set_xlim(10, 20000)
# 3. Guitar amplifier tone stack
print("\nGuitar Amp Tone Stack:")
# Fender-style tone stack simulation
R1, R2, R3 = 250e3, 1e6, 25e3
C1, C2, C3 = 250e-12, 20e-9, 20e-9
# Simplified transfer function (treble at 50%)
H_tone = (R2*C2*s + 1) / ((R1+R2)*C1*s + R2*C2*s + 1)
# Different tone settings
tone_settings = [
("Flat", 1.0),
("Scooped", 0.5),
("Bright", 1.5)
]
for setting, factor in tone_settings:
H_mod = H_tone * factor
H_tone_vals = []
for omega in w:
H_val = complex(H_mod.subs(s, 1j*omega))
H_tone_vals.append(H_val)
H_mag = 20 * np.log10(np.abs(H_tone_vals))
axes[1, 0].semilogx(w/(2*np.pi), H_mag, linewidth=2, label=setting)
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('Magnitude (dB)')
axes[1, 0].set_title('Guitar Amp Tone Stack')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_xlim(20, 20000)
# 4. Room acoustics modeling
print("\nRoom Acoustics Transfer Function:")
# Simple room model with reflections
room_length = 10 # meters
c_sound = 343 # m/s
# Direct path and first reflection
t_direct = room_length / c_sound
t_reflect = 3 * room_length / c_sound
# Transfer function with direct sound and reflection
# H(s) = 1 + 0.7*e^(-s*t_reflect)
H_room = 1 + 0.7 * sp.exp(-s * t_reflect)
# Magnitude response showing comb filtering
H_room_vals = []
for omega in w:
H_val = 1 + 0.7 * np.exp(-1j * omega * t_reflect)
H_room_vals.append(H_val)
H_mag = 20 * np.log10(np.abs(H_room_vals))
axes[1, 1].semilogx(w/(2*np.pi), H_mag, 'purple', linewidth=2)
axes[1, 1].set_xlabel('Frequency (Hz)')
axes[1, 1].set_ylabel('Magnitude (dB)')
axes[1, 1].set_title('Room Transfer Function (Comb Filtering)')
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].set_ylim(-10, 10)
# 5. Dynamic range compression
print("\nDynamic Range Compressor:")
# Attack and release time constants
tau_attack = 0.001 # 1 ms
tau_release = 0.1 # 100 ms
# Time response simulation
t_comp = np.linspace(0, 0.5, 5000)
# Input signal: sudden loud transient
input_signal = np.zeros_like(t_comp)
input_signal[1000:1500] = 2.0 # Loud section
input_signal[2500:3000] = 1.5 # Another loud section
# Envelope follower
envelope = np.zeros_like(t_comp)
for i in range(1, len(t_comp)):
dt = t_comp[i] - t_comp[i-1]
if np.abs(input_signal[i]) > envelope[i-1]:
# Attack
envelope[i] = envelope[i-1] + (np.abs(input_signal[i]) - envelope[i-1]) * (1 - np.exp(-dt/tau_attack))
else:
# Release
envelope[i] = envelope[i-1] * np.exp(-dt/tau_release)
# Compression (soft knee)
threshold = 0.7
ratio = 4
compressed = np.zeros_like(input_signal)
for i in range(len(input_signal)):
if envelope[i] > threshold:
gain_reduction = 1 + (envelope[i] - threshold) * (1 - 1/ratio) / envelope[i]
compressed[i] = input_signal[i] / gain_reduction
else:
compressed[i] = input_signal[i]
axes[2, 0].plot(t_comp*1000, input_signal, 'b-', alpha=0.7, label='Input')
axes[2, 0].plot(t_comp*1000, compressed, 'r-', linewidth=2, label='Compressed')
axes[2, 0].plot(t_comp*1000, envelope, 'g--', label='Envelope')
axes[2, 0].axhline(y=threshold, color='k', linestyle=':', label='Threshold')
axes[2, 0].set_xlabel('Time (ms)')
axes[2, 0].set_ylabel('Amplitude')
axes[2, 0].set_title('Dynamic Range Compression')
axes[2, 0].legend()
axes[2, 0].grid(True, alpha=0.3)
axes[2, 0].set_xlim(0, 500)
# 6. Reverb all-pass filter network
print("\nReverb Design (Schroeder Reverb):")
# All-pass filter for reverb
delay_time = 0.03 # 30 ms
g = 0.7 # Feedback gain
# Transfer function of all-pass filter
# H(s) = (g + e^(-sT)) / (1 + g*e^(-sT))
T = delay_time
# Impulse response simulation
t_reverb = np.linspace(0, 0.2, 2000)
impulse_response = np.zeros_like(t_reverb)
# Generate impulse response with multiple delays
delays = [0.03, 0.037, 0.041, 0.043] # Prime number ratios
gains = [0.7, 0.65, 0.6, 0.55]
for delay, gain in zip(delays, gains):
n_delay = int(delay * len(t_reverb) / t_reverb[-1])
temp_ir = np.zeros_like(t_reverb)
temp_ir[0] = 1
# Recursive all-pass structure
for i in range(n_delay, len(t_reverb)):
temp_ir[i] = -gain * temp_ir[i-n_delay]
if i >= n_delay:
temp_ir[i] += gain * temp_ir[i-n_delay]
impulse_response += temp_ir
# Apply exponential decay
impulse_response *= np.exp(-3 * t_reverb)
axes[2, 1].plot(t_reverb*1000, impulse_response, 'purple', linewidth=1)
axes[2, 1].set_xlabel('Time (ms)')
axes[2, 1].set_ylabel('Amplitude')
axes[2, 1].set_title('Reverb Impulse Response')
axes[2, 1].grid(True, alpha=0.3)
axes[2, 1].set_xlim(0, 200)
plt.tight_layout()
plt.show()
# Summary of transfer functions
print("\nSummary of Audio System Transfer Functions:")
print("==========================================")
print("1. Parametric EQ: Adjustable center frequency, Q, and gain")
print("2. Crossover: Splits audio into frequency bands for speakers")
print("3. Tone Stack: Shapes frequency response in guitar amps")
print("4. Room Acoustics: Models reflections and comb filtering")
print("5. Compressor: Non-linear dynamics processing")
print("6. Reverb: All-pass filters create dense reflections")Key formulas, tables, and resources for Laplace transforms.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# Create comprehensive Laplace transform table
transform_table = [
# Basic functions
["Unit step", "u(t) or H(t)", "1/s", "s > 0"],
["Unit impulse", "δ(t)", "1", "all s"],
["Ramp", "t", "1/s²", "s > 0"],
["Power", "t^n", "n!/s^(n+1)", "s > 0"],
["Exponential", "e^(at)", "1/(s-a)", "s > a"],
["Sine", "sin(ωt)", "ω/(s²+ω²)", "s > 0"],
["Cosine", "cos(ωt)", "s/(s²+ω²)", "s > 0"],
["Hyperbolic sine", "sinh(at)", "a/(s²-a²)", "s > |a|"],
["Hyperbolic cosine", "cosh(at)", "s/(s²-a²)", "s > |a|"],
# Damped functions
["Damped sine", "e^(-at)sin(ωt)", "ω/((s+a)²+ω²)", "s > -a"],
["Damped cosine", "e^(-at)cos(ωt)", "(s+a)/((s+a)²+ω²)", "s > -a"],
# Products with t
["t×exponential", "te^(at)", "1/(s-a)²", "s > a"],
["t×sine", "t·sin(ωt)", "2ωs/(s²+ω²)²", "s > 0"],
["t×cosine", "t·cos(ωt)", "(s²-ω²)/(s²+ω²)²", "s > 0"],
# Special functions
["Error function", "erf(√t)", "1/(s√(s+1))", "s > 0"],
["Bessel J₀", "J₀(t)", "1/√(s²+1)", "s > 0"],
# Discontinuous
["Shifted step", "u(t-a)", "e^(-as)/s", "s > 0"],
["Pulse", "u(t-a)-u(t-b)", "(e^(-as)-e^(-bs))/s", "s > 0"],
["Periodic square", "Square wave (T)", "(1-e^(-sT/2))/(s(1-e^(-sT)))", "s > 0"]
]
# Create DataFrame
df = pd.DataFrame(transform_table, columns=['Function Name', 'f(t)', 'F(s)', 'ROC'])
print("Laplace Transform Reference Table")
print("================================")
print(df.to_string(index=False))
# Visualize some key transform pairs
fig, axes = plt.subplots(3, 3, figsize=(15, 12))
t = np.linspace(0, 5, 1000)
# Define functions for plotting
functions = [
("Unit Step", lambda t: np.ones_like(t), "1/s"),
("Exponential e^(-2t)", lambda t: np.exp(-2*t), "1/(s+2)"),
("Sine (2t)", lambda t: np.sin(2*t), "2/(s²+4)"),
("Cosine (2t)", lambda t: np.cos(2*t), "s/(s²+4)"),
("Ramp t", lambda t: t, "1/s²"),
("t·e^(-t)", lambda t: t*np.exp(-t), "1/(s+1)²"),
("Damped sine", lambda t: np.exp(-t)*np.sin(3*t), "3/((s+1)²+9)"),
("Square pulse", lambda t: np.where((t > 1) & (t < 3), 1, 0), "(e^(-s)-e^(-3s))/s"),
("Triangular", lambda t: np.where(t < 2, t/2, np.where(t < 4, 2-t/2, 0)), "Complex")
]
for idx, (name, func, transform) in enumerate(functions):
row = idx // 3
col = idx % 3
y = func(t)
axes[row, col].plot(t, y, 'b-', linewidth=2)
axes[row, col].set_xlabel('t')
axes[row, col].set_ylabel('f(t)')
axes[row, col].set_title(f'{name}\nF(s) = {transform}')
axes[row, col].grid(True, alpha=0.3)
axes[row, col].set_xlim(0, 5)
plt.tight_layout()
plt.show()print("\nLaplace Transform Properties")
print("===========================")
properties = [
("Linearity", "L{af(t) + bg(t)} = aF(s) + bG(s)"),
("Time shifting", "L{f(t-a)u(t-a)} = e^(-as)F(s)"),
("Frequency shifting", "L{e^(at)f(t)} = F(s-a)"),
("Time scaling", "L{f(at)} = (1/a)F(s/a), a > 0"),
("Differentiation (time)", "L{f'(t)} = sF(s) - f(0)"),
("Differentiation (freq)", "L{tf(t)} = -dF/ds"),
("Integration", "L{∫f(τ)dτ} = F(s)/s"),
("Convolution", "L{f*g} = F(s)G(s)"),
("Initial value", "f(0+) = lim(s→∞) sF(s)"),
("Final value", "f(∞) = lim(s→0) sF(s)")
]
for prop, formula in properties:
print(f"{prop:25} {formula}")
# Common partial fraction patterns
print("\nPartial Fraction Patterns")
print("========================")
print("1. Distinct linear factors:")
print(" F(s) = A/(s-a) + B/(s-b) + ...")
print("\n2. Repeated linear factors:")
print(" F(s) = A₁/(s-a) + A₂/(s-a)² + ... + Aₙ/(s-a)ⁿ")
print("\n3. Quadratic factors:")
print(" F(s) = (As + B)/(s² + ps + q)")
print("\n4. Repeated quadratic factors:")
print(" F(s) = (A₁s + B₁)/(s² + ps + q) + (A₂s + B₂)/(s² + ps + q)² + ...")
# Solving differential equations procedure
print("\nSolving ODEs with Laplace Transforms")
print("===================================")
print("1. Take Laplace transform of the entire equation")
print("2. Apply initial conditions")
print("3. Solve algebraically for Y(s)")
print("4. Use partial fractions if necessary")
print("5. Apply inverse Laplace transform")
print("6. Verify solution satisfies original equation")
# Example workflow
print("\nExample: Solve y'' + 3y' + 2y = e^(-t), y(0) = 0, y'(0) = 1")
print("--------------------------------------------------------")
import sympy as sp
t, s = sp.symbols('t s', real=True)
Y = sp.Function('Y')
# Step 1: Laplace transform
print("Step 1: L{y'' + 3y' + 2y} = L{e^(-t)}")
print(" s²Y(s) - sy(0) - y'(0) + 3[sY(s) - y(0)] + 2Y(s) = 1/(s+1)")
# Step 2: Apply ICs
print("\nStep 2: Apply y(0) = 0, y'(0) = 1")
print(" s²Y(s) - 1 + 3sY(s) + 2Y(s) = 1/(s+1)")
# Step 3: Solve for Y(s)
print("\nStep 3: Y(s)[s² + 3s + 2] = 1/(s+1) + 1")
Y_s = (1/(s+1) + 1) / (s**2 + 3*s + 2)
Y_s = sp.simplify(Y_s)
print(f" Y(s) = {Y_s}")
# Step 4: Partial fractions
Y_s_pf = sp.apart(Y_s, s)
print(f"\nStep 4: Partial fractions: Y(s) = {Y_s_pf}")
# Step 5: Inverse transform
y_t = sp.inverse_laplace_transform(Y_s, s, t)
print(f"\nStep 5: y(t) = {y_t}")
# Verification plot
t_vals = np.linspace(0, 5, 1000)
y_vals = sp.lambdify(t, y_t * sp.Heaviside(t), 'numpy')(t_vals)
plt.figure(figsize=(10, 6))
plt.plot(t_vals, y_vals, 'b-', linewidth=2, label='Solution y(t)')
plt.plot(t_vals, np.exp(-t_vals), 'r--', linewidth=2, label='Forcing e^(-t)')
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Solution to y\'\'\' + 3y\' + 2y = e^(-t)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()