别再死记硬背了！用Python+NumPy手把手复现N-P定理，理解信号检测的本质

用Python实战N-P定理：从数学公式到信号检测的代码实现

第一次接触N-P定理时，我被那一堆积分符号和概率密度函数搞得晕头转向。直到有一天，我决定用Python把整个过程模拟出来，那些抽象的数学概念突然变得清晰可见。这就是我们今天要做的——通过代码让N-P定理"活"起来。

1. 准备工作：理解问题场景

假设你正在设计一个雷达系统，需要从噪声中检测是否有目标存在。这就是典型的二元假设检验问题：

• H₀假设（无信号）：仅观测到噪声
• H₁假设（有信号）：观测到信号加噪声

N-P定理告诉我们：在给定虚警概率(P_FA)的条件下，如何最大化检测概率(P_D)。听起来很抽象？让我们用代码来具象化这个过程。

首先安装必要的库：

pip install numpy matplotlib scipy

然后导入我们将用到的模块：

import numpy as np import matplotlib.pyplot as plt from scipy.stats import norm

2. 生成模拟数据

我们假设信号和噪声都服从高斯分布，这是信号处理中最常见的假设：

# 参数设置 n_samples = 10000 # 样本数量 signal_power = 1 # 信号功率 noise_power = 1 # 噪声功率 snr = signal_power / noise_power # 信噪比 # 生成H0假设下的数据（仅噪声） h0_samples = np.random.normal(0, np.sqrt(noise_power), n_samples) # 生成H1假设下的数据（信号+噪声） h1_samples = np.random.normal(1, np.sqrt(noise_power), n_samples)

可视化这两类数据的分布：

plt.figure(figsize=(10, 6)) plt.hist(h0_samples, bins=50, alpha=0.5, label='H0 (Noise only)') plt.hist(h1_samples, bins=50, alpha=0.5, label='H1 (Signal+Noise)') plt.xlabel('Amplitude') plt.ylabel('Count') plt.legend() plt.title('Distribution of H0 and H1 samples') plt.show()

3. 实现似然比检测

N-P定理的核心是似然比检验(LRT)。对于高斯分布的情况，似然比可以简化为：

def likelihood_ratio(x, signal_mean, noise_variance): """ 计算似然比 L(x) = p(x|H1)/p(x|H0) 对于高斯分布，可以简化为这个形式 """ return np.exp((x * signal_mean - 0.5 * signal_mean**2) / noise_variance)

接下来我们需要确定判决门限γ。根据N-P定理，γ由给定的虚警概率P_FA决定：

def find_threshold(pfa, h0_samples): """ 根据给定的P_FA计算判决门限 """ # 对H0样本计算似然比 lrt_h0 = likelihood_ratio(h0_samples, 1, 1) # 找到使P_FA等于指定值的门限 sorted_lrt = np.sort(lrt_h0) index = int((1 - pfa) * len(sorted_lrt)) return sorted_lrt[index]

4. 性能评估与ROC曲线

现在我们可以计算不同P_FA下的检测概率P_D：

def calculate_pd(pfa, h0_samples, h1_samples): """ 计算给定P_FA下的P_D """ threshold = find_threshold(pfa, h0_samples) lrt_h1 = likelihood_ratio(h1_samples, 1, 1) pd = np.mean(lrt_h1 > threshold) return pd # 计算ROC曲线 pfa_values = np.linspace(0, 1, 100) pd_values = [calculate_pd(pfa, h0_samples, h1_samples) for pfa in pfa_values] # 绘制ROC曲线 plt.figure(figsize=(8, 6)) plt.plot(pfa_values, pd_values) plt.plot([0, 1], [0, 1], 'k--') plt.xlabel('Probability of False Alarm (P_FA)') plt.ylabel('Probability of Detection (P_D)') plt.title('ROC Curve for N-P Detector') plt.grid(True) plt.show()

5. 实际应用中的考量

在实际工程中，我们还需要考虑几个关键因素：

参数估计误差的影响

# 假设我们估计的信号幅度有误差 estimated_signal = 0.9 # 实际为1.0 # 使用错误参数计算似然比 def mismatched_likelihood_ratio(x, estimated_signal, noise_variance): return np.exp((x * estimated_signal - 0.5 * estimated_signal**2) / noise_variance) # 计算性能下降情况 threshold_mismatch = find_threshold(0.1, h0_samples) lrt_h1_mismatch = mismatched_likelihood_ratio(h1_samples, estimated_signal, 1) pd_mismatch = np.mean(lrt_h1_mismatch > threshold_mismatch) print(f"P_D with parameter mismatch: {pd_mismatch:.3f}")

