Joint Radar-Communication Using PMCW and OFDM Waveforms

This example shows how to model a joint radar-communication (JRC) system using the Phase Array System Toolbox. Two modulation schemes for combining radar and communication signals are presented. The first scheme uses the phase-modulated continuous wave (PMCW) waveform. The second scheme uses the orthogonal frequency division multiplexing (OFDM) waveform. For each modulation scheme the example shows how to process the received signal by both the radar and the communication receivers.

JRC System Parameters

JRC is one of many different radar and communication coexistence strategies. Such strategies have emerged recently as a response to an increasingly congested frequency spectrum and a demand of radar and communication systems for a large bandwidth. In a JRC system both the radar and the communication share the same platform and have a common transmit waveform. This example shows two approaches to utilizing a single waveform for both functions. The first approach uses PMCW. Such a JRC system can be considered as radar-centric, since PMCW is a waveform tailored to radar applications. The second approach uses OFDM. This approach is communication-centric since a communication waveform is used to perform the radar function.

Consider a JRC system that operates at a carrier frequency of 24 GHz and has a bandwidth of 100 MHz.

% Set the random number generator for reproducibility
rng('default');

fc = 24e9;                              % Carrier frequency (Hz)
B = 100e6;                              % Bandwidth (Hz)

Set the peak transmit power to 10 mW and the transmit antenna gain to 20 dB.

Pt = 0.01;                              % Peak power (W)
Gtx = 20;                               % Tx antenna gain (dB)

Let the JRC transmitter and the radar receiver be colocated at the origin. Create a phased.Platform object to represent the JRC system.

jrcpos = [0; 0; 0];                     % JRC position
jrcvel = [0; 0; 0];                     % JRC velocity

jrcmotion = phased.Platform('InitialPosition', jrcpos, 'Velocity', jrcvel);

Assume the radar receiver has the following parameters.

Grx = 20;                               % Radar Rx antenna gain (dB)
NF = 2.9;                               % Noise figure (dB)
Tref = 290;                             % Reference temperature (K)

Let the maximum range of interest and the maximum observed relative velocity be 200 m and 60 m/s respectively.

Rmax = 200;                             % Maximum range of interest
vrelmax = 60;                           % Maximum relative velocity

Consider a scenario with three moving radar targets and a single downlink user. Create a phased.Platform object to represent moving targets and a phased.RadarTarget object to represent the target's radar cross sections (RCS).

tgtpos = [110 45 80; -10 5 -4; 0 0 0];  % Target positions
tgtvel = [-15 40 -32; 0 0 0; 0 0 0];    % Target velocities
tgtmotion = phased.Platform('InitialPosition', tgtpos, 'Velocity', tgtvel);

tgtrcs = [4.7 3.1 2.3];                 % Target RCS
target = phased.RadarTarget('Model', 'Swerling1', 'MeanRCS', tgtrcs, 'OperatingFrequency', fc);

Create a position vector to represent the position of the downlink user.

userpos = [100; 20; 0];                  % Downlink user position                   

Visualize the JRC scenario using helperPlotJRCScenario helper function.

helperPlotJRCScenario(jrcpos, tgtpos, userpos, tgtvel);

JRC System Using a PMCW Waveform

The first modulation scheme considered in this example uses a PMCW waveform. A PMCW radar repeatedly transmits a selected phase-coded sequence. The duration of the transmitted sequence is called a modulation period. Since the duty cycle of a PMCW radar is equal to one, the next modulation period starts right after the previous.

This example uses a maximum length pseudo-random binary sequence (PRBS) also frequently referred to as an m-sequence. Maximum length sequences have low autocorrelation sidelobes that result in a good range estimation performance. The communication data is modulated on top of the radar waveform such that each modulation period carries a single bit. This is achieved by simply multiplying the entire modulation period by the corresponding BPSK symbol.

Using the helperMLS helper function create a PRBS.

% Generate an m-sequence. The sequence lengths has to be 2^p-1, where p is
% an integer.
p = 8;
prbs = helperMLS(p);
Nprbs = numel(prbs);                    % Number of chips in PRBS

Given that the chip duration of a phase-modulated waveform is an inverse of the total bandwidth, compute the modulation period.

Tchip = 1/B;                             % Chip duration
Tprbs = Nprbs * Tchip;                   % Modulation period

Generate the communication data to be transmitted. Encode the generated data bits using the BPSK modulation.

