FMCW Radar

In this example we simulate a simple frequency modulated continuous wave (FMCW) radar system.

You can also find this example as a single file in the GitHub repository.

Below we have a very simplistic image of a radar system, it consists of a signal generator, transmitter/reciever chain and directed antennas. The fundamental working principle of a FMCW radar is that if we multiply a chirp signal (linearly frequency modulated sinusoid) with a delayed version of itself, we get a signal that has at each time two frequency components (sum and diffrerence). The low frequency component then is directly proportional to the phase shift between the two signals and this the time delay. This can be used to reconstruct the target distance from the frequency domain spectrum.

The figure below shows the block diagram of a very basic FMCW radar system consisting of a chirp source (sinusoid with a linearly time dependent frequency), a delay (to model the signal propagation), a multiplier (to model a mixer), a low pass filter and of course blocks for time series and frequency domain data recording.

Lets start by importing the requred classes from PathSim:

import numpy as np
import matplotlib.pyplot as plt

# Apply PathSim docs matplotlib style for consistent, theme-friendly figures
plt.style.use('../pathsim_docs.mplstyle')

from pathsim import Simulation, Connection

# From the block library
from pathsim.blocks import Multiplier, Scope, Adder, Spectrum, Delay, ChirpPhaseNoiseSource, ButterworthLowpassFilter

Now lets define the system parameters starting with the ChirpPhaseNoiseSource, which is defined by the starting frequency, the bandwith and the ramp duration.

# Natural constants (approximately)
c0 = 3e8

# Chirp parameters
B, T, f_min = 4e9, 5e-7, 1e9

# Delay for target emulation
tau = 2e-9

# For reference, the expected target distance
R = c0 * tau / 2

# And the corresponding frequency
f_trg = 4 * R * B / (T * c0)

Next we can build the system by defining the blocks and their connections:

Src = ChirpPhaseNoiseSource(f0=f_min, BW=B, T=T)
Add = Adder()
Dly = Delay(tau)
Mul = Multiplier()
Lpf = ButterworthLowpassFilter(f_trg*3, n=2)
Spc = Spectrum(
    freq=np.logspace(6, 10.5, 500),
    labels=["chirp", "delay", "mixer", "lpf"]
)
Sco = Scope(
    labels=["chirp", "delay", "mixer", "lpf"]
)

# Collecting the blocks in a list
blocks = [Src, Add,  Dly, Mul, Lpf, Spc, Sco]

# Connections between the blocks
connections = [
    Connection(Src, Add[0]),
    Connection(Add, Dly, Mul, Sco, Spc),
    Connection(Dly, Mul[1], Sco[1], Spc[1]),
    Connection(Mul, Lpf, Sco[2], Spc[2]),
    Connection(Lpf, Sco[3], Spc[3])
]

Now we are ready to initialize the simulation and run it for some time. Here it makes sense to run it for the duration of one chirp period:

# Initialize simulation
Sim = Simulation(blocks, connections, dt=1e-11, log=True)

# Run simulation for one chirp period
Sim.run(T)

10:58:08 - INFO - LOGGING (log: True)
10:58:08 - INFO - BLOCKS (total: 7, dynamic: 3, static: 4, eventful: 1)
10:58:08 - INFO - GRAPH (nodes: 7, edges: 13, alg. depth: 3, loop depth: 0, runtime: 0.165ms)
10:58:08 - INFO - STARTING -> TRANSIENT (Duration: 0.00s)
10:58:08 - INFO - --------------------   1% | 0.2s<18.1s | 2729.7 it/s
10:58:09 - INFO - #-------------------   6% | 1.2s<16.8s | 2787.5 it/s
10:58:10 - INFO - ##------------------  11% | 2.2s<17.1s | 2575.5 it/s
10:58:11 - INFO - ###-----------------  17% | 3.2s<14.9s | 2779.7 it/s
10:58:12 - INFO - ####----------------  20% | 3.7s<14.5s | 2759.6 it/s
10:58:13 - INFO - #####---------------  25% | 4.7s<14.2s | 2621.7 it/s
10:58:14 - INFO - ######--------------  30% | 5.7s<12.5s | 2767.6 it/s
10:58:15 - INFO - #######-------------  36% | 6.7s<13.1s | 2430.1 it/s
10:58:16 - INFO - ########------------  40% | 7.4s<10.9s | 2754.6 it/s
10:58:17 - INFO - #########-----------  45% | 8.4s<10.1s | 2711.4 it/s
10:58:18 - INFO - ##########----------  50% | 9.4s<9.5s | 2584.1 it/s
10:58:19 - INFO - ###########---------  56% | 10.4s<8.0s | 2748.4 it/s
10:58:19 - INFO - ############--------  60% | 11.1s<7.2s | 2793.8 it/s
10:58:20 - INFO - #############-------  65% | 12.1s<6.2s | 2794.1 it/s
10:58:21 - INFO - ##############------  70% | 13.1s<5.4s | 2726.2 it/s
10:58:22 - INFO - ###############-----  76% | 14.1s<4.4s | 2723.9 it/s
10:58:23 - INFO - ################----  80% | 14.8s<3.7s | 2727.7 it/s
10:58:24 - INFO - #################---  85% | 15.8s<2.7s | 2713.7 it/s
10:58:25 - INFO - ##################--  90% | 16.8s<1.7s | 2768.2 it/s
10:58:26 - INFO - ###################-  96% | 17.8s<0.7s | 2708.4 it/s
10:58:27 - INFO - #################### 100% | 18.6s<--:-- | 2715.6 it/s
10:58:27 - INFO - FINISHED -> TRANSIENT (total steps: 50000, successful: 50000, runtime: 18556.05 ms)

{'total_steps': 50000,
 'successful_steps': 50000,
 'runtime_ms': 18556.04849899828}

Lets have a look at the time series data first. We can do this by calling the plot method of the scope instance. Here we have four traces which we can toggle on and off.

#plot the recording of the scope
Sco.plot()

(<Figure size 960x480 with 1 Axes>, <Axes: xlabel='time [s]'>)

It only really gets interesting in the frequency domain. So lets look at the spectrum block (and scale it logarithmically):

#plot the spectrum
Spc.plot()
Spc.ax.set_xscale("log")
Spc.ax.set_yscale("log")
plt.show()

In the spectrum the trace of interest is the output of the low pass filter (purple trace) which is intended to select the signal component that represents the delay, or radar distance. The position of the peak corresponds directly to the target distance, represented by the delay block.

Isolating the spectrum of the lowpass filter and adding the expected target distance (as a frequency) to the plot shows that the FMCW radar system can indeed correctly resolve the range:

#plot the spectrum
Spc.plot()
Spc.ax.set_xscale("log")
Spc.ax.set_yscale("log")

Spc.ax.axvline(f_trg, ls="--", c="grey")

plt.show()