用Python破解ADC采样迷思：从混叠现象到信号重建的代码实践

在信号处理领域，奈奎斯特定理就像是一把双刃剑——它既是指引采样过程的明灯，也是许多工程师的噩梦。传统教材中复杂的数学推导和抽象的频谱图，往往让学习者陷入公式记忆的泥潭，却难以建立直观理解。本文将带你用Python代码"看见"采样过程，通过可交互的视觉化演示，让那些困扰你的混叠、重建概念变得触手可及。

1. 搭建你的数字信号实验室

1.1 环境配置与基础波形生成

任何实验都需要一个可靠的实验室，我们的数字信号处理实验从配置Python环境开始。确保安装了以下核心库：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

让我们首先生成一个基础的正弦波信号——这是理解采样原理的最佳起点。下面的代码创建了一个1kHz的纯净正弦波，采样率为100kHz（模拟连续信号）：

fs_continuous = 100e3  # 100kHz的"伪连续"采样率
t_continuous = np.arange(0, 0.01, 1/fs_continuous)  # 10ms时间序列
f_signal = 1e3  # 1kHz信号频率
amplitude = 1
continuous_signal = amplitude * np.sin(2 * np.pi * f_signal * t_continuous)

为什么选择100kHz作为"伪连续"采样率？ 因为对于1kHz的信号来说，这已经足够密集，可以近似看作连续信号。我们将在后续步骤中对比不同采样率下的效果。

1.2 可视化基础信号

理解信号的第一步是观察它的时域和频域表现：

plt.figure(figsize=(12, 4))
plt.subplot(121)
plt.plot(t_continuous[:200], continuous_signal[:200])  # 只显示前200个点
plt.title('时域表现 (1kHz正弦波)')
plt.xlabel('时间(s)')
plt.ylabel('幅度')

