常规波束形成——时域波束形成、频域波束形成

最近在写代码的时候发现自己对波束成形的理解依旧不是很透彻，就重新找资料学习了一下，整理一下加深印象并在此做个记录。参考资料如下：
https://zhuanlan.zhihu.com/p/10781326364
https://blog.csdn.net/Xujing1143/article/details/120182410
https://www.sharetechnote.com/html/Handbook_LTE_BeamForming.html

PS：如果大致知道波束形成是什么，但是对它的时延计算、复指数$e^{-jwt}$的符号由来、相位补偿、又或者是频域波束形成有点迷的小伙伴可以从实例开始看~

1. 基本概念

波束形成是一种利用多个传感器（如麦克风或天线）的阵列处理技术。它是指将一定几何形状排列的多元基阵的各个阵元输出，经过处理（如加权、时延、求和等）形成空间指向性的方法。

• 目的：提取特定方向来的信号，抑制其他方向上的信号。因此波束形成器可以等效为一个空间滤波器。
• 实际应用：估计目标信号的存在性，并对信号的入射方位进行估计。
• 核心思想：
  ○ 信号到达不同传感器的时间不同（取决于信号的入射角）；
  ○ 利用这一定向性信息，通过加权和调整延迟来增强某个方向的信号，同时抑制其他方向的信号。
• 分类：一般分为常规波束形成（CBF）和自适应波束形成（ABF）。

2. 常规波束形成

根据实现的途径不同，分为时延波束形成和频移波束形成。

2.1 时延波束形成（Delay-and-Sum Beamforming）

1. 概念

• 直接对接收信号进行时延补偿（对准某个方向的波前），然后对所有传感器的数据求和。
• 在某个方向上完全对齐的信号会相干增强，而其他方向的信号由于相移不同，会部分或完全相互抵消。

2. 公式

设传感器接收到的信号为$s_i(t)$，导向向量为$\mathbf{a}(\theta )$，则波束形成输出为：
$$y(t)={\Sigma }_{i=1}^M{w}_i{s}_i(t)$$
其中：

• $M$是传感器个数。
• $w_i$是波束形成的权重（在常规波束形成中，通常$w_i=1$）。
• 若信号的方向为${\theta }_d$，则需要对$s_i(t)$进行适当的延迟补偿，使得所有信号在该方向上对齐。

3. 缺点

• 不能很好地抑制干扰和噪声，如果有多个信号源，干扰较大。
• 角度分辨率较低。

4. 实例：以均匀线性阵（ULA）为例
   
   （1）阵列接收信号的表示
   假设远场信号为频率为 $f_0$ 的单频信号：$x(t)=\cos (2\pi f_{0}t)$
   相邻阵元间存在时延，对于任意角度 $\theta$，相邻阵元的时延差为：
   $$\tau =\frac{d\sin (\theta )}{c}$$
   其中，$d$ 为阵元间距，$\theta$ 为信号来源方向与基阵法线的夹角。以第一个阵元为参考阵元，则每个阵元相对其的时间差可以表示为：
   $$\Delta t=\left[0,\tau ,2\tau ,\cdots ,(M-1)\tau ]=[0,\frac{d\sin (\theta )}{c},\frac{2d\sin (\theta )}{c},\cdots ,\frac{(M-1)d\sin (\theta )}{c}\right]$$
   相应地，相位差为：
   $$\Delta \varphi =\left[0,2\pi f_0\frac{d\sin (\theta )}{c},2\pi f_0\frac{2d\sin (\theta )}{c},\cdots ,2\pi f_0\frac{(M-1)d\sin (\theta )}{c}\right]$$
   由于这个相位差是阵元间空间位置不同所造成的，因此称之为“空间相位差”。

在某一时刻对所有每个阵元同时采样（一个快拍），将会得到 $M$ 个数据点，将这 $M$ 个数据点排列起来，可以得到这一时刻阵列的接收信号：
$$\begin{array}{rl}\mathbf{y}(n)& =\left[x(n),x(n)e^{-j2\pi f_0\frac{d\sin \left({\theta }_d\right)}{c}},x(n)e^{-j2\pi f_0\frac{2d\sin \left({\theta }_d\right)}{c}},\cdots ,x(n)e^{-j2\pi f_0\frac{(M-1)d\sin \left({\theta }_d\right)}{c}}\right]\\ & =\underset{\mathbf{a}\left({\theta }_d\right)}{\underbrace{\left[1,e^{-j2\pi f_0\frac{d\sin \left({\theta }_d\right)}{c}},e^{-j2\pi f_0\frac{2d\sin \left({\theta }_d\right)}{c}},\cdots ,e^{-j2\pi f_0\frac{(M-1)d\sin \left({\theta }_d\right)}{c}}\right]}}x(n)\\ & =\mathbf{a}\left({\theta }_d\right)\cdot x(n)\end{array}$$
可见，每个快拍阵列接收到的信号 $\mathbf{y}(n)$ 是一个向量 $\mathbf{a}({\theta }_d)$ 乘以标量 $\mathbf{x}(n)$，且 $\mathbf{y}(n)$ 是来波方向 ${\theta }_d$ 的函数。