Mpmcw = 256;                             % Number of transmitted data bits 
dataTx = randi([0, 1], [Mpmcw 1]);       % Binary data
bpskTx = pskmod(dataTx, 2);              % BPSK symbols

Since the communication channel is unknown, the JRC system must perform channel sounding. One way to do it is to send out a single period of the PRBS sequence unmodulated with the communication data. Thus, a transmission of one data bit would require transmitting the PRBS sequence twice. The first time to sound the channel, the second time to transmit the data. Create the PMCW waveform by stacking two modulation periods together, the first period is just the generated PRBS and the second period is the generated PRBS multiplied by the corresponding BPSK symbol.

% Transmit waveform
xpmcw = [prbs * ones(1, Mpmcw); prbs * bpskTx.'];

Thus, the duration of one PMCW block is twice the modulation period.

Tpmcw = 2*Tprbs; 

The block duration can be longer if the same data symbol has to be repeated multiple times in order to improve the signal-to-noise ratio (SNR) by integrating multiple periods at the receiver.

Assume that the coherent processing interval is long enough to coherently process Mpmcw blocks. Compute the duration of a single PMCW frame that carries Mpmcw data bits.

Tpmcw*Mpmcw

ans = 
0.0013

Radar Signal Simulation and Processing

This section of the example shows how to simulate the radar sensing part of the JRC system with the PMCW waveform. First, generate received PMCW blocks.

% Set the sample rate equal to the bandwidth
fs = B;

% Setup the transmitter and the radiator
transmitter = phased.Transmitter('Gain', Gtx, 'PeakPower', Pt);

% Assume the JRC is using an isotropic antenna
ant = phased.IsotropicAntennaElement;
radiator = phased.Radiator('Sensor', ant, 'OperatingFrequency', fc);

% Setup the collector and the receiver
collector = phased.Collector('Sensor', ant, 'OperatingFrequency', fc);
receiver = phased.ReceiverPreamp('SampleRate', fs, 'Gain', Grx, 'NoiseFigure', NF, 'ReferenceTemperature', Tref);

% Setup a free space channel to model the JRC signal propagation from the
% transmitter to the targets and back to the radar receiver
radarChannel = phased.FreeSpace('SampleRate', fs, 'TwoWayPropagation', true, 'OperatingFrequency', fc);

% Preallocate space for the signal received by the radar receiver
ypmcwr = zeros(size(xpmcw));

% Transmit one PMCW block at a time assuming that all Mpmcw blocks are
% transmitted within a single CPI
for m = 1:Mpmcw
    % Update sensor and target positions
    [jrcpos, jrcvel] = jrcmotion(Tpmcw);
    [tgtpos, tgtvel] = tgtmotion(Tpmcw);
    
    % Calculate the target angles as seen from the transmit array
    [tgtrng, tgtang] = rangeangle(tgtpos, jrcpos);

    % Transmit signal
    txsig = transmitter(xpmcw(:, m));

    % Radiate signal towards the targets
    radtxsig = radiator(txsig, tgtang);

    % Apply free space channel propagation effects
    chansig = radarChannel(radtxsig, jrcpos, tgtpos, jrcvel, tgtvel); 

    % Reflect signal off the targets
    tgtsig = target(chansig, false);
    
    % Receive target returns at the receive array
    rxsig = collector(tgtsig, tgtang);
 
    % Add thermal noise at the receiver
    ypmcwr(:, m) = receiver(rxsig);
end

Now the received signals can be processed to obtain estimates of the targets ranges and velocities. The following diagram shows the signal processing steps needed to compute the range-Doppler map.

Start by removing the communication data. Multiply the modulation periods that carry the data by the corresponding BPSK symbols.

ypmcwr1 = ypmcwr(Nprbs+1:end, :) .* (bpskTx.');

The modulation period used for the communication channel sounding can be integrated with the modulation period that carried the data.

ypmcwr1 = ypmcwr1 + ypmcwr(1:Nprbs, :);

It is more computationally efficient to perform matched filtering in the frequency domain. First, compute the frequency domain representation of the integrated signal then, multiply the result by the complex-conjugate version of the DFT of the used PRBS.

% Frequency-domain matched filtering
Ypmcwr = fft(ypmcwr1);
Sprbs = fft(prbs);
Zpmcwr = Ypmcwr .* conj(Sprbs);

Use phased.RangeDopplerResponse to compute and plot the range-Doppler map.

% Range-Doppler response object computes DFT in the slow-time domain and
% then IDFT in the fast-time domain to obtain the range-Doppler map
rdr = phased.RangeDopplerResponse('RangeMethod', 'FFT', 'SampleRate', fs, 'SweepSlope', -B/Tprbs,...
    'DopplerOutput', 'Speed', 'OperatingFrequency', fc, 'PRFSource', 'Property', 'PRF', 1/Tpmcw,...
    'ReferenceRangeCentered', false);

figure;
plotResponse(rdr, Zpmcwr, 'Unit', 'db');
xlim([-vrelmax vrelmax]);
ylim([0 Rmax]);

Communication Signal Simulation and Processing

This example uses a stochastic Rician channel model to simulate the propagation of the JRC signal from the transmitter to the downlink receiver. Create a communication channel object and pass the PMCW signal through the channel.

% Line-of-sight between the JRC transmitter and the downlink user
dlos = vecnorm(jrcpos - userpos);
taulos = dlos/physconst('LightSpeed');

% Delays and gains due to different path in the communication channel
pathDelays = taulos * [1 1.6 1.1 1.45] - taulos;
averagePathGains = [0 -1 2 -1.5];

% Maximum two-way Doppler shift
fdmax = 2 * speed2dop(vrelmax, freq2wavelen(fc));
commChannel = comm.RicianChannel('PathGainsOutputPort', true, 'DirectPathDopplerShift', 0,... 
    'MaximumDopplerShift', fdmax/2, 'PathDelays', pathDelays, 'AveragePathGains', averagePathGains, ...
    'SampleRate', fs);

% Pass the transmitted signal through the channel
ypmcwc = commChannel(xpmcw(:));

For simplicity, model the downlink receiver using additive gaussian noise and assuming that the signal-to-noise ratio (SNR) is 40 dB.

SNRu = 40;                               % SNR at the downlink user's receiver 
ypmcwc = awgn(ypmcwc, SNRu, "measured");
ypmcwc = reshape(ypmcwc, 2*Nprbs, []);

Use the first PRBS repetition in each PMCW block to estimate the channel frequency response.

ysound = ypmcwc(1:Nprbs, :); 
Hpmcw = fft(ysound)./Sprbs;

Once the channel estimates are available, the received signal can be processed to extract the user data. The diagram below show the required processing steps.

Since a radar waveform is used to carry the user data, the matched filtering is the first step in the demodulation process. Perform matched filtering in the frequency domain.

% Compute the frequency domain representation of the received signal
Ypmcwc = fft(ypmcwc(Nprbs + 1:end, :));

% Multiply by the complex-conjugate version of the DFT of the used PRBS
Zpmcwc = Ypmcwc .* conj(Sprbs);

Multiply the frequency representation of the matched filtered signal by the inverse of the estimated channel frequency response to remove the channel effects.

Zpmcwc = Zpmcwc./Hpmcw;

Demodulate the signal to first obtain the BPSK symbols and then the binary data.

zpmcwc = ifft(Zpmcwc);

[~, idx] = max(abs(zpmcwc), [], 'linear');
bpskRx = zpmcwc(idx).';
bpskRx = bpskRx./Nprbs;
dataRx = pskdemod(bpskRx, 2);

refconst = pskmod([0 1], 2);
constellationDiagram = comm.ConstellationDiagram('NumInputPorts', 1, ...
    'ReferenceConstellation', refconst, 'ChannelNames', {'Received PSK Symbols'});
constellationDiagram(bpskRx(:));

Compute the bit error rate.

[numErr, ratio] = biterr(dataTx, dataRx)

numErr = 
0

ratio = 
0

JRC System Using an OFDM Waveform

OFDM is a second modulation scheme considered in this example. OFDM signals are widely used for wireless communication. OFDM is a multicarrier signal composed of a set of orthogonal complex exponentials known as subcarriers. The complex amplitude of each subcarrier can be used to carry communication data allowing for a high data rate. OFDM waveforms were also noted to exhibit Doppler tolerance and a lack of range-Doppler coupling. These OFDM features are attractive to radar applications. This section of the example shows how to use an OFDM waveform to perform both the radar sensing and communication.

Let the number of OFDM subcarriers be 1024. Compute the subcarrier separation given that all subcarriers must be equally spaced in frequency.

Nsc = 1024;                              % Number of subcarriers
df = B/Nsc;                              % Separation between OFDM subcarriers
Tsym = 1/df;                             % OFDM symbol duration

According to the OFDM system design rules, the subcarrier spacing must be at least ten times larger than the maximum Doppler shift experienced by the signal. This is required to ensure the orthogonality of the subcarrier. Verify that this condition is met for the scenario considered in this example.

df > 10*fdmax

ans = logical
   1

In order to avoid intersymbol interference, the OFDM symbols are prepended with cyclic prefix. Compute the cyclic prefix duration based on the maximum range of interest.

Tcp = range2time(Rmax);                  % Duration of the cyclic prefix (CP)
Ncp = ceil(fs*Tcp);                      % Length of the CP in samples
Tcp = Ncp/fs;                            % Adjust duration of the CP to have an integer number of samples

Tofdm = Tsym + Tcp;                      % OFDM symbol duration with CP
Nofdm = Nsc + Ncp;                       % Number of samples in one OFDM symbol

Designate some of the subcarriers to act as guard bands.

nullIdx = [1:9 (Nsc/2+1) (Nsc-8:Nsc)]';  % Guard bands and DC subcarrier
Nscd = Nsc-length(nullIdx);              % Number of data subcarriers

Generate binary data to be transmitted. Since each OFDM subcarrier can carry a complex data symbol, the generated binary data first must be modulated. Use QAM-modulation to modulate the generated binary data. Finally, OFDM-modulate the QAM-modulated symbols.

bps = 6;                                 % Bits per QAM symbol (and OFDM data subcarrier)
K = 2^bps;                               % Modulation order
Mofdm = 128;                             % Number of transmitted OFDM symbols
dataTx = randi([0,1], [Nscd*bps Mofdm]);
qamTx = qammod(dataTx, K, 'InputType', 'bit', 'UnitAveragePower', true);
ofdmTx = ofdmmod(qamTx, Nsc, Ncp, nullIdx);

Assume that the coherent processing interval is long enough to coherently process Mofdm symbols. Compute the duration of a single OFDM frame and a total number of data bits in one frame.

Tofdm*Mofdm                              % Frame duration (s)

ans = 
0.0015

Nscd*bps*Mofdm                           % Total number of bits

ans = 
771840

Notice that the OFDM frame duration is close to the PMCW frame duration, however the total number of transmitted bits is larger by several orders of magnitude.

Radar Signal Simulation and Processing

Use the same system as in the PMCW case to simulate transmission of the OFDM signal by the JRC and the reception of the target returns by the radar receiver.

xofdm = reshape(ofdmTx, Nofdm, Mofdm);

% OFDM waveform is not a constant modulus waveform. The generated OFDM
% samples have power much less than one. To fully utilize the available
% transmit power, normalize the waveform such that the sample with the
% largest power had power of one.
xofdm = xofdm/max(sqrt(abs(xofdm).^2), [], 'all');

% Preallocate space for the signal received by the radar receiver
yofdmr = zeros(size(xofdm));

% Reset the platform objects representing the JRC and the targets before
% running the simulation
reset(jrcmotion);
reset(tgtmotion);

% Transmit one OFDM symbol at a time assuming that all Mofdm symbols are
% transmitted within a single CPI
for m = 1:Mofdm
    % Update sensor and target positions
    [jrcpos, jrcvel] = jrcmotion(Tofdm);
    [tgtpos, tgtvel] = tgtmotion(Tofdm);
    
    % Calculate the target angles as seen from the transmit array
    [tgtrng, tgtang] = rangeangle(tgtpos, jrcpos);

    % Transmit signal
    txsig = transmitter(xofdm(:, m));

    % Radiate signal towards the targets
    radtxsig = radiator(txsig, tgtang);

    % Apply free space channel propagation effects
    chansig = radarChannel(radtxsig, jrcpos, tgtpos, jrcvel, tgtvel); 

    % Reflect signal off the targets
    tgtsig = target(chansig, false);
    
    % Receive target returns at the receive array
    rxsig = collector(tgtsig, tgtang);
 
    % Add thermal noise at the receiver
    yofdmr(:, m) = receiver(rxsig);
end

Process the target returns to obtain a range-Doppler map. The signal processing steps are shown in the block diagram.

Remove the cyclic prefix and compute the DFT in the fast-time domain using the ofdmdemod function.

yofdmr1 = reshape(yofdmr, Nofdm*Mofdm, 1);

% Demodulate the received OFDM signal (compute DFT)
Yofdmr = ofdmdemod(yofdmr1, Nsc, Ncp, Ncp, nullIdx);

The complex amplitudes of the subcarriers that comprise the received OFDM signal are the transmitted QAM symbols. Therefore, to remove the communication data divide the DFT of the received signal by the transmitted QAM symbols.

Zofdmr = Yofdmr./qamTx;

Use phased.RangeDopplerResponse to compute and plot the range-Doppler map.

% Range-Doppler response object computes DFT in the slow-time domain and
% then IDFT in the fast-time domain to obtain the range-Doppler map
rdr = phased.RangeDopplerResponse('RangeMethod', 'FFT', 'SampleRate', fs, 'SweepSlope', -B/Tofdm,...
    'DopplerOutput', 'Speed', 'OperatingFrequency', fc, 'PRFSource', 'Property', 'PRF', 1/Tofdm, ...
    'ReferenceRangeCentered', false);

figure;
plotResponse(rdr, Zofdmr, 'Unit', 'db');
xlim([-vrelmax vrelmax]);
ylim([0 Rmax]);

Communication Signal Simulation and Processing

Use the same stochastic Rician channel model as in the case of the PMCW to model the signal propagation from the JRC transmitter to the downlink user.

[yofdmc, pathGains] = commChannel(ofdmTx);

Add the receiver noise.

yofdmc = awgn(yofdmc, SNRu, "measured");

Compute the OFDM channel response to perform channel equalization.

channelInfo = info(commChannel);
pathFilters = channelInfo.ChannelFilterCoefficients;
toffset = channelInfo.ChannelFilterDelay;
Hofdm = ofdmChannelResponse(pathGains, pathFilters, Nsc, Ncp, ...
    setdiff(1:Nsc, nullIdx), toffset);

Perform demodulation and equalization. The block diagram of the processing steps is shown below.

zeropadding = zeros(toffset, 1);
ofdmRx = [yofdmc(toffset+1:end,:); zeropadding];

qamRx = ofdmdemod(ofdmRx, Nsc, Ncp, Ncp/2, nullIdx);
qamEq = ofdmEqualize(qamRx, Hofdm(:), 'Algorithm', 'zf');

Demodulate the QAM symbols to obtain the received binary data.

dataRx = qamdemod(qamEq, K, 'OutputType', 'bit', 'UnitAveragePower', true);

refconst = qammod(0:K-1, K, 'UnitAveragePower', true);
constellationDiagram = comm.ConstellationDiagram('NumInputPorts', 1, ...
    'ReferenceConstellation', refconst, 'ChannelNames', {'Received QAM Symbols'});

% Show QAM data symbols comprising the first 10 received OFDM symbols
constellationDiagram(qamEq(1:Nscd*10).');

Compute the bit error rate.

[numErr,ratio] = biterr(dataTx, dataRx)

numErr = 
14754

ratio = 
0.0191

Conclusions

This example simulates a JRC system using two different modulation schemes for generating a transmit waveform: PMCW and OFDM. The transmit waveforms are generated such that they can be used to sense radar targets and to carry communication data. The example shows how to simulate transmission, propagation, and reception of these waveforms by the radar receiver and by the downlink user. For each modulation scheme the example shows how to perform radar signal processing to obtain a range-Doppler map of the observed scenario. The example also shows how to process the received signal at the downlink receiver to obtain the transmitted communication data.

References

1. Sturm, Christian, and Werner Wiesbeck. "Waveform design and signal processing aspects for fusion of wireless communications and radar sensing." Proceedings of the IEEE 99, no. 7 (2011): 1236-1259.

2. de Oliveira, Lucas Giroto, Benjamin Nuss, Mohamad Basim Alabd, Axel Diewald, Mario Pauli, and Thomas Zwick. "Joint radar-communication systems: Modulation schemes and system design." IEEE Transactions on Microwave Theory and Techniques 70, no. 3 (2021): 1521-1551.

Supporting Functions

function seq = helperMLS(p)
% Generate a pseudorandom sequence of length N=2^p-1, where p is an
% integer. The sequence is generated using shift registers. The feedback
% coefficients for the registers are obtained from the coefficients of an
% irreducible, primitive polynomial in GF(2p).

    pol = gfprimdf(p, 2);
    seq = zeros(2^p - 1, 1);
    seq(1:p) = randi([0 1], p, 1);
    
    for i = (p + 1):(2^p - 1)
        seq(i) = mod(-pol(1:p)*seq(i-p : i-1), 2);
    end
    
    seq(seq == 0) = -1;
end