别再死记硬背了！用Python+NumPy可视化理解冲激函数如何‘抓取’信号采样点

用Python+NumPy可视化冲激函数：信号采样的数学魔术揭秘

冲激函数（δ函数）在信号处理教材里总是披着数学符号的外衣，让初学者望而生畏。但当我们用Python将其可视化时，这个抽象概念会立刻变得鲜活起来——它本质上是个"数学镊子"，能精准抓取信号在特定时刻的瞬时值。本文将用NumPy和Matplotlib，带您亲手实现从连续信号中"夹取"采样点的完整过程。

1. 冲激函数的可视化破冰

1.1 从理论到屏幕：绘制理想冲激函数

真正的δ函数在数学上是宽度趋近于零、面积保持1的极限概念。在数字世界中，我们用窄脉冲来近似：

import numpy as np import matplotlib.pyplot as plt def plot_impulse(t0=0, width=0.1): t = np.linspace(-1, 1, 1000) impulse = np.where((t >= t0-width/2) & (t <= t0+width/2), 1/width, 0) plt.figure(figsize=(10,4)) plt.stem(t, impulse, use_line_collection=True) plt.title(f"近似冲激函数 @ t={t0} (宽度={width})") plt.ylim(0, 1.2*impulse.max()) plt.grid(True) plt.show() plot_impulse(width=0.05) # 尝试调整width观察变化

关键参数对比：

1.2 冲激函数的采样特性验证

让我们用代码验证δ函数最神奇的特性——筛选积分：

def test_sifting(signal_func, t0=0): t = np.linspace(-2, 2, 10000) width = 0.001 # 极窄的脉冲 # 生成信号和冲激 signal = signal_func(t) impulse = np.where(abs(t-t0)<=width/2, 1/width, 0) # 数值积分 product = signal * impulse integral = np.trapz(product, t) # 可视化 plt.figure(figsize=(12,4)) plt.subplot(131); plt.plot(t, signal); plt.title("原始信号") plt.subplot(132); plt.stem(t, impulse); plt.title("冲激函数") plt.subplot(133); plt.plot(t, product); plt.title(f"乘积 (积分={integral:.3f})") plt.tight_layout() plt.show() test_sifting(lambda t: np.exp(-t**2), t0=0.5) # 高斯函数测试

2. 动态采样过程的可视化

2.1 构建可交互的采样演示

使用IPython的交互控件，创建动态观察窗口：

from ipywidgets import interact def dynamic_sampling(freq=1, sample_rate=10): duration = 2 # 秒 t = np.linspace(0, duration, 5000) signal = np.sin(2*np.pi*freq*t) # 正弦波 # 生成采样点 sample_times = np.arange(0, duration, 1/sample_rate) samples = np.sin(2*np.pi*freq*sample_times) # 绘制 plt.figure(figsize=(12,5)) plt.plot(t, signal, label='连续信号') plt.stem(sample_times, samples, linefmt='r-', markerfmt='ro', basefmt=" ", use_line_collection=True, label='采样点') plt.title(f"信号采样 (f={freq}Hz, 采样率={sample_rate}Hz)") plt.legend() plt.grid(True) plt.show() interact(dynamic_sampling, freq=(0.1, 5, 0.1), sample_rate=(1, 30, 1))

采样率与信号频率的关系：

• 欠采样：当采样率 < 2×信号频率时，会出现混叠现象
• 临界采样：采样率 = 2×信号频率（奈奎斯特极限）
• 过采样：采样率 >> 2×信号频率，保留更多细节

2.2 多信号类型采样对比

不同信号在相同采样策略下的表现：

signals = { "正弦波": lambda t: np.sin(2*np.pi*2*t), "方波": lambda t: np.sign(np.sin(2*np.pi*2*t)), "锯齿波": lambda t: 2*(t % 0.5)/0.5 - 1, "噪声信号": lambda t: np.random.normal(0, 0.3, len(t)) } def compare_sampling(signal_type="正弦波", sample_rate=10): t = np.linspace(0, 1, 1000) signal = signals[signal_type](t) sample_t = np.arange(0, 1, 1/sample_rate) sample_y = signals[signal_type](sample_t) plt.figure(figsize=(12,5)) plt.plot(t, signal, label='原始信号') plt.stem(sample_t, sample_y, linefmt='r-', markerfmt='ro', use_line_collection=True, label='采样点') plt.title(f"{signal_type}采样 (rate={sample_rate}Hz)") plt.legend() plt.grid(True) plt.show() interact(compare_sampling, signal_type=list(signals.keys()), sample_rate=(5, 30, 1))

