4. 频域波束形成

任务说明

本案例实现频域窄带波束形成：将宽带信号经FFT分解为多个频点，对每个频点独立计算阵列流形向量并施加复加权，完成空域滤波后合成输出。最终比较不同频率下的波束方向图，观察主瓣宽度随频率的变化。

频域波束形成原理

对于均匀线阵，第k个频点f_k对应的阵列流形向量为：

a(f_k, θ) = [1, e^{−j2πf_kd·sinθ/c}, …, e^{−j2πf_k(M−1)d·sinθ/c}]^T

常规波束形成（CBF）权向量取w = a(f_k, θ₀)，其中θ₀为期望观测方向。对频点f_k的阵列接收快拍x(f_k)（M×1向量），输出为y(f_k) = w^Hx(f_k)。方向图定义为P(θ) = |w^Ha(f_k, θ)|²，扫描θ得到空间响应。

代码实现

以下Python代码实现频域波束形成核心功能：生成多频点窄带信号、逐频点计算CBF权向量、扫描方向图并绘图比较。

import numpy as np
import matplotlib.pyplot as plt

# 参数配置
M = 8                # 阵元数
d = 0.5              # 阵元间距（以最高频率波长归一化）
c = 1500.0           # 声速 m/s
theta0 = 30.0        # 期望方向（度）
fs = 10000           # 采样率 Hz
N = 1024             # FFT点数

# 多频点设置（窄带假设：各频点独立处理）
freqs = [1000, 2000, 3000]   # 三个频点 Hz
angles_deg = np.linspace(-90, 90, 361)  # 扫描角度
theta0_rad = np.deg2rad(theta0)

# 生成模拟接收数据（每个频点一个复正弦快拍）
# 信号从theta0方向入射，含噪声
np.random.seed(0)
t = np.arange(N) / fs
x_freq = {}   # 存储各频点的频域快拍（M×1复数向量）

for fk in freqs:
    # 时域阵列信号：M个阵元，每个阵元含时延
    X_time = np.zeros((M, N), dtype=complex)
    for m in range(M):
        tau = m * d * np.sin(theta0_rad) / c   # 相对时延
        X_time[m, :] = np.exp(1j * 2 * np.pi * fk * (t - tau))
    # 加噪声
    X_time += 0.1 * (np.random.randn(M, N) + 1j * np.random.randn(M, N))
    # FFT取频点fk处的值
    X_fft = np.fft.fft(X_time, axis=1)
    k_idx = int(fk / fs * N)   # 频点索引
    x_freq[fk] = X_fft[:, k_idx]   # M×1 复数快拍

# 频域波束形成：逐频点计算方向图
def array_manifold(fk, theta_deg):
    """生成频率fk、角度theta_deg下的阵列流形向量"""
    theta_rad = np.deg2rad(theta_deg)
    # 阵元位置（以波长归一化，d固定为0.5*lambda_max，但频率变化时电长度变化）
    # 实际阵元间距物理值 = d * c / max(freqs)，此处用物理间距计算
    d_phys = d * c / max(freqs)   # 物理间距（米）
    positions = np.arange(M) * d_phys
    # 波数 k = 2πf/c
    k_wave = 2 * np.pi * fk / c
    return np.exp(-1j * k_wave * positions * np.sin(theta_rad))

# 计算各频点方向图
patterns = {}
for fk in freqs:
    a_theta0 = array_manifold(fk, theta0)   # 期望方向导向矢量
    w = a_theta0 / M   # CBF权（归一化）
    # 扫描方向图
    P = np.zeros(len(angles_deg))
    for i, ang in enumerate(angles_deg):
        a_scan = array_manifold(fk, ang)
        P[i] = np.abs(np.dot(w.conj(), a_scan))**2
    patterns[fk] = P

# 绘制方向图比较
plt.figure(figsize=(10, 6))
for fk in freqs:
    plt.plot(angles_deg, 10*np.log10(patterns[fk] + 1e-12), 
             label=f'f={fk} Hz')
plt.xlabel('角度 (deg)')
plt.ylabel('归一化功率 (dB)')
plt.title('频域波束形成方向图比较')
plt.legend()
plt.grid(True)
plt.ylim(-60, 5)
plt.xlim(-90, 90)
plt.axvline(theta0, color='gray', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.savefig('frequency_domain_beam_patterns.png', dpi=150)
plt.show()

# 输出主瓣宽度（-3dB）对比
print("各频点主瓣宽度（-3dB）:")
for fk in freqs:
    P_linear = patterns[fk]
    P_db = 10*np.log10(P_linear + 1e-12)
    peak_idx = np.argmax(P_db)
    # 找-3dB交叉点
    half_power = P_db[peak_idx] - 3
    left_idx = peak_idx
    while left_idx > 0 and P_db[left_idx] > half_power:
        left_idx -= 1
    right_idx = peak_idx
    while right_idx < len(P_db)-1 and P_db[right_idx] > half_power:
        right_idx += 1
    bw = angles_deg[right_idx] - angles_deg[left_idx]
    print(f"  f={fk} Hz: 主瓣宽度 ≈ {bw:.1f}°")

关键实现说明

代码中array_manifold函数根据物理阵元间距和频率计算导向矢量，正确反映了频率对阵列响应的影响。CBF权向量取期望方向的导向矢量共轭并除以M做幅度归一化。方向图扫描在−90°到90°范围内进行，功率以dB显示。

运行与验证

运行方式：保存代码为freq_beamforming.py，在装有NumPy和Matplotlib的Python环境中执行python freq_beamforming.py。

预期输出：程序生成方向图比较曲线，并在控制台打印各频点主瓣宽度。根据阵列理论，频率越高（波长越短），相同物理孔径下主瓣越窄。因此预期输出中3000 Hz的主瓣宽度应小于1000 Hz的主瓣宽度。

判定标准：观察生成的frequency_domain_beam_patterns.png，三个频点的方向图主瓣均指向30°，且高频方向图主瓣明显更窄、旁瓣更多；控制台输出的主瓣宽度数值随频率单调递减，即满足f=1000 Hz时宽度最大、f=3000 Hz时宽度最小。

验证要点

运行代码后检查：①方向图主瓣是否对准30°；②三个频点的主瓣宽度是否呈现随频率升高而变窄的趋势；③控制台输出的主瓣宽度数值是否满足单调递减关系。若频率设置改为更高或更低，主瓣宽度变化趋势应保持一致，这验证了频域波束形成对不同频率信号空间分辨率的物理特性。