信号采样入门：用Python和NumPy手把手模拟冲激函数采样（附完整代码）-编程实验室

信号采样实战：用Python代码模拟冲激函数采样全过程

在数字信号处理领域，采样是将连续时间信号转换为离散时间信号的关键步骤。而冲激函数作为理论基石，其数学抽象性常常让初学者感到困惑。本文将通过Python代码实现，带您直观理解采样原理，把抽象公式转化为可运行的实验。

1. 理解冲激函数的本质

冲激函数（δ函数）在数学上被定义为在原点处无限高、无限窄且面积为1的特殊函数。虽然现实中不存在完美的冲激信号，但它的理论模型对采样分析至关重要。

让我们先用NumPy构造一个近似的冲激序列：

import numpy as np import matplotlib.pyplot as plt def impulse(n, k=0): """生成离散冲激序列""" return np.where(n == k, 1.0, 0.0) n = np.arange(-5, 6) # 时间轴 impulse_seq = impulse(n) plt.stem(n, impulse_seq) plt.title('Discrete Impulse Function') plt.xlabel('n') plt.ylabel('δ[n]') plt.grid(True) plt.show()

这段代码生成的图形会显示在n=0处值为1，其余位置为0的离散序列。这就是离散冲激函数δ[n]的可视化表现。

注意：离散冲激函数与连续冲激函数的区别在于幅值——离散情况下为1，而连续情况下理论上是无穷大（实际实现时用极大值近似）

2. 连续信号的采样模拟

采样过程本质上是将连续信号与冲激串相乘。让我们模拟对一个正弦波进行理想采样的过程：

def continuous_signal(t): """定义连续时间信号""" return np.sin(2 * np.pi * 0.5 * t) # 0.5Hz正弦波 # 采样参数 fs = 5 # 采样频率(Hz) T = 1/fs # 采样间隔 duration = 5 # 信号时长(s) # 生成时间轴 t_continuous = np.linspace(0, duration, 1000) # 高密度点模拟连续 t_samples = np.arange(0, duration, T) # 采样时刻 # 生成信号 x_continuous = continuous_signal(t_continuous) x_samples = continuous_signal(t_samples) # 绘制结果 plt.figure(figsize=(10,4)) plt.plot(t_continuous, x_continuous, 'b-', label='Continuous Signal') plt.stem(t_samples, x_samples, 'r-', markerfmt='ro', basefmt=" ", label='Samples') plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.legend() plt.grid(True) plt.title('Signal Sampling Demonstration') plt.show()

这段代码展示了：

蓝色曲线代表原始连续信号
红色圆点代表采样结果
采样间隔T=0.2秒（对应fs=5Hz）

3. 冲激串采样的数学实现

理论上的采样是通过冲激串与信号相乘实现的。让我们用代码精确模拟这一数学过程：

def impulse_train(t, fs): """生成冲激串""" samples = np.zeros_like(t) sample_indices = (t * fs).astype(int) valid_indices = (sample_indices >= 0) & (sample_indices < len(samples)) samples[sample_indices[valid_indices]] = 1 return samples # 生成冲激串 impulse_str = impulse_train(t_continuous, fs) # 数学上的采样过程 sampled_signal = x_continuous * impulse_str # 可视化 plt.figure(figsize=(12,8)) plt.subplot(3,1,1) plt.plot(t_continuous, x_continuous) plt.title('Original Continuous Signal') plt.grid(True) plt.subplot(3,1,2) plt.stem(t_continuous, impulse_str, markerfmt='ro', linefmt='r-', basefmt=" ") plt.title('Impulse Train (Sampling Function)') plt.grid(True) plt.subplot(3,1,3) plt.stem(t_continuous, sampled_signal, markerfmt='go', linefmt='g-', basefmt=" ") plt.title('Sampled Signal (Product of Signal and Impulse Train)') plt.grid(True) plt.tight_layout() plt.show()

这个实现清晰地展示了：

原始信号与冲激串相乘的过程
采样结果只在冲激时刻保留信号值
其他时刻乘积为零

4. 采样定理的验证实验

奈奎斯特采样定理指出，采样频率必须大于信号最高频率的两倍才能避免混叠。让我们通过实验验证这一点：

def plot_sampling_effect(f_signal, f_sample): """观察不同采样频率下的效果""" t = np.linspace(0, 3, 1000) x = np.sin(2 * np.pi * f_signal * t) # 采样时刻 sample_t = np.arange(0, 3, 1/f_sample) sample_x = np.sin(2 * np.pi * f_signal * sample_t) # 重建信号（线性插值） reconstructed = np.interp(t, sample_t, sample_x) plt.figure(figsize=(10,6)) plt.plot(t, x, 'b-', lw=2, label='Original') plt.stem(sample_t, sample_x, 'r-', markerfmt='ro', basefmt=" ", label='Samples') plt.plot(t, reconstructed, 'g--', lw=1.5, label='Reconstructed') plt.title(f'Signal Frequency: {f_signal}Hz, Sampling Frequency: {f_sample}Hz') plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.legend() plt.grid(True) plt.show() # 案例1：满足奈奎斯特条件 (fs > 2*f) plot_sampling_effect(f_signal=2, f_sample=5) # 案例2：不满足奈奎斯特条件 (fs < 2*f) plot_sampling_effect(f_signal=5, f_sample=8) # 看似满足，但实际有高频分量时会出现问题 # 案例3：明显不满足条件 plot_sampling_effect(f_signal=10, f_sample=15)