3. 从数学到代码：采样定理的实践

3.1 冲激串的生成与运用

实际采样使用的是周期性冲激串（Dirac comb）：

def dirac_comb(period, duration=2, width=0.01): t = np.linspace(-duration/2, duration/2, 10000) comb = np.zeros_like(t) for n in range(-int(duration/period), int(duration/period)+1): comb += np.where(abs(t - n*period) < width/2, 1/width, 0) return t, comb t, comb = dirac_comb(0.2) # 每0.2秒一个冲激 plt.figure(figsize=(12,3)) plt.plot(t, comb) plt.title("冲激串 (采样间隔=0.2s)") plt.grid(True) plt.show()

3.2 完整采样过程模拟

将信号与冲激串相乘实现采样：

def full_sampling_process(signal_func, sample_rate): duration = 2 t_cont = np.linspace(0, duration, 10000) # 连续时间轴 signal = signal_func(t_cont) # 生成冲激串 sample_interval = 1/sample_rate t_impulses = np.arange(0, duration, sample_interval) impulse_train = np.zeros_like(t_cont) for ti in t_impulses: idx = np.argmin(abs(t_cont - ti)) impulse_train[idx] = 1/sample_interval # 采样结果 sampled = signal * impulse_train # 绘图 plt.figure(figsize=(12,8)) plt.subplot(311); plt.plot(t_cont, signal); plt.title("原始信号") plt.subplot(312); plt.stem(t_cont, impulse_train); plt.title("冲激串") plt.subplot(313); plt.stem(t_cont, sampled); plt.title("采样结果") plt.tight_layout() plt.show() full_sampling_process(lambda t: np.sin(2*np.pi*3*t) + 0.5*np.sin(2*np.pi*7*t), sample_rate=20)

采样率选择建议表：

4. 常见问题与实战技巧

4.1 数值实现的注意事项

在代码中模拟δ函数时容易遇到的坑：

• 脉冲宽度太窄：会导致数值不稳定，建议保持width > 10*dt
• 采样时间不对齐：离散时间轴可能错过精确的采样点
• 积分误差累积：使用np.trapz比简单求和更准确

改进版的冲激函数生成器：

def robust_impulse(t, t0, width=None): if width is None: width = 3*(t[1]-t[0]) # 自动设置合理宽度 impulse = np.zeros_like(t) idx = np.argmin(np.abs(t - t0)) # 找到最近的点 left = max(0, idx - int(width/2/(t[1]-t[0]))) right = min(len(t), idx + int(width/2/(t[1]-t[0]))) impulse[left:right] = 1/(t[right]-t[left]) return impulse

4.2 信号重建实验

从采样点重建原始信号（带sinc插值）：

def reconstruct(samples, sample_times, t): """使用sinc函数重建信号""" reconstructed = np.zeros_like(t) T = sample_times[1] - sample_times[0] # 采样间隔 for n, (tn, sn) in enumerate(zip(sample_times, samples)): reconstructed += sn * np.sinc((t - tn)/T) return reconstructed # 测试重建 t_highres = np.linspace(0, 1, 5000) original = np.sin(2*np.pi*5*t_highres) + 0.3*np.cos(2*np.pi*12*t_highres) sample_rate = 30 # 必须大于2×12=24 sample_t = np.arange(0, 1, 1/sample_rate) sample_y = original[::int(len(t_highres)/len(sample_t))] reconstructed = reconstruct(sample_y, sample_t, t_highres) plt.figure(figsize=(12,5)) plt.plot(t_highres, original, label='原始信号') plt.stem(sample_t, sample_y, linefmt='r-', markerfmt='ro', use_line_collection=True, label='采样点') plt.plot(t_highres, reconstructed, 'g--', lw=2, label='重建信号') plt.legend() plt.title(f"信号重建 (采样率={sample_rate}Hz)") plt.grid(True) plt.show()

实际项目中，当采样率接近奈奎斯特极限时，需要添加抗混叠滤波器。一个简单的RC滤波器实现：

def rc_filter(signal, sample_rate, cutoff_freq): """一阶RC低通滤波器""" alpha = 1 / (1 + 1/(2*np.pi*cutoff_freq/sample_rate)) filtered = np.zeros_like(signal) filtered[0] = signal[0] for i in range(1, len(signal)): filtered[i] = alpha * signal[i] + (1 - alpha) * filtered[i-1] return filtered