向量 $\mathbf{a}({\theta }_d)$ 包含了不同阵元的空间相位差，这便是常说的导向矢量。

Question：阵元接收信号的复指数为什么是负的？比如是 $x(n)e^{-\mathbf{j}2\mathbf{\pi }{\mathbf{f}}_0\frac{\mathbf{d}\sin \left({\mathbf{\theta }}_{\mathbf{d}}\right)}{\mathbf{c}}}$ 而不是 $x(n)e^{\mathbf{j}2\mathbf{\pi }{\mathbf{f}}_0\frac{\mathbf{d}\sin \left({\mathbf{\theta }}_{\mathbf{d}}\right)}{\mathbf{c}}}$？

这和我们作出的来波示意图有关，下图清晰地解释了该表达式是怎么来的。

（来自单一源的波束到不同天线的路径差表示为 $e^{-jdn\theta }$，此处 $\sin \theta \approx \theta if\theta \to 0$）

（2）波束形成
阵列在某一时刻的接收信号为：
$$\begin{array}{rl}\mathbf{y}(n)& =\left[1,e^{-j2\pi f_0\frac{d\sin \left({\theta }_d\right)}{c}},e^{-j2\pi f_0\frac{2d\sin \left({\theta }_d\right)}{c}},\cdots ,e^{-j2\pi f_0\frac{(M-1)d\sin \left({\theta }_d\right)}{c}}\right]x(n)\\ & =\mathbf{a}\left({\theta }_d\right)\cdot x(n)\end{array}$$

• 入射角估计
  虽然入射角 ${\theta }_d$ 未知，但给定的阵列结构的导向矢量的数学形式是已知的。
  我们可以构造一个方向为 $\alpha$ 的导向矢量，其中 $\alpha$ 是我们猜测的信号的入射角度：
  $$\mathbf{a}(\alpha )=\left[1,e^{-j2\pi f_0\frac{d\sin \left(\alpha \right)}{c}},e^{-j2\pi f_0\frac{2d\sin \left(\alpha \right)}{c}},\cdots ,e^{-j2\pi f_0\frac{(M-1)d\sin \left(\alpha \right)}{c}}\right]$$
  用这个导向矢量与阵列接收信号 $\mathbf{y}(n)$ 做内积：
  $$\begin{array}{rl}& \mathbf{a}(\alpha )\cdot \mathbf{y}(n)\\ & =\mathbf{a}(\alpha {)}^H\cdot \mathbf{a}\left({\theta }_d\right)\cdot x(n)\\ & =\left[1+e^{-j2\pi f_{0}d\frac{|\sin (\alpha )-\sin ({\theta }_d)|}{c}}+e^{-j2\pi f_{0}d\frac{2|\sin (\alpha )-\sin ({\theta }_d)|}{c}}+\cdots +e^{-j2\pi f_{0}d\frac{(M-1)|\sin (\alpha )-\sin ({\theta }_d)|}{c}}\right]x(n)\\ & \le M\cdot x(n)\end{array}$$
  当 $\alpha ={\theta }_d$ 时，取等号。

此处手动强调：使用的是 $\mathbf{a}(\alpha {)}^H$，$\mathbf{a}(\alpha )$ 的共轭转置！！！后文频域波束形成进行的相位补偿就是源于此

该不等式表明，若猜测的方向和信号真实的入射方向一致，此时得到的内积值最大。所以，我们可以通过遍历角度来估计信号的入射方向，最大值对应的角度即为预测的入射方向。可以表述为如下形式：
$$\begin{array}{rl}& for\alpha from-\pi to\pi \\ & S(\alpha )=\mathbf{a}(\alpha {)}^H\cdot \mathbf{y}(n)\\ & {\theta }_d=\underset{\alpha }{\mathrm{argmax}}S(\alpha )\end{array}$$

