Generating IF signal for FMCW radar

Question

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 :)

Answer 1

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