Cooperative Bistatic Radar I/Q Simulation and Processing

Since R2024b

This example uses:

A cooperative bistatic radar is a type of bistatic radar system in which the transmitter and receiver are deliberately designed to collaborate, despite not being collocated. Unlike passive bistatic radar systems, which utilize transmissions from unrelated sources like radio and television broadcasts, cooperative bistatic radars involve transmitters that are specifically engineered to work in conjunction with one or more receivers for the purpose of detecting and tracking targets.

This example demonstrates the simulation of I/Q for a cooperative bistatic system and its subsequent signal processing for a single bistatic transmitter and receiver pair. The cooperative bistatic radar system assumes the:

Transmitter and receiver are time synchronized
Transmitter and receiver positions are known
Transmitted waveform is known
Target is a perfect electrical conductor (PEC) with no motion (e.g., rotation or tumbling)

You will learn in this example how to:

Setup a bistatic scenario
Calculate and import a bistatic radar cross section (RCS) for a cone target using Antenna Toolbox™
Simulate a bistatic datacube
Perform range and Doppler processing
Form surveillance beams and null the direct path
Perform detection and parameter estimation

This example uses Radar Toolbox™ and Antenna Toolbox™.

Set Up the Bistatic Scenario

First initialize the random number generation for reproducible results and define radar parameters for a 300 MHz cooperative bistatic radar system. The maximum range is set to 48 km, and the range resolution is 50 meters.

% Set random number generation for reproducibility
rng('default') 

% Radar parameters
fc         = 3e8;                 % Operating frequency (Hz)
maxrng     = 48e3;                % Maximum range (m)
rngres     = 50;                  % Range resolution (m)
tbprod     = 20;                  % Time-bandwidth product
[lambda,c] = freq2wavelen(fc);    % Wavelength (m)

Create a radar scenario.

% Set up scenario
pri = range2time(maxrng);
scene = radarScenario(UpdateRate=1/pri);

Define the platforms in the scenario. This example includes a transmitter, receiver, and a single target. The stationary transmitter and receiver are separated by a distance of 10 km. The target moves with a velocity of 150 m/s in the y-direction away from the transmitter and receiver.

