Generating IF signal for FMCW radar

I am trying to generate a signal that represents mixer output of FMCW radar. I am using MATLAB and already tried out this tutorial. There, Phased Array System Toolbox is used to showcase how FMCW radar with sawtooth signal works. The part that is of interest for me right now is up until the range-doppler response. Now what I am trying to accomplish is basically getting the same result (as in pass my signal to plotResponse function and get a range-doppler map with a detection of my target), but I want to get what they call "dechirped" signal from the start - without generating transmitted, reflected and received signals. I've been looking through so many papers in order to determine what formula I need to use to get the correct IF signal, but the more I searched, the more different formulas I found. What I've got now is the formula from this TI webinar, which is So with that, I am trying to plot the response, here is my code:

c = 3e8;            %speed of light
range_max = 180;    %max detection range
tm = 6*range2time(range_max,c); %sweep time
%tm is 7.2e-6 s
bw = 200e6;         %sweep bandwidth
sweep_slope = bw/tm;
v_max = 150*1000/3600;  %target max velocity
fc = 77e9;          %radar frequency
lambda = c/fc;      %radar wavelength
fs = 72e6;          %sampling rate
%sampling rate based on ADC datasheet
chirps = 64;        %frame size
samples = ceil(tm*fs);  %samples in one chirp

%% target
R0 = 32; %range in meters
V = 40; %radial velocity, m/s

%%
t = 0; %time
mix = zeros(samples, chirps); %mixer output
for i=1:1:chirps
td = 2 * R0 / c; %round trip delay
phi0 = 4*pi*fc*R0/c; %inital phase

for j=1:1:samples
a = (-2*pi*fc*V*i*tm/c ...    %phase shift
-2*pi*(2*V*(fc+i*bw)/c + sweep_slope*td)*t);   %frequency
mix(j,i) = 0.5*cos(a);
t = t + 1/fs;
end
end

figure(1)
rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,...
'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,...
'RangeMethod','FFT','SweepSlope',sweep_slope,...
'RangeFFTLengthSource','Property','RangeFFTLength',2048,...
'DopplerFFTLengthSource','Property','DopplerFFTLength',256);

clf;
plotResponse(rngdopresp,mix);
axis([-v_max v_max 0 range_max])


The radar that I'm simulating here does't have quadrature channel, so I have to only form real signal (although I tired doing hilbert(...) and cos(...)+1i*sin(...) just to see how that would look). What I assume I get is the Data Cube similar to the one from the tutorial. But no matter how I fiddle with the formula, I can't get the Velocity right. Here's the plot:

The initial phase I tried to form myself, thinking it would be 2*pi*fc*td, where fc is 77 GHz carrier frequency and td is the time delay of 2*R/c.
So what I think is happening is that I don't shift phase properly, since, from my understanding, that is what velocity is primarily estimated upon. The Doppler shift is too small to have significant impact.
Interestingly, I tried to play with the range and velocity values of the target, and noticed, that changing V doesn't affect the target's position on the heatmap too much, it almost doesn't actually. But changing the R0 in increments of 1 makes the target detection on heatmap shift horizontally and cycle to the other side of the map.
So my question is: how do i properly form the IF signal to be able to display it correctly in MATLAB's range-doppler response map?
Also, I acknowledge that I might be plotting or forming the data cube wrong. Or that my math is plain wrong from the beginning. But I've been so desperately trying to solve this by myself for the past week, that I just want someone to show me the light already :)

• Why the sweep time parameter define as 6*(2*range_max/c);
– meni
Mar 22 '21 at 12:02
• @meni "In general, for an FMCW radar system, the sweep time should be at least 5 to 6 times the round trip time. This example uses a factor of 5.5." Source Mar 25 '21 at 4:32

1 Answer

The principle behind FMCW is that you transmit a chirped signal and receive a time delayed version of it after reflecting from a target. After mixing and filtering, the resulting signal is a sinusoid at a frequency that is a function of the target's range. This frequency is known as the "beat" frequency $$f_b$$. Thus, the dechirped signal will have a form of

$$x(t) = e^{j(2{\pi}f_bt + \, \phi)} = e^{j2{\pi}f_bt}\,e^{j{\phi}}$$

Where $$\phi$$ is a general phase term that we'll ignore for now as it won't effect the determination of the beat frequency, and therefore range. Also, let's not worry about Doppler for now...that's an additional phase term we can easily add later. We'll concentrate on the expression for the homodyned (mixed) signal.

Let the chirped signal we transmit be

$$s_{tx}(t) = e^{j\pi\frac{\beta}{\tau}t^2}$$

Where $$\beta$$ is the sweep bandwidth of the chirp and $$\tau$$ is the chirp's length, or pulse width. After reflecting from a target we receive the signal after some delay $$t_d$$, we have

$$s_{rx}(t) = e^{j\pi\frac{\beta}{\tau}(t - t_d)^2} = e^{j\pi\frac{\beta}{\tau}(t^2 - 2tt_d + t_d^2)}$$

After mixing $$s_{rx}(t)$$ with $$s_{tx}(t)$$, which is equivalent to a frequency shift, the higher order term containing $$e^{j\pi\frac{\beta}{\tau}t^2}$$ drops off and we're left with

$$x(t) = e^{j\pi\frac{\beta}{\tau}(-2tt_d + t_d^2)} = e^{-j\pi\frac{\beta}{\tau}2tt_d}\,e^{j\pi\frac{\beta}{\tau}t_d^2}$$

Now compare this with the first equation, paying attention to the first term, again ignoring the constant phase term. We can then equate the phase functions $$-\pi\frac{\beta}{\tau}2tt_d = 2{\pi}f_bt$$

So that then we have

$$f_b = -\frac{\beta}{\tau}t_d$$

