4
$\begingroup$

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 IF signal formula 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:
range-doppler 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 :)

$\endgroup$
2
  • $\begingroup$ Why the sweep time parameter define as 6*(2*range_max/c); $\endgroup$ Commented Mar 22, 2021 at 12:02
  • $\begingroup$ @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 $\endgroup$ Commented Mar 25, 2021 at 4:32

1 Answer 1

4
$\begingroup$

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 enter image description here

R0 = 150 m , v = 40 m/s enter image description here

$\endgroup$
11
  • $\begingroup$ 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. $\endgroup$ Commented Jun 16, 2020 at 1:44
  • 1
    $\begingroup$ 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! $\endgroup$ Commented Jun 17, 2020 at 2:55
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Nov 1, 2021 at 5:09
  • 1
    $\begingroup$ @BandW Adding noise does not make it better. It just seems better because you don't visually see the sidelobes, but that's because the noise is covering them up. $\endgroup$ Commented Nov 3, 2021 at 4:27
  • 1
    $\begingroup$ @Blupon The beat frequency can be negative depending on what conventions are being used to define the waveform and how the mixing is done. Even if the frequency is placed on the negative side of the spectrum, that's still ok as the magnitude value still maps correctly to range. Regarding RVP, SAR techniques are a two-dimensional operation, where the cross-range resolution is achieved by exploiting the changing range between the target and the platform. RVP changes pulse-to-pulse as the platform flies, and this phase error corrupts the imagery if not accounted for. $\endgroup$ Commented Apr 16 at 14:40

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.