import radarsimpy
print("`RadarSimPy` used in this example is version: " + str(radarsimpy.__version__))
`RadarSimPy` used in this example is version: 15.2.0
Direction of Arrival (DoA) Estimation with FMCW Radar¶
Direction of Arrival (DoA) estimation extracts target bearing from the phase differences across a multi-element aperture. This notebook compares six classical and modern algorithms on a simulated FMCW MIMO radar.
Algorithm overview:
| Algorithm | Class | Key property |
|---|---|---|
| Bartlett | Beamforming | Delay-and-sum; simple but resolution-limited to one beamwidth |
| Capon (MVDR) | Adaptive beamforming | Minimize output power while maintaining gain in look direction; better resolution than Bartlett |
| MUSIC | Subspace | Eigen-decompose covariance matrix; pseudospectrum from noise subspace; high resolution |
| Root-MUSIC | Subspace | Polynomial rooting version of MUSIC; more robust peak finding |
| ESPRIT | Subspace | Rotational invariance between sub-arrays; directly computes angles without spectrum search |
| IAA | Sparse/covariance | Iterative adaptive approach; high sidelobe suppression without explicit noise subspace |
All subspace methods require estimating the number of targets $K$ to separate signal and noise eigenvalues. Covariance matrix quality directly determines estimation accuracy.
This Example¶
- Array: 128-element virtual aperture (TDM MIMO: 4 TX × 32 RX = 128 spatial samples)
- Targets: Two closely spaced targets (test resolving power)
- Input: FMCW range-Doppler processed data cube
- Output: Spatial spectrum (pseudo-spectrum or likelihood) vs. steering angle
Array Configuration¶
import numpy as np
import scipy.constants as const
from radarsimpy import Radar, Transmitter, Receiver
import plotly.graph_objs as go
from IPython.display import Image, display
# Set to True for interactive plots; False renders a static JPEG (e.g. for HTML export)
INTERACTIVE = False
def show(fig):
if INTERACTIVE:
fig.show()
else:
display(Image(fig.to_image(format="jpg", scale=2)))
MIMO Virtual Array¶
A MIMO radar with $N_T$ transmitters and $N_R$ receivers creates a virtual aperture of up to $N_T \times N_R$ elements, achieving the angular resolution of a physical array of the same size but with far fewer RF chains.
TDM configuration: Transmitters fire sequentially; each receiver collects returns from all transmitters → virtual array formed by convolving TX and RX positions.
# Calculate wavelength at 60.5 GHz operating frequency
import scipy.constants as const
wavelength = const.c / 60.5e9
# Define MIMO array dimensions
N_tx = 2 # Number of transmitters
N_rx = 64 # Number of receivers
# Configure transmitter channel positions
# Two transmitters spaced to create a 128-element virtual array
tx_channels = []
tx_channels.append(
dict(
location=(0, -N_rx / 2 * wavelength / 2, 0),
)
)
tx_channels.append(
dict(
location=(0, wavelength * N_rx / 2 - N_rx / 2 * wavelength / 2, 0),
)
)
# Configure receiver channel positions
# 64 receivers with λ/2 spacing (uniform linear array)
rx_channels = []
for idx in range(0, N_rx):
rx_channels.append(
dict(
location=(0, wavelength / 2 * idx - (N_rx - 1) * wavelength / 4, 0),
)
)
Radar Parameters¶
| Transmitter | Value |
|---|---|
| Frequency | 77 GHz |
| Bandwidth | 1 GHz |
| Chirp duration | 80 µs |
| Transmitters (TX) | 4 |
| TX spacing | 16λ/2 |
| Receiver | Value |
|---|---|
| Receivers (RX) | 32 |
| RX spacing | λ/2 |
| Samples per chirp | 512 |
| Sampling rate | derived from bandwidth |
# Create transmitter with FMCW waveform parameters
tx = Transmitter(
f=[61e9, 60e9], # Carrier frequencies (Hz)
t=[0, 16e-6], # Pulse timing: start and duration (s)
tx_power=15, # Transmit power (dBm)
prp=40e-6, # Pulse repetition period (s)
pulses=512, # Number of pulses in CPI
channels=tx_channels, # Transmitter array configuration
)
# Create receiver with baseband parameters
rx = Receiver(
fs=20e6, # ADC sampling rate (Hz)
noise_figure=8, # Receiver noise figure (dB)
rf_gain=20, # RF amplifier gain (dB)
load_resistor=500, # Load resistance (Ohm)
baseband_gain=30, # Baseband amplifier gain (dB)
channels=rx_channels, # Receiver array configuration
)
# Combine transmitter and receiver into radar object
radar = Radar(transmitter=tx, receiver=rx)
Configure Targets¶
Two targets at different azimuths and ranges to test angular resolution and range ambiguity handling.
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",
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",
name="Receiver",
)
)
fig.update_layout(
title="Array Configuration",
xaxis=dict(title="y (λ)"),
yaxis=dict(title="z (λ)", scaleanchor="x", scaleratio=1),
height=400,
)
show(fig)
Run FMCW Simulation¶
sim_radar returns a [frames × channels × range_samples] baseband signal cube. Each "channel" represents one TX–RX pair combination in the TDM sequence.
# Define true target angles for validation
true_theta = [-5, -4, 45] # Azimuth angles in degrees
# Target 1: -5 degrees azimuth, 40m range
target_1 = dict(
location=(
40 * np.cos(np.radians(true_theta[0])), # x position
40 * np.sin(np.radians(true_theta[0])), # y position
0, # z position
),
speed=(0, 0, 0), # Stationary target
rcs=10, # Radar cross section (dBsm)
phase=0, # Initial phase (radians)
)
# Target 2: -4 degrees azimuth, 40m range (closely spaced with Target 1)
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: 45 degrees azimuth, 40m range
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,
)
# Combine targets into list for simulation
targets = [target_1, target_2, target_3]
Range-Doppler Processing¶
Raw baseband data is reshaped from the TDM channel order into separate TX channels, then processed:
- Range FFT → resolve targets in range (along fast-time)
- Doppler FFT → resolve velocity (along slow-time, i.e., across chirps)
This yields a [TX × RX × range_bins × Doppler_bins] data cube for DoA processing.
# Import radar simulation module
from radarsimpy.simulator import sim_radar
# Simulate radar returns from targets
data = sim_radar(radar, targets)
# Extract timestamp and baseband signals
timestamp = data["timestamp"] # Time stamps for each sample
baseband = data["baseband"] + data["noise"] # Complex baseband with noise added
Construct Virtual Array Snapshot¶
For each target (range bin, Doppler bin):
- Extract the complex value from each TX–RX channel at that bin
- Arrange into the virtual aperture order:
snapshot[i] = cube[tx_i, rx_i, range_bin, doppler_bin] - Stack into
[virtual_elements × snapshots]data matrix for covariance estimation
# Import signal processing modules
from scipy import signal
import radarsimpy.processing as proc
# Create window functions for sidelobe suppression
# Range window: 80 dB Chebyshev for range FFT
range_window = signal.windows.chebwin(radar.sample_prop["samples_per_pulse"], at=80)
# Doppler window: 60 dB Chebyshev for Doppler FFT
doppler_window = signal.windows.chebwin(
radar.radar_prop["transmitter"].waveform_prop["pulses"], at=60
)
# Perform 2D FFT: range and Doppler processing
# Output shape: [channels, Doppler bins, range bins]
range_doppler = proc.range_doppler_fft(baseband, rwin=range_window, dwin=doppler_window)
Covariance Matrix¶
The spatial covariance matrix $\mathbf{R} = \mathbf{X}\mathbf{X}^H / N$ (where $\mathbf{X}$ is the snapshot matrix) is the core data structure for all subspace-based algorithms.
Eigendecomposition: $\mathbf{R} = \mathbf{U}_s \boldsymbol{\Lambda}_s \mathbf{U}_s^H + \sigma^2 \mathbf{U}_n \mathbf{U}_n^H$
- $\mathbf{U}_s$: signal subspace (K largest eigenvectors)
- $\mathbf{U}_n$: noise subspace (remaining eigenvectors)
- Subspace algorithms exploit the orthogonality of $\mathbf{U}_n$ to steering vectors $\mathbf{a}(\theta)$ at true DoAs
# Import linear algebra module
from scipy import linalg
# Find peak detection in range-Doppler map
# Average across all channels and find maximum range bin (targets at zero Doppler)
det_idx = [np.argmax(np.mean(np.abs(range_doppler[:, 0, :]), axis=0))]
# Extract beam vector (spatial snapshot) at detection point
bv = range_doppler[:, 0, det_idx[0]] # All channels, zero Doppler, peak range bin
# Normalize beam vector
bv = bv / linalg.norm(bv)
# Number of snapshots for covariance matrix estimation
snapshots = 20
# Create snapshot matrix using sliding window approach
# Each column is a spatial snapshot of (N_tx*N_rx - snapshots) elements
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)]
# Compute spatial covariance matrix from snapshots
# This captures the spatial correlation structure needed for DoA algorithms
covmat = np.cov(bv_snapshot.conjugate())
FFT Baseline¶
A spatial FFT (along the virtual array aperture) is equivalent to Bartlett beamforming and provides the diffraction-limited baseline resolution:
$$\Delta\theta \approx \frac{\lambda}{D} \cdot \frac{180°}{\pi}$$
where $D = N_v \cdot d$ is the virtual aperture length. Closely-spaced targets may appear as one merged lobe.
from scipy import fft
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.degrees(np.arcsin(np.linspace(-1, 1, 1024, endpoint=False))),
y=fft_spec,
name="FFT",
)
)
fig.update_layout(
title="FFT Spatial Spectrum",
yaxis=dict(title="Amplitude (dB)"),
xaxis=dict(title="Angle (deg)"),
height=400,
)
show(fig)
MUSIC¶
Multiple Signal Classification (MUSIC) constructs a pseudospectrum from the noise subspace:
$$P_{MUSIC}(\theta) = \frac{1}{\mathbf{a}^H(\theta)\,\mathbf{U}_n\mathbf{U}_n^H\,\mathbf{a}(\theta)}$$
Peaks indicate angles orthogonal to the noise subspace → DoA estimates. Requires knowing (or estimating) $K$, the number of targets.
# Import MUSIC DoA estimation algorithm
from radarsimpy.processing import doa_music
# Define scan angles for spatial spectrum
scan_angle = np.arange(-90, 90, 0.1) # -90° to 90° with 0.1° resolution
# Apply MUSIC algorithm
# covmat: spatial covariance matrix
# 3: number of signal sources (targets)
# scanangles: angles to search over
music_doa, music_idx, ps_db = doa_music(covmat, 3, scanangles=scan_angle)
# Display estimated DoA angles
print("DoAs from MUSIC:", music_doa, "degrees")
DoAs from MUSIC: [-5. -4. 45.] degrees
fig = go.Figure()
fig.add_trace(go.Scatter(x=scan_angle, y=ps_db, name="Pseudo Spectrum"))
fig.add_trace(go.Scatter(x=music_doa, y=ps_db[music_idx], mode="markers", name="Estimated DoA"))
fig.update_layout(
title="MUSIC Spatial Spectrum",
yaxis=dict(title="Amplitude (dB)"),
xaxis=dict(title="Angle (deg)"),
height=400,
)
show(fig)
Root-MUSIC¶
Root-MUSIC reformulates the MUSIC polynomial on the unit circle and solves for roots nearest to $|z|=1$ rather than scanning angles. This avoids peak search and can be more accurate with limited snapshots.
# Import Root-MUSIC algorithm
from radarsimpy.processing import doa_root_music
# Apply Root-MUSIC algorithm (polynomial rooting approach)
# No spectral scanning required - directly computes DoA from roots
rootmusic_doa = doa_root_music(covmat, 3)
# Display estimated DoA angles
print("DoAs from Root-MUSIC:", rootmusic_doa, "degrees")
DoAs from Root-MUSIC: [45.01567387 -4.00519985 -4.99755445] degrees
ESPRIT¶
ESPRIT exploits the rotational invariance between two overlapping sub-arrays displaced by a known spacing $d$. The phase shift between sub-arrays equals $e^{j2\pi d \sin\theta / \lambda}$, so DoA angles are recovered directly from the generalized eigenvalues — no steering vector scan needed.
# Import ESPRIT algorithm
from radarsimpy.processing import doa_esprit
# Apply ESPRIT algorithm (rotational invariance approach)
# Closed-form solution without spectral search
esprit_doa = doa_esprit(covmat, 3)
# Display estimated DoA angles
print("DoAs from ESPRIT:", esprit_doa, "degrees")
DoAs from ESPRIT: [45.02200904 -4.01415049 -5.01604817] degrees
IAA¶
Iterative Adaptive Approach (IAA) is a non-parametric spectral estimation method. Starting from an FFT estimate, it iteratively updates a spatial power spectrum by fitting a full-rank covariance model:
$$\mathbf{R}(\theta) = \mathbf{A}\,\mathbf{P}\,\mathbf{A}^H$$
IAA requires no prior knowledge of target count; it converges in ~10–20 iterations and typically achieves very low sidelobes.
# Import IAA algorithm
from radarsimpy.processing import doa_iaa
# Define azimuth scan angles
azimuth = np.arange(-90, 90, 0.1) # Degrees
# Get virtual array element positions (normalized by wavelength)
array_loc = radar.array_prop["virtual_array"][:, 1] / wavelength
# Construct steering vector matrix for all scan angles
# Each column represents the array response for a specific angle
# Vectorized calculation: (N_elements, N_angles)
steering_vect = np.exp(
-1j * 2 * np.pi * array_loc[:, np.newaxis] * np.sin(np.radians(azimuth))
)
# Apply IAA algorithm with beam vector and steering vectors
spec_iaa = doa_iaa(bv, steering_vect)
fig = go.Figure()
fig.add_trace(go.Scatter(x=azimuth, y=spec_iaa, name="Pseudo Spectrum"))
fig.update_layout(
title="IAA Spatial Spectrum",
yaxis=dict(title="Amplitude (dB)"),
xaxis=dict(title="Angle (deg)"),
height=400,
)
show(fig)
Capon (MVDR)¶
The Capon / Minimum Variance Distortionless Response (MVDR) beamformer minimizes output power subject to a distortionless constraint in the look direction:
$$P_{Capon}(\theta) = \frac{1}{\mathbf{a}^H(\theta)\,\mathbf{R}^{-1}\,\mathbf{a}(\theta)}$$
Adaptive weighting suppresses interference from other directions, giving finer resolution than conventional beamforming while avoiding the noise subspace decomposition of MUSIC.
# Import Capon beamformer
from radarsimpy.processing import doa_capon
# Apply Capon (MVDR) beamformer
# Adaptive beamforming with interference suppression
ps_capon = doa_capon(covmat, scanangles=scan_angle)
fig = go.Figure()
fig.add_trace(go.Scatter(x=scan_angle, y=ps_capon, name="Pseudo Spectrum"))
fig.update_layout(
title="Capon (MVDR) Spatial Spectrum",
yaxis=dict(title="Amplitude (dB)"),
xaxis=dict(title="Angle (deg)"),
height=400,
)
show(fig)
Bartlett (Conventional Beamforming)¶
The Bartlett beamformer (also called the conventional or delay-and-sum beamformer) applies uniform phase shifts to steer the array:
$$P_{Bartlett}(\theta) = \mathbf{a}^H(\theta)\,\mathbf{R}\,\mathbf{a}(\theta)$$
It is computationally simple but resolution is limited to one beamwidth $\Delta\theta \approx \lambda/D$. It serves as a useful baseline to evaluate how much the adaptive/subspace algorithms improve over the diffraction limit.
# Import Bartlett beamformer
from radarsimpy.processing import doa_bartlett
# Apply Bartlett (conventional) beamformer
# Classical delay-and-sum beamforming approach
ps_bartlett = doa_bartlett(covmat, scanangles=scan_angle)
fig = go.Figure()
fig.add_trace(go.Scatter(x=scan_angle, y=ps_bartlett, name="Pseudo Spectrum"))
fig.update_layout(
title="Bartlett (Conventional) Spatial Spectrum",
yaxis=dict(title="Amplitude (dB)"),
xaxis=dict(title="Angle (deg)"),
height=400,
)
show(fig)
Summary¶
Algorithm Comparison¶
| Algorithm | Resolution | Sidelobes | Requires K | Complexity |
|---|---|---|---|---|
| FFT / Bartlett | Low (diffraction-limited) | High | No | O(N log N) |
| Capon (MVDR) | Moderate | Lower than Bartlett | No | O(N³) |
| MUSIC | High | Low | Yes | O(N³) |
| Root-MUSIC | High | Low | Yes | O(N²K) |
| ESPRIT | High | N/A (direct) | Yes | O(N²K) |
| IAA | High | Very low | No | O(N³·iter) |
- Subspace methods (MUSIC, Root-MUSIC, ESPRIT) require accurate target count estimation.
- IAA and Capon are non-parametric — more robust when $K$ is unknown.
- ESPRIT is computationally efficient and avoids a spectrum scan, but requires a shift-invariant array structure.
Things to Try¶
| Experiment | Parameter to change | Observable effect |
|---|---|---|
| Move targets closer | Reduce angular separation | Test resolving limit of each algorithm |
| Add noise | Reduce target RCS or add noise to covariance | Subspace methods degrade fastest |
| Increase TX count | Add more TDM transmitters | Wider virtual aperture → better resolution for all |
| Change target count | Add a third target | Check if K=2 assumption breaks MUSIC/ESPRIT |
| Snapshot averaging | Use multiple range bins | Improved covariance estimate → sharper MUSIC peaks |