别再死记硬背了！用Python+Matplotlib动画演示采样定理与频谱混叠（附代码）

用Python动画拆解采样定理：从频谱混叠到奈奎斯特准则的视觉化实践

在数字信号处理领域，采样定理就像是一把打开模拟世界与数字世界大门的钥匙。但传统教材中复杂的公式推导往往让学习者望而生畏——当我们盯着那些频域搬移的数学表达式时，很难在脑海中形成直观的画面。这就是为什么我们需要换种方式学习：用代码构建信号，用动画观察频谱，让抽象理论变成可交互的视觉实验。本文将带您用Python和Matplotlib搭建一个完整的信号采样模拟系统，通过修改参数实时观察时域采样点与频域频谱的联动变化。

1. 环境配置与基础信号生成

1.1 初始化Python环境

首先确保已安装科学计算三件套：

pip install numpy matplotlib scipy

创建基础正弦波信号生成器：

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def generate_sine_wave(freq, duration=1, sample_rate=44100): t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False) signal = np.sin(2 * np.pi * freq * t) return t, signal

注意：默认采样率44100Hz是CD音质的标准，足够覆盖人耳可闻范围(20-20kHz)

1.2 多频信号合成实战

真实世界的信号很少是单一频率，让我们合成一个包含基波与谐波的复合信号：

def generate_composite_signal(frequencies, amplitudes, duration=1, sample_rate=44100): t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False) signal = np.zeros_like(t) for freq, amp in zip(frequencies, amplitudes): signal += amp * np.sin(2 * np.pi * freq * t) return t, signal

调用示例：

freqs = [50, 120, 300] # 基波50Hz + 两个谐波 amps = [0.6, 0.3, 0.1] t, signal = generate_composite_signal(freqs, amps)

2. 采样过程的可视化实现

2.1 时域采样点提取

关键函数实现采样点抽取：

def sample_signal(t, signal, sample_interval): sample_indices = np.arange(0, len(t), sample_interval) sample_times = t[sample_indices] sample_values = signal[sample_indices] return sample_times, sample_values, sample_indices

2.2 动态采样对比可视化

创建动画展示不同采样率效果：

def animate_sampling(freq=10, duration=1, max_sample_rate=100): fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6)) # 生成原始信号 t_highres = np.linspace(0, duration, 10000) signal = np.sin(2 * np.pi * freq * t_highres) line_orig, = ax1.plot(t_highres, signal, 'b-', alpha=0.3) line_samples, = ax1.plot([], [], 'ro') # 频谱初始化 fft_freq = np.fft.fftfreq(len(t_highres), t_highres[1]-t_highres[0]) fft_vals = np.abs(np.fft.fft(signal)) line_fft, = ax2.plot(fft_freq[:len(fft_freq)//2], fft_vals[:len(fft_freq)//2], 'b-') def update(frame): sample_rate = frame + 1 sample_interval = int(10000 / sample_rate) # 更新采样点 sample_indices = np.arange(0, 10000, sample_interval) line_samples.set_data(t_highres[sample_indices], signal[sample_indices]) # 更新频谱 sampled_signal = np.zeros_like(signal) sampled_signal[sample_indices] = signal[sample_indices] fft_vals = np.abs(np.fft.fft(sampled_signal)) line_fft.set_ydata(fft_vals[:len(fft_freq)//2]) ax1.set_title(f'Sampling Rate: {sample_rate}Hz | Signal Freq: {freq}Hz') return line_samples, line_fft ani = FuncAnimation(fig, update, frames=max_sample_rate, interval=200) plt.tight_layout() return ani

3. 频谱分析与混叠现象观测

3.1 FFT参数设置要点

获取准确频谱的关键参数：

def compute_fft(signal, sample_rate): n = len(signal) fft_vals = np.fft.fft(signal) fft_freq = np.fft.fftfreq(n, 1/sample_rate) # 取单边频谱 positive_freq = fft_freq[:n//2] positive_fft = 2/n * np.abs(fft_vals[:n//2]) return positive_freq, positive_fft

重要提示：FFT结果的幅度需要乘以2/N进行归一化，直流分量(N=0)除外

3.2 混叠实验设计

创建混叠对比实验：

