别再死记硬背采样定理了！用Python+NumPy手动画出频谱混叠全过程

用Python动态演示采样定理：从频谱混叠到信号重建的视觉之旅

当你第一次在《信号与系统》课本上看到采样定理时，是否曾被那些抽象的数学公式和频域图表弄得晕头转向？"采样频率必须大于信号最高频率的两倍"——这条看似简单的规则背后，隐藏着数字信号处理最核心的奥秘。今天，我们将抛弃枯燥的理论推导，用Python代码搭建一个可视化实验室，让你亲眼见证频谱混叠如何发生，以及如何通过调整采样率来避免它。

1. 搭建信号处理的可视化环境

在开始实验前，我们需要配置好Python科学计算的三件套：NumPy用于数值运算，Matplotlib用于可视化，SciPy提供信号处理工具包。打开你的Jupyter Notebook或Colab，运行以下代码块初始化环境：

import numpy as np import matplotlib.pyplot as plt from scipy import signal plt.style.use('seaborn') # 使用更美观的绘图样式

为了更直观地观察信号变化，我们创建一个交互式绘图函数。这个函数将同时显示时域波形和频域频谱，形成"信号显微镜"：

def plot_signal_comparison(original_signal, sampled_signal, t_original, t_sampled, fs, title): fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8)) # 时域波形对比 ax1.plot(t_original, original_signal, 'b-', label='原始信号') ax1.stem(t_sampled, sampled_signal, 'r-', markerfmt='ro', basefmt=" ", label='采样点') ax1.set_xlabel('时间 (s)') ax1.set_ylabel('幅值') ax1.legend() # 频域频谱对比 f_original, X_original = signal.periodogram(original_signal, fs) f_sampled, X_sampled = signal.periodogram(sampled_signal, fs) ax2.semilogy(f_original, X_original, 'b-', label='原始频谱') ax2.semilogy(f_sampled, X_sampled, 'r-', label='采样后频谱') ax2.set_xlabel('频率 (Hz)') ax2.set_ylabel('功率谱密度') ax2.legend() plt.suptitle(title) plt.tight_layout() plt.show()

2. 模拟不同采样率下的信号采集

让我们生成一个简单的正弦波作为测试信号，频率设为5Hz。这个选择是为了后续清晰展示奈奎斯特频率的概念：

f_signal = 5 # 信号频率(Hz) duration = 1 # 信号持续时间(s) t_continuous = np.linspace(0, duration, 5000) # 高密度时间点模拟连续信号 continuous_signal = np.sin(2 * np.pi * f_signal * t_continuous)

2.1 满足奈奎斯特准则的采样

根据采样定理，对于5Hz的信号，采样频率至少需要10Hz。我们先尝试用15Hz采样，这明显高于奈奎斯特频率：

fs_adequate = 15 # 足够的采样频率 t_samples_adequate = np.arange(0, duration, 1/fs_adequate) sampled_signal_adequate = np.sin(2 * np.pi * f_signal * t_samples_adequate) plot_signal_comparison(continuous_signal, sampled_signal_adequate, t_continuous, t_samples_adequate, fs_adequate, "满足奈奎斯特准则的采样 (fs=15Hz > 2×5Hz)")

观察频谱图，你会看到：

• 原始信号在5Hz处有一个清晰的峰值
• 采样后的频谱保持了原始信号的频率成分
• 没有出现其他频率的干扰成分

2.2 临界采样情况

现在将采样率降到刚好等于奈奎斯特频率的10Hz：

fs_critical = 10 # 临界采样频率 t_samples_critical = np.arange(0, duration, 1/fs_critical) sampled_signal_critical = np.sin(2 * np.pi * f_signal * t_samples_critical) plot_signal_comparison(continuous_signal, sampled_signal_critical, t_continuous, t_samples_critical, fs_critical, "临界采样情况 (fs=10Hz = 2×5Hz)")

此时频谱显示：

• 仍然能够识别原始信号频率
• 但频谱峰变得更宽，分辨率下降
• 任何微小的采样时间抖动都可能导致信号无法恢复

2.3 频谱混叠的诞生

最后，我们故意违反采样定理，使用7Hz的采样率（低于10Hz的奈奎斯特频率）：