% Define platforms
txPlat  = platform(scene,Position=[-5e3 0 0],Orientation=rotz(45).');
rxPlat  = platform(scene,Position=[5e3 0 0],Orientation=rotz(135).');
tgtspd = -200;
traj = kinematicTrajectory(Position=8e3*[cosd(60) sind(60) 0],Velocity=[0 150 0]);
tgtPlat = platform(scene,Trajectory=traj);

Visualize the scenario using theaterPlot in the helper helperPlotScenario. The transmitter is shown below as an upward-facing blue triangle. The receiver is shown as a downward-facing blue triangle. The red, green, and blue lines form axes that indicate the orientations of the transmitter and receiver. The target is shown as a red circle with its trajectory shown in gray.

% Visualize the scenario 
helperPlotScenario(txPlat,rxPlat,tgtPlat)

Configure Bistatic Transmitter and Receiver

Now define the bistatic transmitter and receiver.

The transmitter transmits a linear frequency modulated (LFM) waveform with parameters calculated from the previously defined maximum range and range resolution.

% Waveform
pri = range2time(maxrng);
bw  = rangeres2bw(rngres);
fs  = 2*bw;
waveform = phased.LinearFMWaveform(PulseWidth=tbprod/bw,PRF=1/pri,...
    SweepBandwidth=bw,SampleRate=fs);

The peak power of the transmitter is set to 500 Watts with a loss of 3 dB. Define the transmit antenna as a uniform linear array (ULA) with 4 vertically polarized short dipole antenna elements spaced a half-wavelength apart. Attach the array to a phased.Radiator object. This example models a polarimetric scene, so set the Polarization property of the phased.Radiator to 'Combined'.

% Transmitter
transmitter = phased.Transmitter(PeakPower=500,Gain=0,LossFactor=3);
txAntenna   = phased.ShortDipoleAntennaElement(AxisDirection='Z');
txarray     = phased.ULA(4,lambda/2,Element=txAntenna);
radiator    = phased.Radiator(Sensor=txarray,PropagationSpeed=c,...
    Polarization='Combined');

Define the receive antenna as a uniform linear array (ULA) with 16 vertically polarized short dipole antenna elements spaced a half-wavelength apart. Attach the array to a phased.Collector object. Again, set the Polarization property of the phased.Collector to 'Combined'. The receiver has a noise figure of 6 dB.

% Receiver
rxAntenna = phased.ShortDipoleAntennaElement(AxisDirection='Z');
rxarray   = phased.ULA(16,lambda/2,Element=rxAntenna);
collector = phased.Collector(Sensor=rxarray,...
    PropagationSpeed=c,...
    OperatingFrequency=fc,...
    Polarization='Combined');
receiver  = phased.Receiver(Gain=0,SampleRate=fs,NoiseFigure=6);

Define the Cone Target

Define the target structure as a cylinder with radius of 0.05 m at one end and 0.5 meters at the other end. This forms a cone. Set the height of the cone to 2 meters.

% Define cone target
fc         = collector.OperatingFrequency;
cyl        = shape.Cylinder;
cyl.Radius = [0.05 0.5]; 
cyl.Height = 2; 
tgtCone = customAntenna(Shape=cyl,TiltAxis='X',Tilt=90);

Mesh the cone figure with a maximum edge length equal to one-eighth the wavelength.

% Mesh and plot
figure
mesh(tgtCone,MaxEdgeLength=lambda/8);

Figure contains an axes object and an object of type uicontrol. The axes object with title Metal mesh, xlabel x (m), ylabel y (m) contains an object of type patch. This object represents PEC.

stlwrite(tgtCone,'tgtCone.stl');
p = platform;
p.FileName = 'tgtCone.stl';
p.Units = 'm';
pause(0.25)

Initialize a radar target container for the bistatic RCS.

% Target RCS
target = phased.RadarTarget(EnablePolarization=true,Mode='Bistatic',ScatteringMatrixSource='Input port');

Visualize the target in the scenario.

% Visualize target structure in theater plot
hTgtConeAxes = gca; 
helperPlotTarget(hTgtConeAxes,tgtPlat)

Figure contains an axes object. The axes object with title Target in Scenario, xlabel X (m), ylabel Y (m) contains 2 objects of type line, patch. One or more of the lines displays its values using only markers This object represents Target.

Define Bistatic Propagation Paths

Define the propagation paths as free space channels.

There are 3 paths to model in this example:

Transmitter-to-receiver: This is also referred to as the direct path or baseline $L$ .
Transmitter-to-target: This is the first leg of the bistatic path and is denoted as $R_{T}$ in the diagram above.
Target-to-receiver: This is the second leg of the bistatic path and is denoted as $R_{R}$ in the diagram above.

The bistatic range is relative to the direct path and is given by

$R_{Bi} = R_{T} + R_{R} - L$ .

% Direct path 
dpchannel = phased.FreeSpace(SampleRate=fs,...
    OperatingFrequency=fc,PropagationSpeed=c);
smat = sqrt(lambda^2./(4*pi)).*ones(2); % Model direct path return as a target
dp = phased.RadarTarget(EnablePolarization=true,Mode='Bistatic',ScatteringMatrix=smat);

% Transmitter-to-target path
txchannel = phased.FreeSpace(SampleRate=fs,...
    OperatingFrequency=fc,PropagationSpeed=c);

% Target-to-receiver path
rxchannel = phased.FreeSpace(SampleRate=fs,...
    OperatingFrequency=fc,PropagationSpeed=c);

Simulate I/Q Data

First initialize the datacube.

% Initialize datacube 
numPulses   = 65; % Number of pulses to record 
numSamples  = round(fs*pri); % Number of range samples
numElements = collector.Sensor.NumElements; % Number of elements in receive array
y           = zeros(numSamples*numPulses,numElements,like=1i); % Datacube 
endIdx      = 0; % Index into datacube

Advance the scene and calculate the positions of the transmitter, receiver, and target. Assume that the target is stationary within the duration of a coherent processing interval (CPI) for the purposes of calculating the bistatic RCS. If running in a loop, update the RCS for each CPI.

% Advance scene
advance(scene);

% Update positions of transmitter, receiver, and targets
[txpos,txvel,txax] = helperUpdatePosition(scene,txPlat);
[rxpos,rxvel,rxax] = helperUpdatePosition(scene,rxPlat);
[tgtpos,tgtvel,tgtax] = helperUpdatePosition(scene,tgtPlat);

Calculate the scattering matrix of the cone target for the current CPI. The scattering matrix is given by

$S = [\begin{array}{c} S_{HH} & S_{HV} \\ S_{VH} & S_{VV} \end{array}]$ .

Each component of the scattering matrix is calculated separately below.

% Calculate bistatic RCS
[Rtx,fwang] = rangeangle(txpos,tgtpos,tgtax);
[Rrx,bckang] = rangeangle(rxpos,tgtpos,tgtax);
smat(1,1) = rcs(p,fc,TransmitAngle=fwang,ReceiveAngle=bckang,Type="Complex",Polarization="HH",Scale="linear");
smat(2,2) = rcs(p,fc,TransmitAngle=fwang,ReceiveAngle=bckang,Type="Complex",Polarization="VV",Scale="linear");
smat(1,2) = rcs(p,fc,TransmitAngle=fwang,ReceiveAngle=bckang,Type="Complex",Polarization="HV",Scale="linear");
smat(2,1) = rcs(p,fc,TransmitAngle=fwang,ReceiveAngle=bckang,Type="Complex",Polarization="VH",Scale="linear");

Now simulate receiving 65 pulses.

% Simulate CPI 
for ip = 1:numPulses

First transmit the waveform and propagate through the direct path channel $L$ .

    % Direct path
    wav = waveform();
    wav = transmitter(wav);
    [~,dpang] = rangeangle(rxpos,txpos,txax);
    sigdp = radiator(wav,dpang,txax);
    sigdp = dp(sigdp,dpang,dpang,tgtax);
    sigdp = dpchannel(sigdp,txpos,rxpos,txvel,rxvel);

Next, calculate the transmitter-to-target path $R_{T}$ .

    % Transmitter-to-target path
    [~,radang] = rangeangle(tgtpos,txpos,txax); % Target angles as seen by the transmitter
    sigtx = radiator(wav,radang,txax);
    sigtx = txchannel(sigtx,txpos,tgtpos,txvel,tgtvel);

Reflect the signal off the target.

    % Reflect off target
    sigtgt = target(sigtx,fwang,bckang,tgtax,smat);

Then propagate the target signal through the target-to-receive path $R_{R}$ .

    % Target-to-receive path
    sigrx = rxchannel(sigtgt,tgtpos,rxpos,tgtvel,rxvel);

Collect both the target and direct path signals at the collector.

    % Collect direct path and target signals at the collector 
    [~,inangTgt] = rangeangle(tgtpos,rxpos,rxax);
    [dpRng,inangDp] = rangeangle(txpos,rxpos,rxax);    
    sigrx = collector([sigrx sigdp],[inangTgt inangDp],rxax);

Save off the received data and advance the scene.

    % Add data to datacube
    startIdx = endIdx + 1;
    endIdx = startIdx + numSamples - 1;
    y(startIdx:endIdx,:) = receiver(sigrx);

    % Advance scene
    advance(scene);

    % Update positions of transmitter, receiver, and targets
    [txpos,txvel,txax] = helperUpdatePosition(scene,txPlat);
    [rxpos,rxvel,rxax] = helperUpdatePosition(scene,rxPlat);
    [tgtpos,tgtvel,tgtax] = helperUpdatePosition(scene,tgtPlat);
end

Plot the raw received I/Q signal. The direct path is the strongest return in the datacube as expected.

% Visualize raw received signal
figure()
plot(mag2db(abs(sum(y,2))))
grid on
axis tight
xlabel('Samples')
ylabel('Magnitude (dB)')
title('Raw Received Signal')

Figure contains an axes object. The axes object with title Raw Received Signal, xlabel Samples, ylabel Magnitude (dB) contains an object of type line.

Next, plot the raw received I/Q signal for a single pulse repetition interval (PRI). The black dot dash line indicates the start index of the direct path. The red dashed line indicates the index of the target return. The target is not easily visible within a single PRI.

% Calculate direct path start sample 
timeVec = (0:(numSamples - 1))*1/fs;
tDP = range2time(dpRng/2,c);
[~,startIdxDp] = min(abs(timeVec - tDP));

% Calculate target start sample
tgtRng = Rtx + Rrx - dpRng;
tTgt = range2time(tgtRng/2,c);
[~,tgtSampleIdx] = min(abs(timeVec - tTgt));

% Visualize raw received signal including direct path and target sample
% positions
figure()
plot(mag2db(abs(sum(y(startIdx:endIdx,:),2))))
xline(startIdxDp,'-.k',LineWidth=1.5)
xline(tgtSampleIdx + startIdxDp,'--',SeriesIndex=2,LineWidth=1.5)
hold on
grid on
axis tight
xlabel('Samples')
ylabel('Magnitude (dB)')
title('Raw Received Signal, 1 PRI')
legend('Raw Signal','Direct Path','Target Truth')

Figure contains an axes object. The axes object with title Raw Received Signal, 1 PRI, xlabel Samples, ylabel Magnitude (dB) contains 3 objects of type line, constantline. These objects represent Raw Signal, Direct Path, Target Truth.

Create Bistatic Datacube

Build the bistatic datacube. Typically, a bistatic datacube is zero-referenced to the direct path. This means that the target bistatic range and Doppler will also be relative to the direct path. The bistatic radar in this case is stationary, so no zeroing of Doppler is required. Decrease the number of pulses by 1 to prevent an incomplete datacube in the last PRI.

This example assumes that the transmitter and receiver are time synchronized and the positions of the transmitter and receiver are known. All subsequent processing will be performed on the bistatic datacube formed in this section.

% Create bistatic datacube
startIdx = startIdxDp; 
endIdx = startIdx + numSamples - 1;
numPulses = numPulses - 1; 
yBi = zeros(numSamples,numElements,numPulses);
for ip = 1:numPulses
    yBi(:,:,ip) = y(startIdx:endIdx,:);
    startIdx = endIdx + 1;
    endIdx = startIdx + numSamples - 1;
end

Visualize the bistatic datacube by summing over the antenna elements in the receive array. The x-axis is slow time (pulses), and the y-axis is fast-time (range). The strong direct path return is noticeable. The red dashed line indicates the sample index of the target return relative to the direct path.

% Visualize bistatic datacube
helperPlotBistaticDatacube(yBi,tgtSampleIdx)

Figure contains an axes object. The axes object with title Bistatic Datacube Sum over Elements, xlabel Pulse, ylabel Sample contains 2 objects of type image, constantline. This object represents Target Truth.

Process Bistatic Datacube

Matched Filtering and Doppler Processing

Perform matched filtering and Doppler processing using the phased.RangeDopplerResponse object. Oversample in Doppler by a factor of 2 and apply a Taylor window with 40 dB sidelobes to improve the Doppler spectrum.

% Perform matched filtering and Doppler processing
rngdopresp = phased.RangeDopplerResponse(SampleRate=fs,...
    PropagationSpeed=c, ...
    DopplerOutput='Speed', ...
    OperatingFrequency=fc, ...
    DopplerFFTLengthSource='Property', ...
    DopplerWindow='Taylor', ...
    DopplerSidelobeAttenuation=40, ...
    DopplerFFTLength=2*numPulses);
mfcoeff = getMatchedFilter(waveform);
[yBi,rngVec,dopVec] = rngdopresp(yBi,mfcoeff);

The output range vector from the phased.RangeDopplerResponse object assumes a monostatic radar configuration (two-way propagation). To update this range vector for bistatic use, multiply by 2.

% Range vector
rngVec = 2*rngVec;

The dimensions of the range and Doppler processed datacube is bistatic range-by-number of elements-by-bistatic Doppler.

Calculate the expected bistatic Doppler using the radialspeed function for each leg of the bistatic path.

% Calculate bistatic Doppler
tgtSpdTx = radialspeed(tgtpos,tgtvel,txpos,txvel);
tgtSpdRx = radialspeed(tgtpos,tgtvel,rxpos,rxvel);
tgtSpd = (tgtSpdTx + tgtSpdRx)/2;

Plot the range-Doppler response and show the true bistatic range and Doppler on the image. Note that the direct path is still very strong and is localized about 0 Doppler due to the fact that the transmitter and receiver are stationary. The target return is still very difficult to distinguish due to the influence of the direct path. The true target location is shown by the red x.

% Visualize matched filtered and Doppler processed data 
helperPlotRDM('Magnitude (dB)',dopVec,rngVec,yBi,tgtSpd,tgtRng)

Mitigating Direct Path Interference

The direct path is considered to be a form of interference. Typically, the direct path return is removed. Mitigation strategies for direct path removal include:

Adding nulls in the known direction of the transmitter
Using a reference beam or antenna
Adaptive filtering in the spatial/time domain

This example removes the direct path by adding nulls in the known direction of the transmitter.

First, the beamwidth of the receive array is about 6.4 degrees.

% Receiver beamwidth 
rxBeamwidth = beamwidth(collector.Sensor,fc)

rxBeamwidth = 
6.3600

Specify a series of surveillance beams between -45 and 45 degrees. Space the beams 5 degrees apart in order to minimize scalloping losses.

% Surveillance beams
azSpace = 5; 
azAng = -45:azSpace:45;

Null the beams in the direction of the direct path. The angle would be known for this example, since this is a cooperative bistatic case.

% Null angle
nullAng = inangDp;

% Calculate true target angle index
[~,tgtBmIdx] = min(abs(azAng - inangTgt(1)));

Calculate and plot the surveillance beams. The red dashed line indicates the angle of the target. The black dot dash line indicates the null angle in the direction of the direct path. The transmitter and receiver in this example are stationary. Thus, the beamforming weights only need to be calculated once.

% Calculate surveillance beams
pos = (0:numElements - 1).*0.5;
numBeams = numel(azAng);
numDop = numel(dopVec);
yBF = zeros(numBeams,numSamples*numDop);
yBi = permute(yBi,[2 1 3]); % <elements> x <bistatic range> x <bistatic Doppler>
yBi = reshape(yBi,numElements,[]); % <elements> x <bistatic range>*<bistatic Doppler>
afAng = -90:90;
af = zeros(numel(afAng),numBeams); 
for ib = 1:numBeams
    if ~((azAng(ib) - azSpace) < nullAng(1) && nullAng(1) < (azAng(ib) + azSpace))
        wts = minvarweights(pos,azAng(ib),NullAngle=nullAng);

        % Plot the array factor 
        af(:,ib) = arrayfactor(pos,afAng,wts);
        
        hold on

        % Beamform
        yBF(ib,:) = wts'*yBi;
    end
end

Figure contains an axes object. The axes object with title Bistatic Range Doppler Map Sum over Elements, xlabel Bistatic Doppler (m/s), ylabel Bistatic Range (km) contains 2 objects of type image, line. One or more of the lines displays its values using only markers This object represents Target Truth.

% Plot beams
helperPlotBeams(afAng,af,nullAng,inangTgt)

Figure contains an axes object. The axes object with title Beams, xlabel Azimuth Angle (deg), ylabel Magnitude (dB) contains 21 objects of type line, constantline. These objects represent Null, Target Truth.

% Reshape datacube
yBF = reshape(yBF,numBeams,numSamples,numDop); % <beams> x <bistatic range> x <bistatic Doppler>
yBF = permute(yBF,[2 1 3]); %  <bistatic range> x <beams> x <bistatic Doppler>

% Convert to power
yPwr = yBF.*conj(yBF);

Plot the datacube in range-angle space and sum over all Doppler bins. The datacube has now undergone:

Matched filtering
Doppler processing
Beamforming
Nulling of the direct path

The target is now visible. The true target position is indicated by the red x. The deep blue vertical line in the plot below indicates the direct path null. Notice that there is still some leakage from the direct path in the adjacent azimuth bin despite the nulling.

% Plot range-angle after beamforming
helperPlotRAM(azAng,rngVec,yPwr,inangTgt(1),tgtRng)

Figure contains an axes object. The axes object with title Bistatic Range Azimuth Map Sum over Doppler, xlabel Azimuth (deg), ylabel Bistatic Range (km) contains 2 objects of type image, line. One or more of the lines displays its values using only markers This object represents Target Truth.

Now plot the datacube in range-Doppler space and sum over all beams. The target is visible and is located where it is expected.

% Plot range-Doppler after beamforming
helperPlotRDM('Power (dB)',dopVec,rngVec,yPwr,tgtSpd,tgtRng)

Figure contains an axes object. The axes object with title Bistatic Range Doppler Map Sum over Beams, xlabel Bistatic Doppler (m/s), ylabel Bistatic Range (km) contains 2 objects of type image, line. One or more of the lines displays its values using only markers This object represents Target Truth.

Detecting the Target using 2D CFAR

Perform cell-averaging constant false alarm rate (CA CFAR) processing in 2 dimensions using the phased.CFARDetector2D. Set the number of guard cells based on the returns from an example target. The number of training cells ideally totals at least 20, which is easily accomplished in the range dimension. The Doppler dimension has limited bin numbers, so the number of training cells is decreased to 5 on each side.

% Detect using 2D CFAR 
numGuardCells    = [10 5]; % Number of guard cells for each side
numTrainingCells = [10 5]; % Number of training cells for each side
Pfa              = 1e-8;  % Probability of false alarm
cfar = phased.CFARDetector2D(GuardBandSize=numGuardCells, ...
    TrainingBandSize=numTrainingCells, ...
    ProbabilityFalseAlarm=Pfa, ...
    NoisePowerOutputPort=true, ...
    OutputFormat='Detection index');

Visualize the range and Doppler cuts for the target. Note that the training gates do not include contributions from the mainlobe of the target return.

% Visualize range cut
[~,tgtDopIdx] = min(abs(dopVec - tgtSpd)); 
helperPlotCut(yPwr(:,tgtBmIdx,tgtDopIdx),rngVec,tgtRng,numGuardCells(1),numTrainingCells(1),'Range','m')
xlim([(tgtRng - 3e3) (tgtRng + 3e3)])

Figure contains an axes object. The axes object with title Bistatic Range Cut, xlabel Bistatic Range (m), ylabel Power (dB) contains 6 objects of type line, constantline, rectangle. This object represents Target Truth.

% Visualize Doppler cut
helperPlotCut(yPwr(tgtSampleIdx,tgtBmIdx,:),dopVec,tgtSpd,numGuardCells(2),numTrainingCells(2),'Doppler','m/s')

Figure contains an axes object. The axes object with title Bistatic Doppler Cut, xlabel Bistatic Doppler (m/s), ylabel Power (dB) contains 6 objects of type line, constantline, rectangle. This object represents Target Truth.

Perform two-dimensional CFAR over the beams. Assume that the angle of a detection is the beam azimuth angle.

% Perform 2D CFAR detection
[dets,noise] = helperCFAR(cfar,yPwr);

Clustering

Initialize a DBSCAN clustering object. The detections are clustered if they are in adjacent bins.

% Cluster 
clust = clusterDBSCAN('Epsilon',1);
idxClust = clust(dets.');

Plot the clusters. Note that there is a false alarm shown in gray. The false alarm was not clustered, since DBSCAN expects at least 3 detections in each cluster.

% Plot clusters
plot(clust,dets.',idxClust);
xlabel('Range Bin')
ylabel('Beam Bin')
zlabel('Doppler Bin')

Figure Clusters contains an axes object. The axes object with title Clusters, xlabel Range Bin, ylabel Beam Bin contains 4 objects of type line, scatter, text. One or more of the lines displays its values using only markers

Return only the highest SNR detection in each cluster.

% Pare down each cluster to a single detection
[detsClust,snr] = helperDetectionMaxSNR(yPwr,dets,noise,idxClust)

detsClust = 3×2

         168        1878
           3           5
          55          18

snr = 1×2

   24.6559   13.9932

Estimating Target Parameters

Perform range and Doppler estimation using the phased.RangeEstimator and phased.DopplerEstimator objects, respectively. These objects perform a quadratic fit to refine peak estimation of the target.

% Refine detection estimates in range and Doppler
rangeEstimator = phased.RangeEstimator; 
dopEstimator = phased.DopplerEstimator;
rngest = rangeEstimator(yBF,rngVec(:),detsClust);
dopest = dopEstimator(yBF,dopVec,detsClust);

Use the refinepeaks function to improve the azimuth angle estimation.

[~,azest] = refinepeaks(yPwr,detsClust.',azAng,Dimension=2);

Output the true and estimated bistatic target range, Doppler, and azimuth.

% True target location
Ttrue = table(tgtRng*1e-3,tgtSpd,inangTgt(1,:),VariableNames={'Bistatic Range (km)','Bistatic Doppler (m/s)','Azimuth (deg)'});
Ttrue = table(Ttrue,'VariableNames',{'Target Truth'}); 
disp(Ttrue)

                             Target Truth                         
    ______________________________________________________________

    Bistatic Range (km)    Bistatic Doppler (m/s)    Azimuth (deg)
    ___________________    ______________________    _____________
                                                                  
          8.3578                  -119.99               -36.79

% Estimated
Test = table(rngest*1e-3,dopest,azest,VariableNames={'Bistatic Range (km)','Bistatic Doppler (m/s)','Azimuth (deg)'});
Test = table(Test,'VariableNames',{'Target Measured'});
disp(Test)

                           Target Measured                        
    ______________________________________________________________

    Bistatic Range (km)    Bistatic Doppler (m/s)    Azimuth (deg)
    ___________________    ______________________    _____________
                                                                  
          8.3523                  -120.01               -36.273   
          93.875                  -566.84                -25.05

Plot the detections in range-Doppler space. The estimated target location very closely matches the truth.

% Plot detections
helperPlotRDM('Power (dB)',dopVec,rngVec,yPwr,tgtSpd,tgtRng,dopest,rngest)

Figure contains an axes object. The axes object with title Bistatic Range Doppler Map Sum over Beams, xlabel Bistatic Doppler (m/s), ylabel Bistatic Range (km) contains 3 objects of type image, line. One or more of the lines displays its values using only markers These objects represent Target Truth, Target Estimate.

Summary

This example outlines a workflow for creating and processing simulated I/Q signals in a cooperative bistatic radar scenario. It illustrates how to calculate and import a bistatic Radar Cross Section (RCS) for a cone-shaped target using the Antenna Toolbox™. The process involves simulating a bistatic datacube, executing range and Doppler processing, mitigating the direct path interference, and conducting target detection and parameter estimation. The demonstrated I/Q simulation and processing techniques can be iteratively applied within a loop to simulate target returns over extended periods.

Helpers

helperPlotScenario

function helperPlotScenario(txPlat,rxPlat,tgtPlat)
% Helper to visualize scenario geometry using theater plotter
% 
% Plots:
%    - Transmitter and receiver orientation
%    - Position of all platforms
%    - Target trajectory 

% Setup theater plotter
tp = theaterPlot(XLim=[-10 10]*1e3,YLim=[-10 10]*1e3,ZLim=[-10 10]*1e3);

% Setup plotters:
%   - Orientation plotters
%   - Platform plotters
%   - Trajectory plotters 
opTx = orientationPlotter(tp,DisplayName='Tx Orientation',...
    LocalAxesLength=4e3,MarkerSize=1);
opRx = orientationPlotter(tp,DisplayName='Rx Orientation',...
    LocalAxesLength=4e3,MarkerSize=1);
plotterTx = platformPlotter(tp,DisplayName='Transmitter',Marker='^',MarkerSize=8,MarkerFaceColor=[0 0.4470 0.7410],MarkerEdgeColor=[0 0.4470 0.7410]);
plotterRx = platformPlotter(tp,DisplayName='Receiver',Marker='v',MarkerSize=8,MarkerFaceColor=[0 0.4470 0.7410],MarkerEdgeColor=[0 0.4470 0.7410]);
plotterTgt = platformPlotter(tp,DisplayName='Target',Marker='o',MarkerSize=8,MarkerFaceColor=[0.8500 0.3250 0.0980],MarkerEdgeColor=[0.8500 0.3250 0.0980]);
trajPlotter = trajectoryPlotter(tp,DisplayName='Target Trajectory',LineWidth=3);

% Plot transmitter and receiver orientations
plotOrientation(opTx,txPlat.Orientation(3),txPlat.Orientation(2),txPlat.Orientation(1),txPlat.Position)
plotOrientation(opRx,rxPlat.Orientation(3),rxPlat.Orientation(2),rxPlat.Orientation(1),rxPlat.Position)

% Plot platforms
plotPlatform(plotterTx,txPlat.Position,txPlat.Trajectory.Velocity,{'Tx'});
plotPlatform(plotterRx,rxPlat.Position,rxPlat.Trajectory.Velocity,{'Rx'});
plotPlatform(plotterTgt,tgtPlat.Position,tgtPlat.Trajectory.Velocity,{'Target'});

% Plot target trajectory 
plotTrajectory(trajPlotter,{(tgtPlat.Position(:) + tgtPlat.Trajectory.Velocity(:).*linspace(0,5000,3)).'})
grid('on')
title('Scenario')
end

helperPlotTarget

function helperPlotTarget(hTgtConeAxes,tgtPlat)
% Helper to plot the target structure in the radar scenario

% Setup theater plot
tp = theaterPlot(XLim=(tgtPlat.Position(1) + [-5 5]),YLim=(tgtPlat.Position(2) + [-5 5]),ZLim=(tgtPlat.Position(3) + [-5 5]));

% Plot platform
plotterTgt = platformPlotter(tp,DisplayName='Target',Marker='o',MarkerSize=8,MarkerFaceColor=[0.8500 0.3250 0.0980],MarkerEdgeColor=[0.8500 0.3250 0.0980]);
plotPlatform(plotterTgt,tgtPlat.Position,tgtPlat.Trajectory.Velocity);
hold on 

% Plot cone structure
f = hTgtConeAxes.Children(1).Children.Faces;
verts = hTgtConeAxes.Children(1).Children.Vertices;
patch(tp.Parent,Faces=f,Vertices=verts + tgtPlat.Position,FaceColor=[0.9290 0.6940 0.1250],EdgeColor='none');
grid on

% Label
title('Target in Scenario')
end

helperUpdatePosition

function [pos,vel,ax] = helperUpdatePosition(scene,plat)
% Helper to get the position, velocity, and local axes for the input
% platform in the global coordinates

poses = platformPoses(scene,'rotmat');
idx = find(plat.PlatformID == [poses.PlatformID]);
pos = poses(idx).Position(:); vel = poses(idx).Velocity(:); ax = poses(idx).Orientation.';
end

helperPlotBistaticDatacube

function helperPlotBistaticDatacube(yBi,tgtSampleIdx)
% Helper plots the bistatic datacube

numSamples = size(yBi,1);
numPulses = size(yBi,3); 

% Plot bistatic datacube 
figure()
imagesc(1:numPulses,1:numSamples,mag2db(abs(squeeze(sum(yBi,2)))))
axis xy
hold on

% Plot target truth line 
yline(tgtSampleIdx,'--',SeriesIndex=2,LineWidth=1.5)
ylim([0 tgtSampleIdx + 100])

% Label
xlabel('Pulse')
ylabel('Sample')
title(sprintf('Bistatic Datacube\nSum over Elements'))
legend('Target Truth')

% Colorbar
hC = colorbar;
hC.Label.String = 'Magnitude (dB)';
end

helperPlotRDM

function helperPlotRDM(labelChar,dopVec,rngVec,resp,tgtSpd,tgtRng,dopest,rngest)
% Helper to plot the range-Doppler map
%   - Includes target range and Doppler
%   - Optionally includes estimated range and Doppler
%
% Sums over elements or beams depending on input label

% Plot RDM
figure
imagesc(dopVec,rngVec*1e-3,mag2db(abs(squeeze(sum(resp,2)))))
axis xy
hold on

% Plot target truth
plot(tgtSpd,tgtRng*1e-3,'x',SeriesIndex = 2,LineWidth=2)

% Colorbar
hC = colorbar;
hC.Label.String = labelChar;

% Labels
xlabel('Bistatic Doppler (m/s)')
ylabel('Bistatic Range (km)')

% Plot title 
if contains(labelChar,'Magnitude')
    title(sprintf('Bistatic Range Doppler Map\nSum over Elements'))
else
    title(sprintf('Bistatic Range Doppler Map\nSum over Beams'))
end

% Optionally plot estimated range and Doppler 
if nargin == 8
    plot(dopest,rngest*1e-3,'sw',LineWidth=2,MarkerSize=10)
    ylim([tgtRng - 1e3 tgtRng + 1e3]*1e-3)
    legend('Target Truth','Target Estimate',Color=[0.8 0.8 0.8])
else
    ylim([0 tgtRng + 2e3]*1e-3)
    legend('Target Truth');
end
end

helperPlotBeams

function helperPlotBeams(afAng,af,nullAng,inangTgt) 
% Helper to plot:
%    - Surveillance beams
%    - Null angle
%    - Target angle

% Plot array factor 
figure;
hArrayAxes = axes; 
plot(hArrayAxes,afAng,mag2db(abs(af)),SeriesIndex=1)
hold on
grid on
axis tight
ylim([-60 5])

% Label
xlabel('Azimuth Angle (deg)')
ylabel('Magnitude (dB)')
title('Beams')

% Null angle in direction of direct path 
h(1) = xline(nullAng(1),'-.k',LineWidth=1.5);

% Target angle
h(2) = xline(inangTgt(1),'--',LineWidth=1.5,SeriesIndex = 2);
legend(h,{'Null','Target Truth'},Location='southwest')
end

helperPlotRAM

function helperPlotRAM(azAng,rngVec,yPwr,tgtAng,tgtRng)
% Helper to plot the range-angle map
%   - Includes target range and Doppler
%
% Sums over Doppler

% Plot RAM
figure
imagesc(azAng,rngVec*1e-3,pow2db(abs(squeeze(sum(yPwr,3)))))
axis xy
hold on

% Plot target truth
plot(tgtAng,tgtRng*1e-3,'x',SeriesIndex = 2,LineWidth=2)
ylim([0 tgtRng + 2e3]*1e-3)

% Labels
xlabel('Azimuth (deg)')
ylabel('Bistatic Range (km)')
title(sprintf('Bistatic Range Azimuth Map\nSum over Doppler'))
legend('Target Truth')

% Colorbar
hC = colorbar;
hC.Label.String = 'Power (dB)';
end

helperPlotCut

function helperPlotCut(sig,vec,tgtPos,numGuardCells,numTrainingCells,labelChar,unitsChar)
% Helper to plot range and Doppler cuts

% Plot signal
figure
plot(vec,pow2db(abs(squeeze(sig))))
hold on
grid on

% Plot truth
h = xline(tgtPos,'--',SeriesIndex = 2,LineWidth=1.5);

% Guard cells
spacing = vec(end) - vec(end - 1); 
ylims = ylim; 
x = tgtPos - numGuardCells*spacing;
y = ylims(1); 
width = numGuardCells*spacing;
height = ylims(end) - ylims(1); 
rectangle(Position=[x y width height],FaceColor=[0.6350 0.0780 0.1840],FaceAlpha=0.3,EdgeColor='none'); 
x = tgtPos;
rectangle(Position=[x y width height],FaceColor=[0.6350 0.0780 0.1840],FaceAlpha=0.3,EdgeColor='none'); 

% Training cells
x = tgtPos - (numTrainingCells + numGuardCells)*spacing;
width = numTrainingCells*spacing;
rectangle(Position=[x y width height],FaceColor=[0.4660 0.6740 0.1880],FaceAlpha=0.3,EdgeColor='none'); 
x = tgtPos + numGuardCells*spacing;
rectangle(Position=[x y width height],FaceColor=[0.4660 0.6740 0.1880],FaceAlpha=0.3,EdgeColor='none'); 

% Labels
title(['Bistatic ' labelChar ' Cut'])
xlabel(['Bistatic ' labelChar ' (' unitsChar ')'])
ylabel('Power (dB)')
legend(h,{'Target Truth'})
axis tight
end

helperCFAR

function [dets,noise] = helperCFAR(cfar,yPwr)
% Helper to determine CUT indices, as well as perform 2D CFAR detection over
% beams

numGuardCells = cfar.GuardBandSize;
numTrainingCells = cfar.TrainingBandSize;

% Determine cell-under-test (CUT) indices
[numSamples,numBeams,numDop] = size(yPwr); 
numCellsRng       = numGuardCells(1) + numTrainingCells(1); 
idxCUTRng         = (numCellsRng+1):(numSamples - numCellsRng);  % Range indices of cells under test (CUT) 
numCellsDop       = numGuardCells(2) + numTrainingCells(2); 
idxCUTDop         = (numCellsDop + 1):(numDop - numCellsDop); % Doppler indices of cells under test (CUT) 
[idxCUT1,idxCUT2] = meshgrid(idxCUTRng,idxCUTDop);
idxCUT = [idxCUT1(:) idxCUT2(:)].';

% Perform CFAR detection over beams 
dets = []; 
noise = []; 
for ib = 1:numBeams % Perform detection over beams 
    [detTmp,noiseTmp] = cfar(squeeze(yPwr(:,ib,:)),idxCUT);
    beamTmp = ib.*ones(1,size(detTmp,2)); 
    detTmp = [detTmp(1,:); beamTmp; detTmp(2,:)];
    dets = [dets detTmp]; %#ok<AGROW>
    noise = [noise noiseTmp]; %#ok<AGROW>
end
end

helperDetectionMaxSNR

function [detsClust,snr] = helperDetectionMaxSNR(yPwr,dets,noise,idxClust)
% Helper that outputs the maximum SNR detection in each cluster

idxU = unique(idxClust);
detsClust = []; 
snr = [];
for iu = 1:numel(idxU)
    if idxU(iu) ~= -1
        idx = find(idxClust == idxU(iu)); 

        % Find max SNR
        linidx = sub2ind(size(yPwr),dets(1,idx),dets(2,idx),dets(3,idx));
        snrs = yPwr(linidx)./noise(:,idx);
        [~,idxMax] = max(snrs); 

        % Save dets 
        detsClust = [detsClust [dets(1,idx(idxMax));dets(2,idx(idxMax));dets(3,idx(idxMax))]]; %#ok<AGROW>
        snr = [snr yPwr(linidx(idxMax))./noise(:,idx(idxMax))]; %#ok<AGROW>
    end
end
snr = pow2db(snr); 
end