In [1]:
import radarsimpy
radarsimpy.__version__
Out[1]:
'10.2.0'

FMCW Radar DoA Estimation¶



This is an example of using DoA algorithms (MUSIC & Esprit) with RadarSimPy.

Setup FMCW radar¶

Transmitter¶

Setup the basic transmitter parameters through Transmitter module.

The following table lists the parameters in this example.

Parameter Variable in Transmitter Value
Frequency ($f$) f [60e9, 61e9] GHz
Time ($T$) t 16e-6 s
Transmitted power ($P_t$) tx_power 15 dBm
Pulse repetition period ($PRP$) prp 40 us
Number of pulses pulses 512

Receiver¶

Setup the receiver parameters through Receiver module.

The parameters of the receiver are listed in the table below.

Parameter Variable in Receiver Value
Sampling rate ($f_s$) fs 20 Msps
Noise figure ($NF$) noise_figure 8 dB
RF gain/loss ($G_{rf}$) rf_gain 20 dB
Load resistor ($R_L$) load_resistor 500 $\Omega$
Baseband voltage gain ($G_{BB}$) baseband_gain 30 dB
In [2]:
import numpy as np
from radarsimpy import Radar, Transmitter, Receiver
wavelength = 3e8 / 60.5e9

channels = []

N_tx = 2
N_rx = 64

channels.append(
    dict(
        location=(0, -N_rx/2*wavelength/2, 0),
    ))

channels.append(
    dict(
        location=(0, wavelength*N_rx/2-N_rx/2*wavelength/2, 0),
    ))

tx = Transmitter(f=[61e9, 60e9],
                 t=[0, 16e-6],
                 tx_power=15,
                 prp=40e-6,
                 pulses=512,
                 channels=channels)

channels = []
for idx in range(0, N_rx):
    channels.append(
        dict(
            location=(0, wavelength/2*idx-(N_rx-1) * wavelength/4, 0),
        ))

rx = Receiver(fs=20e6,
              noise_figure=8,
              rf_gain=20,
              load_resistor=500,
              baseband_gain=30,
              channels=channels)

radar = Radar(transmitter=tx, receiver=rx)

There are 2 transmitters and 64 receivers, the synthesized array is a 1 x 128 uniform linear array.

Note: The modulation is not considered here. Assume the 2 transmitters are orthogonal.

In [3]:
import plotly.graph_objs as go
from IPython.display import Image

fig = go.Figure()
fig.add_trace(
    go.Scatter(x=radar.transmitter.locations[:, 1]/wavelength,
               y=radar.transmitter.locations[:, 2]/wavelength,
               mode='markers',
               name='Transmitter',
               opacity=0.7,
               marker=dict(size=10)))

fig.add_trace(
    go.Scatter(x=radar.receiver.locations[:, 1]/wavelength,
               y=radar.receiver.locations[:, 2]/wavelength,
               mode='markers',
               opacity=1,
               name='Receiver'))

fig.update_layout(
    title='Array configuration',
    xaxis=dict(title='y (λ)'),
    yaxis=dict(title='z (λ)',
               scaleanchor="x",
               scaleratio=1),
)

# fig.show()
Image(fig.to_image(format="jpg", scale=2))
Out[3]:

Targets¶

Create 3 targets at 40 m. The azimuth angles of the 3 targets are -5, -4, and 45 degrees relative to the radar.

In [4]:
true_theta = [-5, -4, 45]

target_1 = dict(location=(40*np.cos(np.radians(true_theta[0])), 40*np.sin(np.radians(true_theta[0])),
                0), speed=(0, 0, 0), rcs=10, phase=0)
target_2 = dict(location=(40*np.cos(np.radians(true_theta[1])), 40*np.sin(np.radians(true_theta[1])),
                0), speed=(0, 0, 0), rcs=10, phase=0)
target_3 = dict(location=(40*np.cos(np.radians(true_theta[2])), 40*np.sin(np.radians(true_theta[2])),
                0), speed=(0, 0, 0), rcs=10, phase=0)

targets = [target_1, target_2, target_3]

Simulate Baseband Signals¶

Use the simulator module to simulate the baseband samples. The user can choose between Python engine simpy or C++ engine simc.

The output baseband data is a 3-D matrix:

$[channels, pulses, ADC~samples]$

In [5]:
from radarsimpy.simulator import simc

bb_data = simc(radar, targets, noise=True)
time_matrix = bb_data['timestamp']
baseband = bb_data['baseband']

Radar signal processing¶

Range-Doppler processing¶

In [6]:
from scipy import signal
import radarsimpy.processing as proc

range_window = signal.chebwin(radar.samples_per_pulse, at=80)
doppler_window = signal.chebwin(radar.transmitter.pulses, at=60)

range_doppler = proc.range_doppler_fft(
    baseband, rwin=range_window, dwin=doppler_window)

Range-Doppler average map

In [7]:
max_range = (3e8 * radar.receiver.fs *
             radar.transmitter.pulse_length /
             radar.transmitter.bandwidth / 2)
range_axis = np.flip(np.linspace(
    0, max_range, radar.samples_per_pulse, endpoint=False))

unambiguous_speed = 3e8 / radar.transmitter.prp[0] / \
    radar.transmitter.fc_vect[0] / 2
doppler_axis = np.linspace(
    -unambiguous_speed, 0, radar.transmitter.pulses, endpoint=False)

fig = go.Figure()
fig.add_trace(go.Surface(x=range_axis, y=doppler_axis, z=20 *
              np.log10(np.mean(np.abs(range_doppler), axis=0)), colorscale='Rainbow'))

fig.update_layout(
    title='Range Doppler',
    height=600,
    scene=dict(
        xaxis=dict(title='Range (m)'),
        yaxis=dict(title='Velocity (m/s)'),
        zaxis=dict(title='Amplitude (dB)'),
    ),
    margin=dict(l=0, r=0, b=60, t=100),
    legend=dict(orientation='h'),
)

# fig.show()
Image(fig.to_image(format="jpg", scale=2))
Out[7]:

Find the beam vector of the peak and create the covariance matrix.

In [8]:
from scipy import linalg

det_idx = [np.argmax(np.mean(np.abs(range_doppler[:, 0, :]), axis=0))]

bv = range_doppler[:, 0, det_idx[0]]
bv = bv/linalg.norm(bv)

snapshots = 20

bv_snapshot = np.zeros((N_tx*N_rx-snapshots, snapshots), dtype=complex)

for idx in range(0, snapshots):
    bv_snapshot[:, idx] = bv[idx:(idx+N_tx*N_rx-snapshots)]

covmat = np.cov(bv_snapshot.conjugate())

FFT of the beam vector. Two peaks can be seen, and it is not able to discriminate the targets at -5 and -4 degrees.

In [9]:
from scipy import fft
from scipy import signal

fft_spec = 20 * \
    np.log10(np.abs(fft.fftshift(fft.fft(bv.conjugate(), n=1024))))

fig = go.Figure()

fig.add_trace(go.Scatter(x=np.arcsin(np.linspace(-1, 1, 1024, endpoint=False))/np.pi*180,
                         y=fft_spec,
                         name='FFT')
              )

fig.update_layout(
    title='FFT',
    yaxis=dict(title='Amplitude (dB)'),
    xaxis=dict(title='Angle (deg)'),
    margin=dict(l=10, r=10, b=10, t=40),
)

# fig.show()
Image(fig.to_image(format="jpg", scale=2))
Out[9]: