📚 课程目录第1章基于互相关的布里渊频移提取

1. 基于互相关的布里渊频移提取

任务说明

本案例实现从布里渊散射谱中提取布里渊频移（BFS）的互相关算法。输入为一段离散的布里渊增益谱（BGS）数据，包含已知的中心频率（即真实BFS）与噪声；输出为算法估计的BFS值，并通过与真实值的偏差评估精度，通过不同信噪比下的估计结果评估抗噪性。

算法原理与实现

互相关法提取BFS的核心思想：将实测BGS与一组已知中心频率的参考洛伦兹谱模板逐一计算互相关，互相关值最大时对应的模板中心频率即为BFS估计值。相比直接寻峰，互相关对噪声和谱形畸变具有更好的鲁棒性。

代码实现包含三个核心函数：lorentzian()生成洛伦兹谱模板；cross_correlate()计算两序列的互相关值；extract_bfs()在预设频率搜索范围内扫描模板，返回使互相关最大的频率。主程序生成含噪BGS，调用提取函数并输出估计结果。

import numpy as np

def lorentzian(f, f0, width=30e6, peak=1.0):
    """洛伦兹线型谱模板
    f: 频率数组 (Hz)
    f0: 中心频率 (Hz)
    width: 半高全宽 (Hz)
    peak: 峰值幅度
    """
    return peak * (width/2)**2 / ((f - f0)**2 + (width/2)**2)

def cross_correlate(x, y):
    """计算两个等长序列的互相关值（归一化）"""
    x = x - np.mean(x)
    y = y - np.mean(y)
    denom = np.sqrt(np.sum(x**2) * np.sum(y**2))
    if denom == 0:
        return 0.0
    return np.sum(x * y) / denom

def extract_bfs(freq, bgs, search_range, step=1e6):
    """通过互相关扫描提取BFS
    freq: 频率数组 (Hz)
    bgs: 实测布里渊增益谱
    search_range: (fmin, fmax) 搜索范围 (Hz)
    step: 扫描步长 (Hz)
    """
    fmin, fmax = search_range
    scan_freqs = np.arange(fmin, fmax + step, step)
    best_corr = -1.0
    best_f0 = None
    for f0 in scan_freqs:
        template = lorentzian(freq, f0)
        corr = cross_correlate(bgs, template)
        if corr > best_corr:
            best_corr = corr
            best_f0 = f0
    return best_f0

# ===== 参数设置 =====
freq = np.linspace(10.0e9, 11.0e9, 2001)      # 频率扫描范围 10~11 GHz
true_bfs = 10.580e9                            # 真实布里渊频移
snr_db = 15                                    # 信噪比 (dB)

# 生成无噪BGS（洛伦兹线型）
clean_bgs = lorentzian(freq, true_bfs)

# 添加高斯白噪声
signal_power = np.mean(clean_bgs**2)
noise_power = signal_power / (10**(snr_db/10))
noise = np.random.normal(0, np.sqrt(noise_power), size=freq.shape)
noisy_bgs = clean_bgs + noise

# 互相关提取BFS
search_range = (10.4e9, 10.7e9)               # 搜索范围包含真实值
estimated_bfs = extract_bfs(freq, noisy_bgs, search_range, step=0.5e6)

# 输出结果
error_mhz = (estimated_bfs - true_bfs) / 1e6
print(f"真实BFS: {true_bfs/1e9:.4f} GHz")
print(f"估计BFS: {estimated_bfs/1e9:.4f} GHz")
print(f"误差: {abs(error_mhz):.2f} MHz")

抗噪性评估

为定量评估算法在不同噪声水平下的表现，以下代码在多个SNR下重复估计，统计误差变化趋势：

# ===== 抗噪性测试 =====
snr_list = [5, 10, 15, 20, 25, 30]
errors = []
for snr in snr_list:
    noise_power = signal_power / (10**(snr/10))
    trial_errors = []
    for _ in range(20):  # 每个SNR重复20次
        noise = np.random.normal(0, np.sqrt(noise_power), size=freq.shape)
        noisy = clean_bgs + noise
        est = extract_bfs(freq, noisy, search_range, step=0.5e6)
        trial_errors.append(abs(est - true_bfs)/1e6)
    errors.append(np.mean(trial_errors))

print("\nSNR(dB)   平均误差(MHz)")
for snr, err in zip(snr_list, errors):
    print(f"{snr:6d}    {err:10.3f}")

落地与验证

运行方式：将两段代码合并保存为 bfs_cross_correlation.py，在装有NumPy的Python 3环境中执行 python bfs_cross_correlation.py。

验证方法：观察第一段输出中“估计BFS”与“真实BFS”的偏差，在SNR=15dB时误差应小于5MHz；第二段输出应显示SNR从5dB增至30dB时，平均误差总体呈下降趋势，且SNR≥15dB时误差通常小于3MHz。

判定标准：程序无异常退出；估计值落在搜索范围内且与真实值偏差在合理范围；抗噪性测试中误差随SNR提高而减小，即证明互相关算法能有效提取BFS并具备抗噪能力。