Phase gradient autofocus

import torch
import torchbp
import matplotlib.pyplot as plt
from scipy.signal import get_window
import numpy as np

if torch.cuda.is_available():
    device = "cuda"
else:
    device = "cpu"
print("Device:", device)

Device: cpu

Generate synthetic radar data

nr = 100 # Range points
ntheta = 128 # Azimuth points
nsweeps = 128 # Number of measurements
fc = 6e9 # RF center frequency
bw = 100e6 # RF bandwidth
tsweep = 100e-6 # Sweep length
fs = 1e6 # Sampling frequency
nsamples = int(fs * tsweep) # Time domain samples per sweep

# Imaging grid definition. Azimuth angle "theta" is sine of radians. 0.2 = 11.5 degrees.
grid_polar = {"r": (90, 110), "theta": (-0.2, 0.2), "nr": nr, "ntheta": ntheta}

target_pos = torch.tensor([[100, 0, 0], [105, 10, 0], [97, -5, 0], [102, -10, 0], [95, 5, 0]], dtype=torch.float32, device=device)
target_rcs = torch.tensor([1,1,1,1,1], dtype=torch.float32, device=device)
pos = torch.zeros([nsweeps, 3], dtype=torch.float32, device=device)
pos[:,1] = torch.linspace(-nsweeps/2, nsweeps/2, nsweeps) * 0.25 * 3e8 / fc

# Oversampling input data decreases interpolation errors
oversample = 3

# Modulation frequency in range direction to center the spectrum at DC
# for more accurate interpolation.
data_fmod = -torch.pi * (1 - (oversample-1) / oversample)

data = torchbp.util.generate_fmcw_data(target_pos, target_rcs, pos, fc, bw, tsweep, fs)
# Apply windowing function in range direction
wr = torch.tensor(get_window(("taylor", 3, 30), data.shape[-1])[None,:], dtype=torch.float32, device=device)
wa = torch.tensor(get_window(("taylor", 3, 30), data.shape[0])[:,None], dtype=torch.float32, device=device)
data = torch.fft.ifft(data * wa * wr, dim=-1, n=nsamples * oversample)

data_fmod_f = torch.exp(1j*data_fmod*torch.arange(data.shape[-1], device=device))[None,:]
data = data * data_fmod_f

data_db = 20*torch.log10(torch.abs(data)).detach()
m = torch.max(data_db)

plt.figure()
plt.imshow(data_db.cpu().numpy(), origin="lower", vmin=m-30, vmax=m, aspect="auto")
plt.xlabel("Range samples")
plt.ylabel("Azimuth samples");

Focused image without motion error

r_res = 3e8 / (2 * bw * oversample) # Range bin size in input data

# dealias=True removes range spectrum aliasing
img = torchbp.ops.backprojection_polar_2d(data, grid_polar, fc, r_res, pos, dealias=True, data_fmod=data_fmod)
img = img.squeeze(0) # Removes singular batch dimension
# Backprojection image has spectrum with DC at zero index.
# Shifting the spectrum shifts the DC to center bin.
# This makes the solved phase to have same order as the position vector
# Without shifting of the image, fftshift needs to be applied to
# the solved phase for it to be in the same order as the position vector.
# This doesn't affect the absolute value of the image.
img = torchbp.util.shift_spectrum(img)

img_db = 20*torch.log10(torch.abs(img)).detach()

m = torch.max(img_db)

extent = [*grid_polar["r"], *grid_polar["theta"]]

plt.figure()
plt.imshow(img_db.cpu().numpy().T, origin="lower", vmin=m-40, vmax=m, extent=extent, aspect="auto")
plt.xlabel("Range (m)")
plt.ylabel("Angle (sin radians)");

Create corrupted image

phase_error = torch.exp(1j*2*torch.pi*torch.linspace(-3, 3, ntheta, dtype=torch.float32, device=device)[None,:]**2)

plt.figure()
plt.plot(torch.angle(phase_error.squeeze()).cpu().numpy())
plt.xlabel("Azimuth sample")
plt.ylabel("Phase error (radians)")

img_corrupted = torch.fft.ifft(torch.fft.fft(img, dim=-1) * phase_error, dim=-1)

plt.figure()
plt.imshow(20*torch.log10(torch.abs(img_corrupted)).cpu().numpy().T, origin="lower", vmin=m-40, vmax=m, extent=extent, aspect="auto")
plt.xlabel("Range (m)")
plt.ylabel("Angle (sin radians)");

Phase gradient autofocus with phase difference estimator

img_pga, phi = torchbp.autofocus.pga(img_corrupted, remove_trend=False, estimator="pd")

plt.figure()
plt.imshow(20*torch.log10(torch.abs(img_pga)).cpu().numpy().T, origin="lower", vmin=m-40, vmax=m, extent=extent, aspect="auto")
plt.xlabel("Range (m)")
plt.ylabel("Angle (sin radians)");

plt.figure()
plt.plot(torch.angle(torch.exp(1j*phi)).cpu().numpy())
plt.xlabel("Azimuth samples")
plt.ylabel("Phase error (radians)");

Apply maximum likelihood phase gradient autofocus

img_pga, phi = torchbp.autofocus.pga(img_corrupted, remove_trend=False, estimator="ml")

plt.figure()
plt.imshow(20*torch.log10(torch.abs(img_pga)).cpu().numpy().T, origin="lower", vmin=m-40, vmax=m, extent=extent, aspect="auto")
plt.xlabel("Range (m)")
plt.ylabel("Angle (sin radians)");

plt.figure()
plt.plot(torch.angle(torch.exp(1j*phi)).cpu().numpy())
plt.xlabel("Azimuth samples")
plt.ylabel("Phase error (radians)");

Multiplying the solved phase with FFT of the corrupted image gives the focused image and taking inverse FFT gives the focused image. This should be identical to the image returned by pga_ml.

img_focused = torch.fft.ifft(torch.fft.fft(img_corrupted, dim=-1) * torch.exp(-1j*phi), dim=-1)

plt.figure()
plt.imshow(20*torch.log10(torch.abs(img_focused)).cpu().numpy().T, origin="lower", vmin=m-40, vmax=m, extent=extent, aspect="auto")
plt.xlabel("Range (m)")
plt.ylabel("Angle (sin radians)");