fs_aliasing = 7 # 不足的采样频率 t_samples_aliasing = np.arange(0, duration, 1/fs_aliasing) sampled_signal_aliasing = np.sin(2 * np.pi * f_signal * t_samples_aliasing) plot_signal_comparison(continuous_signal, sampled_signal_aliasing, t_continuous, t_samples_aliasing, fs_aliasing, "频谱混叠现象 (fs=7Hz < 2×5Hz)")

这个结果令人惊讶：

• 采样点连成的波形看起来是一个2Hz的低频信号（7Hz-5Hz）
• 频谱图上出现了5Hz原始频率和2Hz的混叠频率
• 这就是为什么欠采样会导致信号失真的直观证明

3. 混叠现象的深入解析

频谱混叠不是简单的数据丢失，而是频率成分的"伪装"。当采样率不足时，高频信号会"假扮"成低频信号混入基带。这种现象可以通过以下公式解释：

f_alias = |f_signal - n × fs|

其中n是使f_alias落在0到fs/2范围内的整数。在我们的例子中：

f_alias = |5 - 1×7| = 2Hz

为了更全面地观察这一现象，我们可以创建一个交互式可视化工具，实时显示不同采样率下的效果：

from ipywidgets import interact, FloatSlider def interactive_aliasing_demo(f_signal=5, fs=15): t = np.linspace(0, 1, 1000) signal_cont = np.sin(2*np.pi*f_signal*t) sample_times = np.arange(0, 1, 1/fs) signal_samples = np.sin(2*np.pi*f_signal*sample_times) plt.figure(figsize=(12, 6)) plt.plot(t, signal_cont, 'b-', alpha=0.5, label='连续信号') plt.stem(sample_times, signal_samples, 'r-', markerfmt='ro', basefmt=" ", label='采样点') if fs < 2*f_signal: f_alias = abs(f_signal - fs*round(f_signal/fs)) plt.title(f'发生混叠！{f_signal}Hz信号被采样为{f_alias:.1f}Hz (fs={fs}Hz)') else: plt.title(f'正确采样：{f_signal}Hz信号被准确采样 (fs={fs}Hz)') plt.xlabel('时间 (s)') plt.ylabel('幅值') plt.legend() plt.grid(True) plt.show() interact(interactive_aliasing_demo, f_signal=FloatSlider(min=1, max=20, step=0.5, value=5), fs=FloatSlider(min=1, max=30, step=0.5, value=15))

通过拖动滑块，你可以实时观察到：

• 当fs > 2×f_signal时，采样点准确重建原始波形
• 当fs < 2×f_signal时，采样点形成完全不同的低频波形
• 混叠频率随采样率变化而跳变的规律

4. 从采样到重建：完整的信号处理流程

理解了混叠现象后，让我们完成信号处理的闭环——从采样到重建。我们将演示如何通过理想低通滤波器从采样数据中恢复原始信号。

4.1 采样信号的频谱特性

采样过程在数学上等价于原始信号与脉冲序列的乘积。在频域中，这表现为原始频谱的周期性复制：

def plot_sampling_spectrum(f_signal=5, fs=15): # 生成信号 t = np.linspace(0, 1, 10000) signal_cont = np.sin(2*np.pi*f_signal*t) # 计算频谱 f, Pxx = signal.periodogram(signal_cont, fs=1000) # 绘制 plt.figure(figsize=(12, 6)) plt.semilogy(f, Pxx, 'b-', label='原始信号频谱') # 标记奈奎斯特频率 nyquist = fs/2 plt.axvline(x=nyquist, color='r', linestyle='--', label=f'奈奎斯特频率 ({nyquist}Hz)') # 标记混叠频率 if fs < 2*f_signal: f_alias = abs(f_signal - fs*round(f_signal/fs)) plt.axvline(x=f_alias, color='g', linestyle=':', label=f'混叠频率 ({f_alias}Hz)') plt.title(f'信号频谱分析 (f={f_signal}Hz, fs={fs}Hz)') plt.xlabel('频率 (Hz)') plt.ylabel('功率谱密度') plt.legend() plt.grid(True) plt.xlim(0, 30) plt.show() interact(plot_sampling_spectrum, f_signal=FloatSlider(min=1, max=20, step=0.5, value=5), fs=FloatSlider(min=1, max=30, step=0.5, value=15))

4.2 信号重建的滤波器设计

要从采样信号中恢复原始信号，我们需要一个理想低通滤波器，其截止频率为fs/2。虽然现实中不存在完美的滤波器，但我们可以用数字滤波器近似实现：