plt.subplot(122)
freq = np.fft.fftfreq(len(continuous_signal), 1/fs_continuous)
fft_result = np.abs(np.fft.fft(continuous_signal))
plt.plot(freq[:len(freq)//2], fft_result[:len(freq)//2])
plt.title('频域表现')
plt.xlabel('频率(Hz)')
plt.ylabel('幅度')
plt.tight_layout()
plt.show()

这段代码会生成一个双面板图：左侧显示信号的时域波形，右侧显示其频谱。纯净正弦波在频谱上应该只显示单一的1kHz峰值。

2. 采样率与混叠现象的直观演示

2.1 实施不同采样率的采样过程

现在，让我们模拟真实ADC的采样过程。我们将尝试三种不同的采样率：远高于奈奎斯特频率、等于奈奎斯特频率、以及低于奈奎斯特频率的情况。

def sample_signal(original_signal, original_fs, new_fs):
    """对信号进行采样率转换"""
    ratio = int(original_fs / new_fs)
    sampled_signal = original_signal[::ratio]
    sampled_time = t_continuous[::ratio]
    return sampled_time, sampled_signal

# 三种采样率场景
fs_high = 10e3  # 10kHz (>> 2kHz Nyquist)
fs_nyquist = 2e3  # 2kHz (刚好满足Nyquist)
fs_low = 1.5e3  # 1.5kHz (< Nyquist)

2.2 可视化混叠现象

关键的时刻到了——让我们看看不同采样率下会发生什么：

plt.figure(figsize=(15, 10))

# 高采样率情况
t_high, sig_high = sample_signal(continuous_signal, fs_continuous, fs_high)
plt.subplot(321)
plt.plot(t_continuous, continuous_signal, 'b-', alpha=0.3)
plt.stem(t_high, sig_high, 'r', markerfmt='ro', basefmt=" ")
plt.title(f'高采样率 ({fs_high/1e3}kHz) 时域')

# 奈奎斯特采样率情况
t_nyq, sig_nyq = sample_signal(continuous_signal, fs_continuous, fs_nyquist)
plt.subplot(323)
plt.plot(t_continuous, continuous_signal, 'b-', alpha=0.3)
plt.stem(t_nyq, sig_nyq, 'r', markerfmt='ro', basefmt=" ")
plt.title(f'奈奎斯特采样率 ({fs_nyquist/1e3}kHz) 时域')

# 低采样率情况
t_low, sig_low = sample_signal(continuous_signal, fs_continuous, fs_low)
plt.subplot(325)
plt.plot(t_continuous, continuous_signal, 'b-', alpha=0.3)
plt.stem(t_low, sig_low, 'r', markerfmt='ro', basefmt=" ")
plt.title(f'低采样率 ({fs_low/1e3}kHz) 时域')

# 频谱分析
def plot_spectrum(signal, fs, subplot_num, title):
    n = len(signal)
    freq = np.fft.fftfreq(n, 1/fs)
    fft_result = np.abs(np.fft.fft(signal))
    plt.subplot(subplot_num)
    plt.stem(freq[:n//2], fft_result[:n//2], markerfmt=' ')
    plt.title(title)
    plt.xlabel('频率(Hz)')
    plt.ylabel('幅度')

plot_spectrum(sig_high, fs_high, 322, '高采样率频谱')
plot_spectrum(sig_nyq, fs_nyquist, 324, '奈奎斯特采样率频谱')
plot_spectrum(sig_low, fs_low, 326, '低采样率频谱')

plt.tight_layout()
plt.show()

这个对比图会清晰地展示：当采样率低于奈奎斯特频率时，频谱中会出现虚假的频率成分——这就是混叠现象。特别有趣的是，1.5kHz采样率下，原本1kHz的信号会"伪装"成500Hz的信号。

3. 信号重建的艺术与技术

3.1 零阶保持重建

实际DAC中最常用的重建方法是零阶保持(ZOH)。让我们实现一个简单的ZOH重建：

def zoh_reconstruction(sampled_time, sampled_signal, output_fs):
    """零阶保持重建"""
    reconstruction = np.zeros(int(sampled_time[-1] * output_fs))
    time_vector = np.arange(len(reconstruction)) / output_fs
    
    sample_indices = (sampled_time * output_fs).astype(int)
    for i in range(len(sample_indices)-1):
        start = sample_indices[i]
        end = sample_indices[i+1]
        reconstruction[start:end] = sampled_signal[i]
    
    return time_vector, reconstruction

output_fs = fs_continuous  # 重建到原始"连续"采样率
time_zoh, signal_zoh = zoh_reconstruction(t_high, sig_high, output_fs)

3.2 理想重建与ZOH对比

理想重建理论上需要sinc函数插值，而ZOH则是一种实用的近似：

plt.figure(figsize=(12, 6))

# 理想重建（使用scipy的resample函数）
num_samples = int(t_high[-1] * output_fs)
signal_ideal = signal.resample(sig_high, num_samples)

plt.plot(t_continuous, continuous_signal, 'b-', label='原始信号', alpha=0.3)
plt.plot(time_zoh, signal_zoh, 'r-', label='ZOH重建')
plt.plot(time_zoh, signal_ideal, 'g--', label='理想重建')
plt.title('不同重建方法对比')
plt.xlabel('时间(s)')
plt.ylabel('幅度')
plt.legend()
plt.grid()
plt.show()

这个对比清晰地展示了ZOH重建的阶梯特性与理想重建的平滑曲线之间的差异。在实际工程中，我们会在ZOH之后使用模拟低通滤波器来平滑这些阶梯。

4. 量化效应与噪声分析

4.1 模拟ADC的量化过程

量化是ADC的另一个核心环节。让我们模拟一个3位ADC的量化过程：

def quantize(signal, bits):
    """模拟量化过程"""
    max_val = np.max(np.abs(signal))
    quantized = np.round(signal / (2*max_val) * (2**bits - 1)) 
    quantized = quantized / (2**bits - 1) * (2*max_val)
    return quantized

# 对采样后的信号进行量化
bits = 3
quantized_signal = quantize(sig_high, bits)

4.2 量化误差分析

量化会引入误差，我们可以计算并可视化这种误差：

quantization_error = sig_high - quantized_signal

plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.stem(t_high, sig_high, 'b', markerfmt='bo', label='原始采样', basefmt=" ")
plt.step(t_high, quantized_signal, 'r', where='post', label='量化后')
plt.title('3位ADC量化效果')
plt.legend()

plt.subplot(212)
plt.stem(t_high, quantization_error, 'g', markerfmt='go', basefmt=" ")
plt.title('量化误差')
plt.xlabel('时间(s)')
plt.ylabel('误差幅度')

plt.tight_layout()
plt.show()

# 计算SQNR
signal_power = np.mean(sig_high**2)
noise_power = np.mean(quantization_error**2)
sqnr = 10 * np.log10(signal_power / noise_power)
print(f'实测SQNR: {sqnr:.2f} dB')
print(f'理论SQNR (6.02B + 1.76): {6.02*bits + 1.76:.2f} dB')

这个分析展示了量化如何引入误差，以及如何计算信号量化噪声比(SQNR)。你会发现实测结果与理论公式6.02B + 1.76非常接近。

5. 综合实验：从采样到重建的完整流程

现在，让我们把这些知识整合到一个完整的ADC/DAC模拟流程中：

# 1. 原始信号
fs_signal = 5e3  # 5kHz信号
t_signal = np.arange(0, 0.01, 1/fs_continuous)
signal = 0.5 * np.sin(2 * np.pi * 1e3 * t_signal) + 0.3 * np.sin(2 * np.pi * 3e3 * t_signal)

# 2. 抗混叠滤波 (模拟实际ADC前的滤波器)
b, a = signal.butter(4, 2e3, 'lowpass', fs=fs_continuous)
filtered_signal = signal.filtfilt(b, a, signal)

# 3. 采样
fs_adc = 5e3  # 5kHz采样率
t_samples, sampled = sample_signal(filtered_signal, fs_continuous, fs_adc)

# 4. 量化
quantized = quantize(sampled, 8)  # 8位ADC

# 5. DAC重建 (ZOH)
t_dac, dac_output = zoh_reconstruction(t_samples, quantized, fs_continuous)

# 6. 重建滤波
b_recon, a_recon = signal.butter(4, 2e3, 'lowpass', fs=fs_continuous)
reconstructed = signal.filtfilt(b_recon, a_recon, dac_output)

# 可视化整个流程
plt.figure(figsize=(15, 10))
plt.subplot(411)
plt.plot(t_signal, signal)
plt.title('原始信号 (1kHz + 3kHz)')

plt.subplot(412)
plt.plot(t_signal, filtered_signal)
plt.title('抗混叠滤波后')

plt.subplot(413)
plt.stem(t_samples, sampled, 'r', markerfmt='ro', basefmt=" ")
plt.plot(t_samples, quantized, 'g', linestyle='steps-post')
plt.title('采样与量化 (红色=采样, 绿色=量化后)')

plt.subplot(414)
plt.plot(t_dac, dac_output, 'r', alpha=0.5)
plt.plot(t_signal, reconstructed, 'b')
plt.title('DAC输出 (红色) 与重建滤波后 (蓝色)')

plt.tight_layout()
plt.show()

这个完整的流程展示了实际ADC/DAC系统中的每个关键环节：抗混叠滤波、采样、量化、零阶保持重建以及最终的平滑滤波。通过调整代码中的参数，你可以直观地观察每个环节对最终信号质量的影响。