# How to simulate multiple targets in FMCW radar simulation in python?

I have simulated a FMCW radar for a single stationary target, now the next step in my project is to simulate multiple moving targets. For this I had thought of using doppler shift to show a difference in received frequecy, but for some reason after calculating the effective frequency the values of frequency for different targets remains the same. even though the doppler shift values are different.

Below is the code for Simulation of FMCW radar used by me(single target)

import numpy as np
import matplotlib.pyplot as plt
from math import pi

chirp_duration = 60e-6
f0 = 1e6
f1 = 10e6
no_samp = 8001
target_range = 3e3
bandwidth = f1 - f0
c = 3e8

f_step = (f1-f0)/no_samp

rf_axis = np.arange(f0, f1, f_step)
rt_axis = np.linspace(0, chirp_duration, len(rf_axis))

rchirp = np.cos(2 * np.pi * rf_axis * rt_axis)

plt.figure(figsize=(12, 3))
plt.xlabel("Time")
plt.ylabel("Frequency")
plt.plot(rt_axis, rchirp)
plt.show()

# transmtted signal
slope = bandwidth/no_samp
time_delay = 2*target_range/c
change_in_frequency = slope * time_delay
frequency_at_time = f0 + change_in_frequency

# Define f_axis
tf_axis = np.linspace(frequency_at_time, f1, int(f_step))

# Define t_axis
tt_axis = np.linspace(0, chirp_duration, int(f_step))
changed_axis = tt_axis + time_delay

# Create chirp
tchirp = np.cos(2 * np.pi * tf_axis * tt_axis)

# Plot chirp
plt.figure(figsize=(12, 3))
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.title('Transmitted Chirp')
plt.plot(changed_axis, tchirp)
plt.show()

For simulation of FMCW radar for multiple targets I have onl used the received chirp part as that's where the changes will happen.

import numpy as np
import matplotlib.pyplot as plt

chirp_duration = 60e-6
f0 = 1e6
f1 = 10e6
no_samp = 8001
c = 3e8

f_step = (f1 - f0) / no_samp
rf_axis = np.arange(f0, f1, f_step)
rt_axis = np.linspace(0, chirp_duration, len(rf_axis))

# Define target parameters
target_ranges = [3e3, 5e3, 8e3]  # Distances to the targets in meters
target_velocities = [10, 5, 15]  # Velocities of the targets in m/s

plt.figure(figsize=(12, 9))

# Plot the received chirp signal for each target
for i in range(len(target_ranges)):
# Calculate the Doppler shift for the current target
doppler_shift = 2 * target_velocities[i] / c * rf_axis

# Apply the Doppler shift to the chirp signal for the current target
rchirp = np.cos(2 * np.pi * (rf_axis + doppler_shift) * rt_axis)

# Create a subplot for the current target
plt.subplot(len(target_ranges), 1, i + 1)
plt.plot(rt_axis, rchirp)
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.title(f"Received Chirp for Target {i + 1}")

# Adjust the spacing between subplots
plt.tight_layout()

# Show the plot
plt.show()

The updated code is given below.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
chirp_duration = 60e-6
f0 = 1e6
f1 = 10e6
no_samp = 8001
c = 3e8

# Frequency axis
f_step = (f1 - f0) / no_samp
rf_axis = np.arange(f0, f1, f_step)
rt_axis = np.linspace(0, chirp_duration, len(rf_axis))

# Define target parameters
target_ranges = [3e3, 5e3, 8e3]  # Distances to the targets in meters
target_velocities = [30, 40, 50]  # Velocities of the targets in m/s

plt.figure(figsize=(12, 9))

# Plot the received chirp signal for each target
for i in range(len(target_ranges)):
# Calculate the phase shift for the current target
phase_shift = -2 * np.pi * (f1-f0) * 2 * target_velocities[i] / c

# Apply the phase shift to the received chirp signal for the
current target
rchirp = np.cos((2 * np.pi * rf_axis * rt_axis) + phase_shift)
print('rchirp',rchirp)

# Create a subplot for the current target
plt.subplot(len(target_ranges), 1, i + 1)
plt.plot(rt_axis, rchirp)
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.title(f"Received Chirp for Target {i + 1}")

# Adjust the spacing between subplots
plt.tight_layout()

# Show the plot
plt.show()

• @V.V.T or I thought that instead I'll just simulate by using the phase shift caused ny caused by the doppler shift due to change in frequency...new code is added in the question Dec 12, 2023 at 10:53
• You never calculate the received frequency; you simply plot the transmitted and received chirps. I guess you hope to see the Doppler shift by inspecting and comparing visually the graphs. Visual comparison is hardly possible even with the velocities of your updated code. With the target velocity of 50 m/s, $v=50$, the Doppler shift is $f_d = f_{received}-f_{transmitted} = 2v·f_{transmitted}/(c/n-v) \approx 2·50/(299792548/1.000293-50) \approx 100/299704594.54 \approx 0.00000033366$ Dec 13, 2023 at 5:20
• Can you notice the frequency variation as small as a 1/3 millionth part in the plot with bitmap resolution of even as great as 4K? To simulate the Doppler shift measurement, use triangular modulation. Plot the graph of chirp frequency vs. time and use frequencies at the rising/falling edges of this graph to compute. See formulas and discussion of what modulation type is required for the Doppler shift measurement, for example, in Frequency-Modulated Continuous-Wave Radar (FMCW Radar) tutorial radartutorial.eu/02.basics/… . Dec 13, 2023 at 5:24
• with the updated code the shift in the frequency is not quite visible but the phase shift is quite clearly visible Dec 13, 2023 at 6:46
• @V.V.T int he new code I didn't even calculate the doppler shift in frequency I directly calculated the phase shift yhat would be caused by the doppler shift. Now after reading your previous comments I am starting to think that the method followed by me is not right?? Dec 13, 2023 at 7:01

