OpenRadar DOA函数 Bartlett/CBF和Capon使用

Bartlett / CBF原理看这里
Capon原理看这里

这里只讲怎么调用openradar提供的aoa_bartlett和aoa_capon函数：
一些吐槽：虽然看起来openradar的作者代码水平很高，但里面有很多匪夷所思的写法，比如他demo里的解析文件格式就很迷啊等等…

1. 波束成形 或者 Angle-FFT 之前要先做 Range-FFT

先对adc_sample维度做距离FFT，一般是最后一个维度，得到相位信息。

2. 获取矢量矩阵 (steering vector):

比如果你想得到[-90°,90°]，分辨率为1°(就有181个角度点)，8根虚拟天线的矢量矩阵，那么就会生成一个shape为[181, 8]的ndarray (也就是下面那行代码的steering_vec)。num_vec得到的是一共有多少角度（上面那个例子的话就是181），没啥用，可以不管。

num_vec, steering_vec = dsp.gen_steering_vec(ANGLE_RANGE, ANGLE_RES, VIRT_ANT)

By the way, dsp.gen_steering_vec里面用了np.complex, 高版本的numpy可能不支持，报错说和python自带的complex有歧义，把np去掉，变成complex即可。

ANGLE_RANGE 是在一侧，你要得到多少角度信息，那么一共就有 2 * ANGLE_RANGE + 1个角度信息，+1是因为0°。
ANGLE_RES 是你需要的角度分辨率
VIRT_ANT 是虚拟天线数
以下举例

ANGLE_RANGE = 90
ANGLE_RES = 1
VIRT_ANT = 8
BINS_PROCESSED = 112

意思想得到的角度为[-90°, 90°]，因为有0，所以有2*ANGLE_RANGE + 1 = 181个角度格子，想要的角度分辨率为1°，处理前112个range-bin。

3.1 Bartlett/CBF:

可以看到源码aoa_bartlett这边已经帮我们对range-angle图取模并且平方放大了。

下面这段是我使用openradar提供的aoa_bartlett API
radar_cube shape : (128, 8, 128) - (chirp, vx, adc_samples)

ra = np.zeros((ANGLE_BINS, BINS_PROCESSED), dtype=np.complex128)
# 只取前BIN_PROCESSED个range-bin处理
ra = dsp.aoa_bartlett(steering_vec, radar_cube[:, :, :BINS_PROCESSED], axis=1) 
ra = ra.sum(0) # 对所有chirp的得到的角度信息求和
plt.title('bartlet sum all chirps')
plt.imshow(ra) # !!! 源码里面帮我们对复数取模并且平方放大了，就不用我们自己做啦。下面那段复现代码是自己取模放大

下面这段是我尝试跟着Bartlett的原理自己写了一个函数，实测效果一致。
只取了第一个chirp，没像上面那样求和，实测这两种方法差不多，对所有chirp得到的角度信息求和（上面那种）的话效果会好一点点点点点。
radar_cube2 shape : (8, 128) - (vx, adc_samples)

ra = np.zeros((ANGLE_BINS, BINS_PROCESSED), dtype=np.complex128)
steering_vec = np.conjugate(steering_vec) # 要取共轭！
for i in range(BINS_PROCESSED): # 遍历每一个range-bin，得到每个range-bin上的角度信息
    tmp = steering_vec @ radar_cube2[:, i]
    ra[:, i] = tmp
plt.title('CBF conjugate')
# openradar bartlett源码就是把得到的每个角度的复数取模然后再平方，这样写是对的，我对比过用dsp.aoa_bartlett的range-angle图，效果一样
plt.imshow(np.abs(ra) ** 2)  # 仿照源码 取模 平方放大信号

3.2 Capon

radar_cube shape : (128, 8, 256) - (chirps, vx, adc_samples)

range_azimuth = np.zeros((ANGLE_BINS, BINS_PROCESSED), dtype=np.complex128)
for i in range(BINS_PROCESSED):
    range_azimuth[:, i], beamWeights[:,i] = dsp.aoa_capon(radar_cube[:, :, i].T, steering_vec, magnitude=False)
plt.title('256 chirps Capon')
plt.imshow(np.abs(range_azimuth))

此外，如果你的chirp小于8，比如1个frame里只有1个chirp，那么就要把openradar的dsp.angle_estimation.cov_matrix里的这个判断删了：

if x.shape[0] > x.shape[1]:
   warnings.warn("cov_matrix input should have Vrx as rows. Needs to be transposed", RuntimeWarning)
   x = x.T

大概意思就是他觉得你输错了，还帮你改正了，但实际上你没输错，他帮了倒忙，就是要处理一帧里面，1个chirp，8根虚拟天线的数据，anyway，如果一帧里只有小于虚拟天线个数的chirp数的话，把这段判断删了或注释掉就行。

下面的视频是我对这Bartlett/CBF和Capon和Angle-FFT在openradar的circle.bin这个数据上的效果对比：
视频链接