通过这三个案例，可以直观观察到：

当采样频率足够高时，重建信号与原始信号几乎一致
当采样频率不足时，会出现频率混叠现象
重建信号会表现出低于原始信号的频率成分

5. 实际应用中的采样技巧

在实际工程中，理想采样难以实现。以下是几个实用技巧：

抗混叠滤波器的模拟实现

from scipy import signal def apply_antialiasing(x, fs, cutoff_ratio=0.4): """模拟抗混叠滤波器""" nyquist = fs / 2 b, a = signal.butter(4, cutoff_ratio * nyquist / nyquist, 'low') return signal.filtfilt(b, a, x) # 带高频噪声的信号 noisy_signal = np.sin(2*np.pi*2*t_continuous) + 0.3*np.random.randn(len(t_continuous)) # 应用抗混叠滤波器 filtered_signal = apply_antialiasing(noisy_signal, fs=10) plt.figure(figsize=(10,4)) plt.plot(t_continuous, noisy_signal, 'b-', alpha=0.5, label='Noisy Signal') plt.plot(t_continuous, filtered_signal, 'r-', lw=2, label='Filtered Signal') plt.title('Anti-aliasing Filter Demonstration') plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.legend() plt.grid(True) plt.show()

采样保持电路的模拟

def sample_and_hold(t_continuous, t_samples, x_samples): """模拟采样保持电路""" return np.interp(t_continuous, t_samples, x_samples) # 生成采样数据 t_samples_sh = np.arange(0, duration, T) x_samples_sh = continuous_signal(t_samples_sh) # 采样保持 sh_signal = sample_and_hold(t_continuous, t_samples_sh, x_samples_sh) # 绘图 plt.figure(figsize=(10,4)) plt.plot(t_continuous, continuous_signal(t_continuous), 'b-', label='Original') plt.stem(t_samples_sh, x_samples_sh, 'r-', markerfmt='ro', basefmt=" ", label='Samples') plt.plot(t_continuous, sh_signal, 'g-', lw=2, alpha=0.7, label='Sample & Hold') plt.title('Sample-and-Hold Circuit Simulation') plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.legend() plt.grid(True) plt.show()

6. 完整采样系统模拟

将上述组件整合成一个完整的采样系统模拟：

def full_sampling_system(f_signal=1.0, fs=10.0, duration=3.0, noise_level=0.1): """完整的采样系统模拟""" # 1. 生成原始信号（含噪声） t = np.linspace(0, duration, 1000) original = np.sin(2 * np.pi * f_signal * t) + noise_level * np.random.randn(len(t)) # 2. 抗混叠滤波 nyquist = fs / 2 cutoff = 0.4 * nyquist # 保守的截止频率 b, a = signal.butter(4, cutoff / nyquist, 'low') filtered = signal.filtfilt(b, a, original) # 3. 采样 sample_t = np.arange(0, duration, 1/fs) sample_x = np.interp(sample_t, t, filtered) # 4. 采样保持（零阶保持） reconstructed = np.interp(t, sample_t, sample_x) # 计算误差 error = filtered - reconstructed rmse = np.sqrt(np.mean(error**2)) # 绘图 plt.figure(figsize=(12,8)) plt.subplot(3,1,1) plt.plot(t, original, 'b-', label='Original (with noise)') plt.plot(t, filtered, 'r-', lw=2, label='After anti-aliasing filter') plt.title('Input Signal Processing') plt.legend() plt.grid(True) plt.subplot(3,1,2) plt.plot(t, filtered, 'b-', label='Filtered Signal') plt.stem(sample_t, sample_x, 'r-', markerfmt='ro', basefmt=" ", label='Samples') plt.title(f'Sampling Process (fs={fs}Hz, f_signal={f_signal}Hz)') plt.legend() plt.grid(True) plt.subplot(3,1,3) plt.plot(t, filtered, 'b-', label='Original Filtered') plt.plot(t, reconstructed, 'g--', lw=2, label='Reconstructed') plt.title(f'Reconstruction (RMSE: {rmse:.4f})') plt.legend() plt.grid(True) plt.tight_layout() plt.show() # 运行完整系统 full_sampling_system(f_signal=2, fs=10) full_sampling_system(f_signal=4, fs=8) # 边缘情况 full_sampling_system(f_signal=6, fs=10) # 不满足奈奎斯特条件

这个完整系统展示了从原始信号到重建信号的完整流程，包括：