Since we know our pulse travels at the speed of light $$c$$, we can rewrite the target's delay in terms of range $$R$$ and yield the mapping between target range and its beat frequency

$$t_d = \frac{2R}{c} => f_b = -\frac{2R\beta}{c\tau}$$

So generating a dechirped signal is straight forward since it's simply a sinusoid at some beat frequency $$f_b$$.

Note that these equations apply only to upchirps and downchirps. The negative sign takes care of itself in either case. Triangular and more exotic chirps yield additional frequency terms but this process can be expanded to cover that as well.

To add Doppler, you can add a constant phase term that is updated as you collect pulses to form a range-Doppler map. You can actually start at zero phase for the first pulse and progress from there for the purpose of a simulation. Your additional phase term would look something like

$$e^{j2{\pi}f_dnT_c}$$

Where $$n$$ is the current pulse number starting at 0 and $$T_c$$ is your equivalent pulse repetition interval (PRI) which is akin to your sweep time for FMCW radars.

EDIT: After having some time to look at your code directly, I found a few issues.

First, you are missing a factor of two in the Doppler component of the phase.

Second, without getting into theory, your particular system supports a Doppler span wider than what you force the horizontal axis to. This will erroneously change where you perceive the target to be.

Third, the time vector you are using to generate the beat frequency must be reset to 0 after every pulse. This is because the time vector needs to be relative to the time delay of the target $$t_d$$.

Here is your modified code. I currently do not have the Phased Array toolbox to generate and display the range-Doppler map, so I did it manually.

c = 3e8;            %speed of light
range_max = 180;    %max detection range
tm = 6*(2*range_max/c); %sweep time
%tm is 7.2e-6 s
bw = 200e6;         %sweep bandwidth
sweep_slope = bw/tm;
v_max = 150*1000/3600;  %target max velocity
fc = 77e9;          %radar frequency
lambda = c/fc;      %radar wavelength
fs = 72e6;          %sampling rate
%sampling rate based on ADC datasheet
chirps = 64;        %frame size
samples = ceil(tm*fs);  %samples in one chirp

%% target
R0 = 20; %range in meters
V = 40; %radial velocity, m/s

%%
t = 0; %time
mix = zeros(samples, chirps); %mixer output
for i=1:1:chirps
td = 2 * R0 / c; %round trip delay
phi0 = 4*pi*fc*R0/c; %inital phase
t = 0; % Reset

for j=1:1:samples
a = (-2*pi*fc*2*V*i*tm/c ...    %phase shift
-2*pi*(2*V*(fc+i*bw)/c + sweep_slope*td)*t);   %frequency
mix(j,i) = 0.5*cos(a);
t = t + 1/fs;
end
end

%% Form the range-Doppler map (RDM)

% RDM axes
rangeBinAxis = (0:samples-1).*c/(2*bw);
dopplerBinSize = (1/tm)/chirps;
velocityBinAxis = (-chirps/2:chirps/2-1).*dopplerBinSize*lambda/2;

% 2D FFT to perform range and Doppler compression (i.e. form the RDM)
rdm = fftshift(fft2(mix), 2);

% Plot the RDM for the valid ranges of interest - targets ahead of you
figure;
surf(velocityBinAxis, rangeBinAxis(1:ceil(samples/2)), 20*log10(abs(rdm(1:ceil(samples/2), :))));
% surf(velocityBinAxis, rangeBinAxis, 20*log10(abs(rdm)));  % See the entire spectrum
xlabel("Range (m)");
ylabel("Velocity (m/s)");
axis tight;
shading flat;
view(0, 90);
colorbar;

% figure(1)
% rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,...
%     'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,...
%     'RangeMethod','FFT','SweepSlope',sweep_slope,...
%     'RangeFFTLengthSource','Property','RangeFFTLength',2048,...
%     'DopplerFFTLengthSource','Property','DopplerFFTLength',256);
%
% clf;
% plotResponse(rngdopresp,mix);
% axis([-v_max v_max 0 range_max])


Some examples to show the target being positioned appropriately in the RDM (scale is in dB):

R0 = 32 m , v = 40 m/s

R0 = 150 m , v = 40 m/s

• So, seeing that the equation that I used in my simulation was partially right, I suppose I was on the right track with the form of the signal. But after leaving only the range term and the constant phase shift term and starting from zero phase, I yet again see a target at the right range but with a completely wrong velocity on my range-Doppler map in matlab. I now think that my struggle has more to do with the specific implementation in matlab. Though it bothers me that changing range by a small amount drastically affects the velocity on the map. Jun 16 '20 at 1:44
• @narcoleptic Updated the post. Note that I used a 2D FFT of the same dimensions as the collected data. You can modify the code so that it takes larger FFTs as desired. This should clear up your issues! Jun 16 '20 at 22:39
• The need to reset the time for each chirp is an eye-opener, it would have taken me ages to realize that. I suppose it renders the time variable mostly useless, since I can just use the (sample iterator - 1) / fs instead. Also, manually plotting the map was another thing I was interested in, since I couldn't get my contourf of fft2 to look very presentable, but I didn't want to clutter my post with too many questions. So you helped me with even more stuff, than I could have hoped for. Thank you! Jun 17 '20 at 2:55
• @Envidia .. why do RD maps have end to end streaks along the range,velcoty calculated values ? In the first figure,, why not just a point(32,40) ... and what about the image streaks ? how can we remove them just to get only one point ? Oct 29 '21 at 14:14
• @BandW Because the RDM is created using a finite length DFT in both dimensions, it creates sidelobes which is what you are seeing. You will never get just a point in the response, since that requires that you observe and sample the signal over infinite time, which is not practically possible. By the way, this example did not have any noise, so the sidelobes are more apparent. Nov 1 '21 at 5:09