radarsimpy 是一个用于雷达仿真的 Python 库项目。功能齐全，文档完善，很好用。
本文将使用 radarsimpy 对TDM MIMO FCMW 雷达进行仿真，并实践 Range-Doppler 和 Angle-FFT。
创建雷达系统
创建一个 2T4R 的雷达系统。
- 中心频率 \(f_c\) 为 60.5e
- 单 chirp 时长 \(Tc\) 为 16e-6
- 整个 chirp 周期时长为 40e-6
- 频率调制范围在 [61e9, 60e9]
- 一帧数量为 512
- 采样率 \(f_s\) 为 20e6
关于雷达系统的位置和脉冲信息，可以看后文的图示。
import numpy as np
from radarsimpy import Radar, Transmitter, Receiver
from scipy.constants import speed_of_light
wavelength = speed_of_light / 60.5e9
N_tx = 2
N_rx = 4
tx = []
tx.append(
dict(
location=(0, -N_rx / 2 * wavelength, 0),
delay=0,
)
)
tx.append(
dict(
location=(0, 0, 0),
delay=20e-6,
)
)
rx = [dict(location=(0, wavelength / 2 * idx, 0)) for idx in range(N_rx)]
tx = Transmitter(f=[61e9, 60e9], t=16e-6, tx_power=15, prp=40e-6, pulses=512, channels=tx)
rx = Receiver(fs=20e6, noise_figure=8, rf_gain=20, load_resistor=500, baseband_gain=30, channels=rx)
radar = Radar(transmitter=tx, receiver=rx)
显示阵列的位置信息。
import plotly.graph_objs as go
from IPython.display import Image
fig = go.Figure()
fig.add_trace(
go.Scatter(
x=radar.radar_prop["transmitter"].txchannel_prop["locations"][:, 1] / wavelength,
y=radar.radar_prop["transmitter"].txchannel_prop["locations"][:, 2] / wavelength,
mode='markers',
name='Transmitter',
opacity=0.7,
marker=dict(size=10),
)
)
fig.add_trace(
go.Scatter(
x=radar.radar_prop["receiver"].rxchannel_prop["locations"][:, 1] / wavelength,
y=radar.radar_prop["receiver"].rxchannel_prop["locations"][:, 2] / wavelength,
mode='markers',
opacity=1,
name='Receiver',
)
)
fig.update_layout(
title='Array configuration',
height=400,
width=700,
xaxis=dict(title='y (λ)'),
yaxis=dict(title='z (λ)', scaleanchor="x", scaleratio=1),
)
fig.show()
创建目标，获得回波数据
创建一个目标，位于 30m 处、10° 的位置，其相对雷达的速度为 25m/s。
true_theta = [10]
speed = [25]
ranges = [30]
targets = []
for theta, v, r in zip(true_theta, speed, ranges):
target = dict(
location=(r * np.cos(np.radians(theta)), r * np.sin(np.radians(theta)), 0),
speed=(v * np.cos(np.radians(theta)), v * np.sin(np.radians(theta)), 0),
rcs=10,
phase=0,
)
targets.append(target)
计算回波数据。
from radarsimpy.simulator import simc
data = simc(radar, targets, noise=False)
timestamp = data['timestamp']
baseband = data['baseband']
现在让我们看下两个 Tx 天线发射的信号：
# reshaped_timestamp 维度：Tx，Rx，脉冲，采样点
reshaped_timestamp = timestamp.reshape((N_tx, N_rx, 512, -1))
fig = go.Figure()
bandwidth = radar.radar_prop["transmitter"].waveform_prop["bandwidth"]
center_freq = 60.5e9
start_freq = center_freq - bandwidth / 2
stop_freq = center_freq + bandwidth / 2
samples_per_pulse = radar.sample_prop["samples_per_pulse"]
for Tx in range(2):
x_data = reshaped_timestamp[Tx, 0, 0, :] * 1e6 # x 轴。单位为微秒
y_data = np.linspace(start_freq, stop_freq, samples_per_pulse, endpoint=False) / 1e9
trace_name = f'Tx {Tx}, Chirp 0'
fig.add_trace(go.Scatter(x=x_data, y=y_data, name=trace_name))
fig.update_layout(
title='Transmitter chirps',
height=400,
width=700,
yaxis=dict(tickformat='.2f', title='Frequency (GHz)'),
xaxis=dict(tickformat='.2f', title='Time (µs)'),
)
fig.show()
进行 Range-Doppler
首先定义函数：
def range_fft(X, n=None, window=None):
if window is not None:
shape = [1] * X.ndim
shape[-1] = X.shape[-1]
X = X * window.reshape(shape)
output = np.fft.fft(X, n, axis=-1)
return output
def doppler_fft(X, n=None, window=None):
if window is not None:
shape = [1] * X.ndim
shape[-2] = X.shape[-2]
X = X * window.reshape(shape)
output = np.fft.fftshift(np.fft.fft(X, n, axis=-2), axes=-2)
return output
然后进行 Range-Doppler 处理：
from scipy import signal
range_window = signal.windows.chebwin(radar.sample_prop["samples_per_pulse"], at=80)
dop_window = signal.windows.chebwin(radar.radar_prop["transmitter"].waveform_prop["pulses"], at=60)
ranged = range_fft(baseband.conj(), window=range_window)
range_doppler = doppler_fft(ranged, window=dop_window)
绘制处理结果：
from scipy.constants import speed_of_light
max_range = ( # 计算最大探测距离
speed_of_light # 光速
* radar.radar_prop["receiver"].bb_prop["fs"] # 采样率
* radar.radar_prop["transmitter"].waveform_prop["pulse_length"] # 脉冲时长 Tc
/ radar.radar_prop["transmitter"].waveform_prop["bandwidth"]
/ 2
)
range_axis = np.linspace(0, max_range, radar.sample_prop["samples_per_pulse"], endpoint=False)
doppler_axis = np.linspace(
-1 / radar.radar_prop["transmitter"].waveform_prop["prp"][0] / 2,
1 / radar.radar_prop["transmitter"].waveform_prop["prp"][0] / 2,
radar.radar_prop["transmitter"].waveform_prop["pulses"],
endpoint=False,
)* speed_of_light / 60.5e9 / 2
fig = go.Figure()
fig.add_trace(
go.Heatmap(
z=20 * np.log10(np.mean(np.abs(range_doppler), axis=0)).T,
x=doppler_axis,
y=range_axis,
colorscale='Viridis',
zmin=-40,
zmax=0,
)
)
fig.update_layout(title='Range-Doppler map', xaxis=dict(title='Doppler (m/s)'), yaxis=dict(title='Range (m)'))
fig.show()
Angle-FFT
在进行 Angle-FFT 之前，先进行多普勒补偿。
# 获取 dop_bin 长度
dop_bin = np.arange(range_doppler.shape[-2])
# 计算相位差
delta_phi = np.exp(-1j * np.pi * (dop_bin - len(dop_bin) / 2) / len(dop_bin))
# 进行补偿
range_doppler[8:, :, :] *= delta_phi[None, :, None]
然后进行 Angle-FFT 处理：
def angle_fft(X, n=None, window=None):
if window is not None:
shape = [1] * X.ndim
shape[0] = X.shape[0]
X = X * window.reshape(shape)
output = np.fft.fftshift(np.fft.fft(X, n, axis=0), axes=0)
return output
angled = angle_fft(range_doppler, n=128)
绘制一个 Speed-Azimuth 图看效果。下图主要看 y 轴对不对，x 轴我标歪了。
azimuth = np.linspace(-1, 1, 128)
azimuth = np.arcsin(azimuth) / np.pi * 180
fig = go.Figure()
z = 20 * np.log10(np.abs(angled[:, :, 0]) + 0.001)
fig.add_trace(go.Heatmap(x=doppler_axis, y=azimuth, z=z, colorscale='Rainbow'))
fig.update_layout(
title='Speed-Azimuth Map',
height=400,
width=700,
scene=dict(
xaxis=dict(title='Doppler (m/s)'),
yaxis=dict(title='Azimuth (deg)'),
zaxis=dict(title='Amplitude (dB)'),
aspectmode='data',
),
margin=dict(l=0, r=0, b=0, t=40),
)
fig.show()