def reconstruct_signal(sampled_signal, fs, cutoff_factor=0.8): # 设计低通滤波器 nyquist = fs / 2 cutoff = cutoff_factor * nyquist b, a = signal.butter(4, cutoff/(fs/2), 'low') # 上采样 upsampling_factor = 20 upsampled = signal.resample(sampled_signal, len(sampled_signal)*upsampling_factor) # 滤波 reconstructed = signal.filtfilt(b, a, upsampled) return reconstructed # 示例：重建欠采样信号 reconstructed_aliasing = reconstruct_signal(sampled_signal_aliasing, fs_aliasing) t_reconstructed = np.linspace(0, duration, len(reconstructed_aliasing)) plt.figure(figsize=(12, 6)) plt.plot(t_continuous, continuous_signal, 'b-', label='原始信号') plt.plot(t_reconstructed, reconstructed_aliasing, 'r--', label='重建信号') plt.title("欠采样信号的重建尝试 (fs=7Hz < 2×5Hz)") plt.xlabel('时间 (s)') plt.ylabel('幅值') plt.legend() plt.grid(True) plt.show()

这个实验清楚地展示了：

• 即使使用复杂的重建算法，欠采样导致的混叠也无法消除
• 重建出的信号是混叠频率(2Hz)而非原始信号(5Hz)
• 这验证了采样定理不是理论游戏，而是数字信号处理的根本限制

5. 实际工程中的采样策略

理解了基本原理后，我们需要讨论实际工程中如何应用采样定理。以下是几个关键考虑因素：

在真实系统中，采样前必须使用抗混叠滤波器限制信号带宽。例如，音频CD采用44.1kHz采样率处理20kHz音频信号，这提供了：

安全余量 = (44.1 - 2×20)/20 = 10.5%

这种设计考虑了：

• 模拟滤波器的滚降特性
• 采样时钟的微小抖动
• 后续数字处理的灵活性

以下是一个模拟完整采样系统的Python示例：

def full_sampling_system(f_signal=5, fs=15, filter_order=4, ripple_db=1): # 1. 生成原始信号(含高频成分模拟真实信号) t = np.linspace(0, 1, 10000) signal_cont = np.sin(2*np.pi*f_signal*t) + 0.2*np.sin(2*np.pi*3*f_signal*t) # 2. 模拟抗混叠滤波器 nyquist = fs/2 cutoff = 0.9 * nyquist # 预留10%过渡带 b, a = signal.butter(filter_order, cutoff/(10000/2), 'low') filtered_signal = signal.filtfilt(b, a, signal_cont) # 3. 采样 sample_indices = (np.arange(0, 1, 1/fs) * 10000).astype(int) sampled_signal = filtered_signal[sample_indices] sample_times = t[sample_indices] # 4. 重建 reconstructed = reconstruct_signal(sampled_signal, fs) t_reconstructed = np.linspace(0, 1, len(reconstructed)) # 绘图 plt.figure(figsize=(12, 8)) plt.plot(t, signal_cont, 'b-', alpha=0.3, label='原始信号') plt.plot(t, filtered_signal, 'g-', alpha=0.5, label='抗混叠滤波后') plt.stem(sample_times, sampled_signal, 'r-', markerfmt='ro', basefmt=" ", label='采样点') plt.plot(t_reconstructed, reconstructed, 'm--', linewidth=2, label='重建信号') title = f'完整采样系统模拟 (f={f_signal}Hz, fs={fs}Hz)' if fs < 2*f_signal: f_alias = abs(f_signal - fs*round(f_signal/fs)) title += f' - 混叠至{f_alias:.1f}Hz' plt.title(title) plt.xlabel('时间 (s)') plt.ylabel('幅值') plt.legend() plt.grid(True) plt.show() interact(full_sampling_system, f_signal=FloatSlider(min=1, max=20, step=0.5, value=5), fs=FloatSlider(min=1, max=30, step=0.5, value=15), filter_order=[2, 4, 6, 8], ripple_db=FloatSlider(min=0.1, max=5, step=0.1, value=1))

通过这个完整模拟，你会发现：

• 即使采样率略低于2×f_max，好的抗混叠滤波器也能减轻混叠
• 滤波器阶数越高，过渡带越陡峭，但会引入相位失真
• 工程实践总是在各种因素间寻找平衡点