def aliasing_experiment(signal_freq=30, sample_rates=[100, 60, 40]): duration = 1 t_highres = np.linspace(0, duration, 10000) signal = np.sin(2 * np.pi * signal_freq * t_highres) fig, axes = plt.subplots(len(sample_rates), 2, figsize=(12, 8)) for i, sr in enumerate(sample_rates): # 时域采样 interval = int(10000 / (sr * duration)) sample_indices = np.arange(0, 10000, interval) # 频域分析 sampled_signal = np.zeros_like(signal) sampled_signal[sample_indices] = signal[sample_indices] freq, fft = compute_fft(sampled_signal, sr) # 绘图 axes[i,0].plot(t_highres, signal, 'b-', alpha=0.3) axes[i,0].plot(t_highres[sample_indices], signal[sample_indices], 'ro') axes[i,0].set_title(f'Sample Rate: {sr}Hz') axes[i,1].stem(freq, fft, 'r', markerfmt=" ") axes[i,1].set_xlim(0, max(sample_rates)) axes[i,1].axvline(signal_freq, color='b', linestyle='--') # 标记混叠频率 if sr < 2 * signal_freq: alias_freq = abs(sr - signal_freq) axes[i,1].axvline(alias_freq, color='g', linestyle=':') axes[i,1].text(alias_freq+2, max(fft), f'Alias: {alias_freq}Hz', color='g') plt.tight_layout() return fig

4. 奈奎斯特准则的交互验证

4.1 动态阈值测试工具

构建可交互的采样率测试界面：

from ipywidgets import interact, FloatSlider def interactive_sampling_test(signal_freq=10): def plot_sampling(sample_rate=20.0): duration = 1 t = np.linspace(0, duration, 1000) signal = np.sin(2 * np.pi * signal_freq * t) # 采样 sample_interval = int(1000 / sample_rate) sample_indices = np.arange(0, 1000, sample_interval) # 重建信号 reconstructed = np.zeros_like(signal) for i in sample_indices: reconstructed += signal[i] * np.sinc((t - t[i]) * sample_rate) # 绘图 plt.figure(figsize=(10,4)) plt.plot(t, signal, 'b-', label='Original', alpha=0.3) plt.plot(t[sample_indices], signal[sample_indices], 'ro', label='Samples') plt.plot(t, reconstructed, 'g--', label='Reconstructed', alpha=0.7) plt.title(f'Sample Rate: {sample_rate:.1f}Hz | Nyquist: {2*signal_freq}Hz') plt.legend() plt.grid(True) interact(plot_sampling, sample_rate=FloatSlider(min=1, max=50, step=0.5, value=20))

4.2 实际应用中的采样率选择

常见场景的采样率规范：

经验法则：实际采样率通常设为最高频率的2.56-4倍，为抗混叠滤波器留出过渡带

5. 脉冲编码调制(PCM)的Python实现

5.1 量化与编码过程

def pcm_encode(signal, bits=8): # 归一化到[0,1]范围 normalized = (signal - signal.min()) / (signal.max() - signal.min()) # 量化 quantized = np.round(normalized * (2**bits - 1)) # 转换为整数类型 return quantized.astype(int) def pcm_decode(encoded, bits=8, original_min=0, original_max=1): # 反量化 normalized = encoded / (2**bits - 1) # 恢复原始范围 return normalized * (original_max - original_min) + original_min

5.2 完整PCM流程演示

def demo_pcm_workflow(freq=5, sample_rate=30, bits=4): # 生成信号 t = np.linspace(0, 1, sample_rate) analog = np.sin(2 * np.pi * freq * t) # PCM编码 digital = pcm_encode(analog, bits) # PCM解码 reconstructed = pcm_decode(digital, bits, analog.min(), analog.max()) # 绘图 plt.figure(figsize=(10,4)) plt.plot(t, analog, 'b-o', label='Analog', alpha=0.5) plt.step(t, reconstructed, 'r--', where='post', label='PCM') plt.title(f'PCM with {bits}-bit Quantization') plt.legend() plt.grid(True)

在Jupyter notebook中运行这些代码，您会看到当采样率低于奈奎斯特频率时，那些高频分量会"伪装"成低频信号出现在频谱中——这正是导致电话中声音失真的根本原因。通过调整代码中的参数，您可以直观观察到：当采样率提高到信号最高频率的两倍以上时，这些"幽灵频率"会立即消失。