For students with special needs who prefer to learn via coding (like me and maybe you), the perfect option is to find a Jupyter notebook on FMCW radar basics. Having failed to find one, I'm writing an answer under the guise of Jupyter, with intermittent pieces of Python code, formulas, and explanatory text. Hope this helps.

Let us start with the problems which a FMCW radar is intended to solve. When characterizing detection-and-ranging technologies like lidar, radar, sonar, 'ranging' names a technique to measure distances. Regardless of medium (outer space, air, water) and carrier (light, radio waves, sound waves) used, we have a template code to process the sensor data:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

c = 10.0 # speed of carrier waves (light, radio, sound) in medium (outer space, air, water)
samplerate = 2000
# search web with 'radar equation' keyword to understand the meaning of the expression below
radar_eq_coeff = 1 #[squaremeter] = [antenna gain] * [cross section] * [antenna area] * [efficiency]

delayval = 2 * np.int_(distance * samplerate / c)
prefix = np.zeros(delayval)

duration = 5
time = np.linspace(0,duration,samplerate*duration, endpoint= False)

sqeaktime = 0.5
dist = 1.25

squeak = np.zeros(len(time))
# send 0.1 s pulse
for i in range(samplerate//10):
squeak[np.int_(samplerate*sqeaktime) + i] = 1.0

echo = recvd_sig(squeak, dist)[:len(time)]

plt.figure('Range: sqeak, echo')
plt.legend()
plt.show()

The only specification in this code is the radar equation coefficient used to give an amplitude of the signal reflected from the target. It may become important when studying radars, but for now this amplitude can be safely set to unity.

The parameter values are purely arbitrary and chosen only to provide signal waveforms conveniently shown with plots. The function recvd_sig(radar_sends, distance) is defined with the only purpose to represent the physical arrangement in the program flow: the radar transmitter sends the pulse (squeak) which is reflected from the target (echo = recvd_sig(squeak, dist)), the distance between the radar and the target dist = 1.25 is passed as a parameter in the recvd_sig function call. The coder is free to unwrap this call to inline operators in the course of further code development.

As this code is a learning aid, I ain't attaching the plot output.

Given the transmitted pulse and the received echo timing data, the distance is calculated as the delay between these signals times the speed of carrier wave. The 'ranging' part can be solved in this way.

Two distance values dist(t1), dist(t2) measured at times t1, t2 seem to enable us to calculate the velocity using the velocity definition: v = (dist(t2)-dist(t1))/(t2-t1). However, the precision of this velocity measurement technique is poor, and the velocity value measured in this way is not instantaneous by definition (is it v(t1)? or v(t2)? or v((t1+t2)/2)?). Detection-and-ranging techniques enable us to measure the target velocity directly, using the Doppler effect.

See a template code where the sensor data is processed for velocity measurement using the Doppler effect:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

c = 10.0
samplerate = 2000

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

c = 10.0
samplerate = 2000

# target with POSITIVE velocity WITHDRAWS from the radar
# I omit the echo amplitude change here
# for light waves, this is only valid for target_velocity << c
# the exact Doppler shift for light waves is sqrt((1+beta)/(1-beta)); beta = target_velocity / c
dopplershifted_sample_count = np.int_(sample_count * (c + target_velocity) / (c - target_velocity))

freq = 10
duration = 1
time = np.linspace(0,duration,samplerate*duration, endpoint= False)

wave = np.sin(2 * np.pi * freq * time)
velocity = -0.1 # target with negative velocity APPROACHES the radar

echo = recvd_sig(wave, velocity)[:len(time)]
# set signal durations to min(len(wave),len(echo))
if (len(echo) < len(time)):
time = time[:len(echo)]
wave = wave[:len(echo)]

# plot wave and echo
plt.figure('Doppler: wave, echo' )
plt.title('Doppler shift, (target velocity)/c = %f' % (velocity/c))
plt.legend()
plt.show()

With this great ratio value of 'target velocity'-to-'carrier wave speed' (abs(velocity/10) = 0.1/10 = 0.01) as in this piece of code, had c be the light speed, velocity would have to be 0.01*300,000km/s = 3,000km/s, that is, 375 times the velocity of an LEO satellite. Of course, the frequency variation for this great target velocity is visible by inspection in the 'Doppler shift, target velocity = 0.1' plot, but for Earth-bound targets you cannot expect that the effect visibility is also pronounced in the graphs, and it is not.

The comment near the recvd_sig function def is a stipulation about a relativistic Doppler effect. This comment can be ignored for measurements with Earth-bound objects, but it is necessary to respect possible rigorous reviewers of this post.

Also, I omit the radar-equation factor in the amplitude of the reflected wave which is not essential for our consideration; the amplitude of the reflected wave is unity.

So, we need a continuous wave ('CW' in the 'FMCW' abbreviation) to measure target velocities, but the pure-sine continuous wave cannot be used to measure distances ('ranging'). The solution is to use a MODULATED continuous wave.

Had we have a means to attach a parameter to the continuous wave which makes it possible to carry this data with the wave, we would be able to calculate the delay between the moment t1, when this piece of data is transmitted, and the moment t2, when this piece of data is received after performing a round trip from the radar (lidar, sonar) to the target and back to the radar (lidar, sonar). At the same time, the wave, being a continuous wave, makes it possible to measure the instantaneous velocity of the target.

As it turns out, the necessary parameter can be attached in the form of instantaneous frequency, using a FREQUENCY MODULATION technique ('FM' in the 'FMCW' abbreviation). It is what the 'chirps' in your code are for.

To be sure that marking the continuous wave followed by mark reading is possible, examine the code where an CW wave is first modulated (xmit_wave = np.sin(np.pi*2*(fcarrier*time + message_integral/samplerate))) and then demodulated (analytic_xmit = signal.hilbert(xmit_wave); inst_freq = np.diff(np.unwrap(np.angle(analytic_xmit)))*samplerate/(2*np.pi)), thus recovering the sawtooth waveform (message) used to frequency-modulate a sine wave sin(np.pi*2*(fcarrier*time):

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

samplerate = 1600
fcarrier = 15
duration = 5
fmax = fcarrier + 20
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
# a message used to frequency-modulate the continuous wave
message = (fmax-fcarrier)*(signal.sawtooth(2*np.pi*time)+1)/2.0

# integrate the frequency modulating signal (message) over time
# to arrive at the variable phase following the principle of frequency modulation
message_integral = np.cumsum(message)

# FM wave; chirp from fcarrier to fmax
xmit_wave = np.sin(np.pi*2*(fcarrier*time + message_integral/samplerate))

plt.figure('chirp')
plt.plot(time, xmit_wave)
plt.show()

# recover instantaneous frequency
# first make the signal analytic
analytic_xmit = signal.hilbert(xmit_wave)
# instantaneous frequency is the derivative of analytic signal's phase ('angle')
inst_freq = np.diff(np.unwrap(np.angle(analytic_xmit)))*samplerate/(2*np.pi)

plt.figure('instantaneous frequency')
plt.title('recovered msg; chirp from f_carrier = %g, to f_max = %g' % (fcarrier, fmax))
plt.plot(time, message, label='message')
plt.plot(time[:len(time)-1], inst_freq, label='inst_freq = message + f_carrier')
plt.legend(loc="right")
plt.show()

This code processes the transmitted wave and therefore serves only as a feasibility proof. In the working program, you have to process the signal from the radar receiver.

In the explanation that follows, I assume that FM modulation/demodulation is performed elsewhere, and the function recvd_signal(radar_sends, target_distance, target_velocity) returns a decoded FM-demodulated message which the transmitter has encoded in the radar-emitted chirp.

This code is intended to perform ranging and velocity measurement in one batch. Only the speed-of-light parameter value, 300,000 km/s, is real; selection of the other parameter values is dictated by the convenience considerations of graphic presentation. In reality, the carrier frequency is much greater than 1MHz. However, using this code, the learner can see how the target distance is calculated from measurement data of how the message sawtooth is shifted in a horizontal direction, and the target velocity is calculated from measurement data of how the message sawtooth is shifted in a vertical direction:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

c = 3E5 # 300,000 km/s
samplerate = 300000
carrier_frequency = 1E6 # 1MHz
# target with POSITIVE velocity WITHDRAWS from the radar
# for light waves, this Doppler shift formula is approximate and only valid for target_velocity << c
# the exact Doppler shift for light waves is sqrt((1+beta)/(1-beta)); beta = target_velocity / c
# Also, I omit the echo amplitude change here
delayval = 2 * np.int_(target_distance * samplerate / c)
reflected_message_prefix = np.zeros(delayval)
message_shifted = np.append(reflected_message_prefix, message_truncated) - carrier_frequency * 2 * target_velocity / (c - target_velocity)
return message_shifted

duration = 0.001 # 1ms
distance = 3 # 3km
time = np.linspace(0,duration,np.int_(samplerate*duration), endpoint= False)
sawtooth_freq = 5000 # 5KHz

message = (signal.sawtooth(2*np.pi*time*sawtooth_freq)+1)/2.0
velocity = -0.01 # -10m/s. Target with negative velocity APPROACHES the radar

echo = recvd_signal(message, distance, velocity)

# plot wave and echo
plt.figure('Ranging, Doppler measurements' )
plt.title('distance = %g km, velocity = %g km/s' % (distance, velocity))