• 波束响应
  要使目标方向形成最大响应，我们可以对各阵元补偿指定方向 ${\theta }_d$ 的相位，构建最优加权向量：
  $$w=\mathbf{a}({\theta }_d)=\left[1,e^{-j2\pi f_0\frac{d\sin \left({\theta }_d\right)}{c}},e^{-j2\pi f_0\frac{2d\sin \left({\theta }_d\right)}{c}},\cdots ,e^{-j2\pi f_0\frac{(M-1)d\sin \left({\theta }_d\right)}{c}}\right]$$
  则输出信号为：
  $$Y(n)=w\cdot \mathbf{y}(n)$$
  输出的波束响应：
  $$B(\alpha )=w={\Sigma }_{i=1}^M{w}_i\cdot \mathbf{a}(\alpha )$$
  $B(\alpha )$ 即为各方向的响应。

2.2 频域波束形成

频域波束形成类似，主要步骤如下：
（1）对阵列接收信号 $\mathbf{y}(n)$ 进行快速傅里叶变换，将其从时域转换为频域：
$$X=FFT(\mathbf{y}(n))$$
其中：

• $X$为$M\times N$的二维向量；
• $M$是阵元个数；
• $N$是频域点数。

（2）对各阵元信号进行指定方向 $\alpha$ 的相位补偿：
$$\begin{array}{rl}& forfinfreqs:\\ & forkfrom0to(M-1):\\ & phase\_shift\left[k,f\right]=e^{-j2\pi f\frac{kd\sin \left(\alpha \right)}{c}}\end{array}$$

补偿后的信号：
$$\begin{array}{rl}phase\_shift& =phase\_shift.\mathbf{conj}()\\ X^{\prime }& =phase\_shift\cdot X\end{array}$$

（3）对各个阵元信号的某一个频率进行求和，再对各阵元信号的所有频率对应的信号值进行求和，可以得到该角度的响应幅度。

（4）对给定角度范围进行遍历，可以得到不同角度下的响应。

• 螺旋阵单声源的运行效果如下：
  
  PS：！！！！频域补偿的符号一定要注意啊啊啊啊。调试n天的血泪教训 … 做出来的结果终于对了！

• 实现代码如下：

def get_steer_vector(positions, theta_rad, freq, c=340):
    u = np.array([np.cos(theta_rad), np.sin(theta_rad)])
    tau = positions.T @ u / c  # 正确时延计算
    return np.exp(-1j * 2 * np.pi * freq * tau)

def freq_beamforming(waveform, fs, theta_range, array_geometry, c=340, phi=0):
    num_sensors = waveform.shape[0]
    
    # 对信号进行FFT 计算频谱
    n_fft = waveform.shape[1]
    fft_values = np.fft.fft(waveform, n_fft, axis=1)  # (num_sensors, freqs)
    fft_freqs = np.fft.fftfreq(n_fft, 1 / fs) 
    # 只取前半部分频率
    half_n = len(fft_freqs) // 2

    power = []
    for theta in theta_range:
        svs = []
        for i in range(half_n):
            f = fft_freqs[i]
            sv = get_steer_vector(array_geometry, theta, f, c) # 计算相位补偿
            svs.append(sv)
        svs = np.array(svs).reshape(-1, num_sensors)
        phase_shift = svs.conj().T          # 注意相位补偿的符号
        sum_signal = np.sum(fft_values[:, :half_n] * phase_shift, axis=0) 
        power.append(np.sum(np.abs(sum_signal) ** 2))
    # 归一化
    power /= np.max(power)

    plt.figure(figsize=(8, 5))
    plt.plot(np.rad2deg(theta_range), power)
    plt.xlabel("Angle (degrees)")
    plt.ylabel("Normalized Energy")
    plt.legend()
    plt.title("Frequency Domain Beamforming")
    plt.show()
    return power

wav_path = './record/2kHz_sine_25degree.wav'  # 根据自己的待测音频路径进行修改
waveform1, fs = torchaudio.load(wav_path)

# 读取阵元位置文件
arrangement_path = './DoA/coordinate.csv'
df = pd.read_csv(arrangement_path)
arrangement = df[['X', 'Y']].to_numpy().T   # 单位为 m
# arrangement 根据自己数据采集阵列的位置参数进行修改，也可以构建 ULA 来模拟

theta_scan = np.linspace(- np.pi, np.pi, 360, endpoint=False)
power = freq_beamforming(signal, fs, theta_scan, arrangement)