不同信噪比下的性能比较

snr_values = [0.5, 1, 2, 5] roc_curves = {} for snr in snr_values: # 生成不同SNR的数据 h1_samples_snr = np.random.normal(snr, 1, n_samples) pd_values_snr = [calculate_pd(pfa, h0_samples, h1_samples_snr) for pfa in pfa_values] roc_curves[f"SNR={snr}dB"] = pd_values_snr # 绘制不同SNR的ROC曲线 plt.figure(figsize=(10, 6)) for label, pd_vals in roc_curves.items(): plt.plot(pfa_values, pd_vals, label=label) plt.xlabel('Probability of False Alarm (P_FA)') plt.ylabel('Probability of Detection (P_D)') plt.title('ROC Curves for Different SNR Values') plt.legend() plt.grid(True) plt.show()

6. 扩展到多维情况

现实中的信号往往是多维的。让我们考虑二维信号的检测问题：

# 生成二维数据 n_samples_2d = 5000 h0_2d = np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], n_samples_2d) h1_2d = np.random.multivariate_normal([1, 1], [[1, 0], [0, 1]], n_samples_2d) # 二维似然比函数 def likelihood_ratio_2d(x, signal_mean, noise_cov): inv_noise_cov = np.linalg.inv(noise_cov) exponent_h0 = -0.5 * x @ inv_noise_cov @ x.T exponent_h1 = -0.5 * (x - signal_mean) @ inv_noise_cov @ (x - signal_mean).T return np.exp(exponent_h1 - exponent_h0) # 计算二维情况下的ROC曲线 lrt_h0_2d = np.array([likelihood_ratio_2d(x, [1,1], [[1,0],[0,1]]) for x in h0_2d]) lrt_h1_2d = np.array([likelihood_ratio_2d(x, [1,1], [[1,0],[0,1]]) for x in h1_2d]) pfa_values_2d = np.linspace(0, 1, 50) pd_values_2d = [] for pfa in pfa_values_2d: threshold = np.percentile(lrt_h0_2d, 100*(1-pfa)) pd = np.mean(lrt_h1_2d > threshold) pd_values_2d.append(pd) # 绘制二维ROC曲线 plt.figure(figsize=(8, 6)) plt.plot(pfa_values_2d, pd_values_2d) plt.plot([0, 1], [0, 1], 'k--') plt.xlabel('P_FA (2D case)') plt.ylabel('P_D (2D case)') plt.title('ROC Curve for 2D Signal Detection') plt.grid(True) plt.show()

7. 工程实践中的优化技巧

在实际系统中，我们还需要考虑计算效率和实时性。以下是几个实用技巧：

对数似然比简化计算

def log_likelihood_ratio(x, signal_mean, noise_variance): """ 使用对数似然比避免数值下溢 """ return (x * signal_mean - 0.5 * signal_mean**2) / noise_variance # 比较两种计算方式的数值稳定性 x_large = 10 print(f"Original LRT: {likelihood_ratio(x_large, 1, 0.1)}") print(f"Log-LRT: {log_likelihood_ratio(x_large, 1, 0.1)}")

滑动窗口检测实现

def sliding_window_detection(signal, window_size, threshold): """ 实现滑动窗口检测 """ n = len(signal) decisions = np.zeros(n) for i in range(n - window_size + 1): window = signal[i:i+window_size] llr = np.mean(log_likelihood_ratio(window, 1, 1)) if llr > threshold: decisions[i:i+window_size] = 1 return decisions # 生成含脉冲信号的测试数据 test_signal = np.random.normal(0, 1, 1000) test_signal[300:350] += 1 # 加入信号脉冲 # 执行检测 detection_result = sliding_window_detection(test_signal, 20, 0.5) # 可视化结果 plt.figure(figsize=(12, 6)) plt.plot(test_signal, label='Signal') plt.plot(detection_result, label='Detection Result', alpha=0.7) plt.legend() plt.title('Sliding Window Detection Example') plt.show()

通过这次代码实践，N-P定理从一个抽象的数学概念变成了可以直观感受的工程工具。在雷达系统中，我们可能会用C++实现更高效的版本；在医疗信号处理中，可能需要考虑更复杂的噪声模型。但核心思想不变：在可控的虚警概率下，最大化我们的检测能力。