Analysis and Simulation of a Low Probability of Intercept Radar System

Since R2025a

This example models a monostatic radar system and a non-cooperative receiver that aims to intercept signals transmitted by the radar. The example begins with a power-level analysis to determine conditions under which a radar system has low probability of intercept (LPI) characteristics and when it becomes susceptible to interception. Next, an equivalent waveform-level simulation that generates I/Q signals validates the predictions from the power-level analysis. This demonstrates that regardless of the transmitted waveform type, the radar system can always be intercepted as long as a sufficient amount of its transmitted energy is within the interceptor's bandwidth. Finally, the example shows that frequency hopping is an effective strategy to counteract interception.

Introduction

The goal of an LPI radar is to evade detection by non-cooperative intercept receivers. To achieve this, LPI radars transmit waveforms with large time-bandwidth products and maintain low peak transmit power. This is typically accomplished by using continuous-wave (CW) waveforms, such as frequency-modulated continuous-wave (FMCW) waveform. However, the performance of an LPI radar also depends on the characteristics of the intercept receiver. A radar may exhibit LPI properties against one type of interceptor but not against another [1].

This example begins by defining the parameters of a CW radar and calculating its maximum detection range. Next, it defines the parameters of an intercept receiver and introduces the concept of intercept receiver range advantage, a key performance metric for LPI radars. The intercept receiver range advantage is calculated as the ratio of the interceptor range to the radar's maximum detection range. The example shows how to determine the intercept receiver's required processing gain to ensure that the range advantage exceeds one, indicating that the radar system can be intercepted beyond its maximum detection range. Additionally, the example explains how to compute the non-coherent integration time necessary for the intercept receiver to achieve this processing gain.

The remainder of the example involves waveform-level I/Q simulations, which confirm that the interceptor achieves the required probability of detecting the radar waveform within the calculated non-coherent integration window, regardless of the transmitted waveform type. Finally, the example shows that the radar waveform interception can be avoided if the radar employs frequency agility.

Radar System

Consider a continuous wave (CW) radar system located at the origin.

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

Radar.Position = [0; 0; 0];
Radar.Velocity = [0; 0; 0];    % The radar is static

The radar system parameters used in this example are selected based on the PILOT "quiet" radar developed by the Philips Research Laboratory [2]. Let the radar have a peak transmit power of 1 W, operating frequency of 9.375 GHz, and a sweep bandwidth of 50 MHz.

Radar.PeakPower = 1;
Radar.OperatingFrequency = 9.375e9;
Radar.Wavelength = freq2wavelen(Radar.OperatingFrequency);
Radar.Bandwidth = 50e6;

Set the sweep time to 1 ms.

Radar.SweepTime = 1e-3;

Compute the time-bandwidth product of the transmitted waveform.

timeBandwidthProduct = Radar.SweepTime*Radar.Bandwidth;
fprintf('Time-bandwidth product: %.2f\n', timeBandwidthProduct);

Time-bandwidth product: 50000.00

Set the radar receiver noise figure to 5 dB and the reference noise temperature to 290 K.

Radar.NoiseFigure = 5.0;
Radar.ReferenceTemperature = 290;

Compute the equivalent system noise temperature and the noise power.

Radar.SystemTemperature = systemp(Radar.NoiseFigure, Radar.ReferenceTemperature);
Radar.NoisePower = noisepow(Radar.Bandwidth, Radar.NoiseFigure, Radar.ReferenceTemperature);

The helperGetAntennaPattern helper function returns an antenna pattern with a gain of 32 dB and a sidelobe level of approximately -30 dB, which is close to the parameters of the referenced PILOT radar. Use the generated pattern and the phased.CustomAntennaElement to model the radar antenna.

[patmag, pataz, patel] = helperGetAntennaPattern(Radar.OperatingFrequency);
Radar.Antenna = phased.CustomAntennaElement(AzimuthAngles=pataz, ElevationAngles=patel, MagnitudePattern=patmag, PhasePattern=zeros(size(patmag)));

Compute the antenna gain.

Gtx = directivity(Radar.Antenna, Radar.OperatingFrequency, [0; 0]);
Grx = Gtx;
fprintf('Tx/Rx antenna gain: %.2f dB\n', Gtx);

Tx/Rx antenna gain: 32.00 dB

Compute the maximum range of this radar system, assuming it is designed to detect Swerling 1 targets with a radar cross section of 10 m², with a detection probability of 0.9 and the probability of false alarm of $10^{- 6}$ .

Target = phased.RadarTarget(OperatingFrequency=Radar.OperatingFrequency);
Target.MeanRCS = 10;
Target.Model = "Swerling1";

Radar.Pd = 0.9;
Radar.Pfa = 1e-6;

Given the selected radar system parameters and the target detection requirements, the maximum detection range can be computed using the radar equation [3]

$R_{m a x} = {[\frac{P_{t} τ G_{t x} G_{r x} λ^{2} σ}{(4 π)^{3} D_{x} k T_{s r} L_{r}}]}^{1 / 4}$

where $P_{t}$ is the peak transmit power, $τ$ is the coherent processing interval (CPI), $G_{tx}$ is the transmit antenna gain, $G_{rx}$ is the receive antenna gain, $λ$ is the wavelength, $σ$ is the target radar cross section (RCS), $k$ is the Boltzmann's constant, $T_{sr}$ is the radar receiver system temperature, $D_{x}$ is the detectability factor, and $L_{r}$ is the cumulative loss term. This equation can be modified to explicitly include the radar bandwidth $B_{r}$ and the receiver processing gain $β_{r}$ .

$R_{m a x} = {[\frac{P_{t} β_{r} G_{t x} G_{r x} λ^{2} σ}{(4 π)^{3} D_{x} k T_{s r} B_{r} L_{r}}]}^{1 / 4}$

In the case when the CPI is equal to the sweep time, the radar processing gain is equal to the time-bandwidth product.

betar = timeBandwidthProduct;

Assuming the cumulative losses are 0 dB ( $L_{r} = 1$ ), compute the maximum detection range $R_{\max}$ .

% Detectability factor
Dx = detectability(Radar.Pd, Radar.Pfa, 1, Target.Model);

% Radar maximum range
Rmax = radareqrng(Radar.Wavelength, Dx, Radar.PeakPower, Radar.SweepTime,...
    "Gain", [Gtx Grx], "RCS", Target.MeanRCS, "Ts", Radar.SystemTemperature);
fprintf('Radar maximum detection range: %.2f km\n', Rmax*1e-3);

Radar maximum detection range: 9.42 km

Intercept Receiver

Let the intercept receiver be located at an angle of 12 degrees in azimuth with respect to the radar.

interceptorAngle = [12; 0];

Compute the transmit antenna gain $\tilde{G_{tx}}$ in the direction of the intercept receiver.

Gtxi = directivity(Radar.Antenna, Radar.OperatingFrequency, interceptorAngle);
fprintf('Transmit antenna gain in the direction of the interceptor is %.2f dB or %.2f dB below the mainlobe peak\n', Gtxi, (Gtx-Gtxi));

Transmit antenna gain in the direction of the interceptor is 0.83 dB or 31.17 dB below the mainlobe peak

To investigate how the probability of intercepting a radar waveform changes with distance between the radar and the interceptor, consider intercept receiver located at ranges that vary from two kilometers before $R_{\max}$ to four kilometers beyond $R_{\max}$ .

interceptorRange = 1e3*(-2:0.5:4) + Rmax;
numInterceptorRanges = numel(interceptorRange);

Assume the intercept receiver platform is static.

Interceptor.Velocity = [0; 0; 0];

Visualize this scenario using the helperPlotLPIScenario helper function.

helperPlotLPIScenario(Radar, Rmax, interceptorRange, interceptorAngle);
xlim([-0.1 1.5]*Rmax*1e-3);
ylim([-0.4 1.2]*Rmax*1e-3);

$Figure contains an axes object. The axes object with title LPI Scenario, xlabel X (km), ylabel Y (km) contains 5 objects of type line. One or more of the lines displays its values using only markers These objects represent Bearing towards intercept receiver, Intercept receiver positions, R_{max}, Radar antenna pattern, Radar.$

In general, a non-cooperative intercept receiver has no knowledge of the radar operating frequency and bandwidth. Let the interceptor have an operating frequency of 9.39 GHz, which is 15 MHz higher than the radar operating frequency.

Interceptor.OperatingFrequency = 9.39e9;

The intercept receiver also must have a much wider bandwidth compared to the radar. Let the interceptor bandwidth be 150 MHz.

% Interceptor bandwidth is larger than the radar bandwidth
Interceptor.Bandwidth = 150e6;

Let the intercept receiver have a noise figure of 5 dB and a reference noise temperature of 290 K.

Interceptor.NoiseFigure = 5.0;
Interceptor.ReferenceTemperature = 290;

Compute the equivalent system noise temperature and the noise power.

Interceptor.SystemTemperature = systemp(Interceptor.NoiseFigure, Interceptor.ReferenceTemperature);
Interceptor.NoisePower = noisepow(Interceptor.Bandwidth, Interceptor.NoiseFigure, Interceptor.ReferenceTemperature);

For simplicity, assume the intercept receiver has an omnidirectional antenna.

Interceptor.Antenna = phased.IsotropicAntennaElement;

% Interceptor antenna gain
Gi = 0;

The equation for determining the maximum intercept receiver range resembles the radar equation, but with a key difference: due to the one-way propagation, it features an exponent of (1/2) instead of (1/4). Additionally, it excludes the target's radar cross-section (RCS) because it models a direct path between the radar and the intercept receiver.

$R_{i} = {[\frac{P_{t} β_{i} {G_{}^{\sim}}_{t x} G_{i} λ^{2}}{(4 π)^{2} S N R_{i} k T_{s i} B_{i} L_{i}}]}^{1 / 2}$

Here ${SNR}_{i}$ is a signal-to-noise ratio required to detect a presence of a radar waveform with required probabilities of detection and false alarm, ${\tilde{G}}_{rx}$ is the radar transmit antenna gain in the direction of the interceptor, $β_{i}$ is the intercept receiver processing gain, $B_{i}$ is the interceptor bandwidth, $T_{si}$ is the intercept receiver system noise temperature, and $L_{i}$ is the cumulative losses at the intercept receiver.

Since the interceptor is trying to detect the presence of a radar signal arriving through the direct path, the signal-plus-noise voltage can be modeled using the Rician distribution. This means that the detectability equations for a non-fluctuating target case can be used to compute the required SNR at the interceptor. Assume that the desired probability of detecting the radar waveform is 0.9 and the required probability of false alarm is $10^{- 6}$ . Compute the required SNR.

Interceptor.Pd = 0.9;
Interceptor.Pfa = 1d-6;
SNRi = detectability(Interceptor.Pd, Interceptor.Pfa, 1);

The processing gain $β_{i}$ depends on the signal processing algorithm used by the interceptor. Potentially, $β_{i}$ could be as large as the waveform's time-bandwidth product. However, since the interceptor has no knowledge of the transmitted radar waveform, in practice $β_{i}$ is considerably smaller. If the interceptor does not perform any integration and tries to detect the presence of the radar waveform at a single time instance, then $β_{i} = 1$ .

Intercept Receiver Range Advantage

Whether or not a radar system is LPI can be assessed by computing the ratio of the maximum interceptor range to the maximum detector range. This ratio is called the intercept receiver range advantage [1, 4].

$\frac{R_{i}}{R_{m a x}} = R_{m a x} {[\frac{4 π}{σ} \frac{1}{δ} \frac{{G_{}^{\sim}}_{t x} G_{i}}{G_{t x} G_{r x}} \frac{L_{r}}{L_{i}}]}^{1 / 2}$

where $δ$ is a ratio of the intercept receiver sensitivity to the radar receiver sensitivity

$δ = \frac{k T_{s i} S N R_{i}}{k T_{s r} D_{x}} \cdot \frac{B_{i}}{B_{r}} \cdot \frac{β_{r}}{β_{i}}$

Note, that $δ$ is larger for radar receivers with higher sensitivity and a greater processing gain. It also increases when the intercept receiver bandwidth increases, which reduces the intercept receiver sensitivity. Use the helperInterceptorRangeAdvantage helper function to plot the intercept receiver range advantage against the radar's maximum range for different values of $δ$ .

% Radar maximum range values
Rr = [0.1 (1:100)]*1e3;

% Receiver sensitivity ratio
delta = [1 10 100];
helperInterceptorRangeAdvantage(Rr, Target.MeanRCS, delta, Gtx, Gtxi, Gi);

The shaded area in the plot corresponds to the region where the intercept receiver range advantage is less than one. This means that the interceptor range is shorter than the radar's maximum range, and the interceptor can potentially be detected by the radar before the radar can be intercepted. In other words, if the interceptor is inside the shaded area, the radar system is LPI.

This diagram illustrates that a radar's status as an LPI system is heavily influenced by the intercept receiver to the radar receiver sensitivity ratio $δ$ , as well as whether the interception occurs via the mainlobe or sidelobe.

For a radar system considered in this example, an intercept receiver that does not perform any integration ( $β_{i} = 1$ ) is well inside the LPI area (the intercept receiver range advantage is much less than one). In this case the receiver sensitivity ratio is very large indicating the radar receiver advantage over the interceptor.

% Intercept receiver processing gain
betai1 = 1;

% Compute receiver sensitivity ratio using the
% helperReceiverSensitivityRatio helper function
delta1 = helperReceiverSensitivityRatio(SNRi, Dx, Interceptor.Bandwidth, Radar.Bandwidth,...
    betai1, betar, Interceptor.SystemTemperature, Radar.SystemTemperature);

% Intercept receiver range advantage
rangeAdvantage1 = Rmax*sqrt((4*pi/delta1) * (1/Target.MeanRCS) * db2pow(Gtxi + Gi - (Gtx + Grx)));

plot(Rmax, rangeAdvantage1, 'pentagram', LineWidth=2, MarkerEdgeColor='auto',....
    DisplayName=sprintf('\\delta = %.2f, \\beta_i = %.2f dB', delta1, pow2db(betai1)));

To determine the minimum processing gain needed for the intercept receiver range advantage to be greater than one, set $R_{i} = R_{\max}$ in the interceptor range equation and solve for $β_{i}$ .

$β_{i} = \frac{(4 π)^{2} S N R_{i} k T_{s i} B_{i} R_{m a x}^{2}}{P_{t} {G_{}^{\sim}}_{t x} G_{i} λ^{2}}$

betai = (4*pi)^2*Interceptor.NoisePower*db2pow(SNRi)*Rmax^2./(Radar.PeakPower*db2pow(Gtxi + Gi)*Radar.Wavelength^2);
fprintf('Required processing gain at interceptor assuming %.2f sidelobe intercept: %.2f \n', -(Gtx-Gtxi), betai);

Required processing gain at interceptor assuming -31.17 sidelobe intercept: 440.27

% Corresponding receiver sensitivity ratio
delta = helperReceiverSensitivityRatio(SNRi, Dx, Interceptor.Bandwidth, Radar.Bandwidth,...
    betai, betar, Interceptor.SystemTemperature, Radar.SystemTemperature);

% Intercept receiver range advantage
rangeAdvantage = Rmax*sqrt((4*pi/delta) * (1/Target.MeanRCS) * db2pow(Gtxi + Gi - (Gtx + Grx)));
plot(Rmax, rangeAdvantage, 'pentagram', LineWidth=2, MarkerEdgeColor='auto',...
    DisplayName=sprintf('\\delta = %.2f, \\beta_i = %.2f dB', delta, pow2db(betai)));

$Figure contains an axes object. The axes object with title Intercept Receiver Range Advantage vs. Radar Maximum Detection Range, xlabel R indexOf max baseline blank (m), ylabel R indexOf i baseline /R indexOf max baseline contains 9 objects of type line, patch. One or more of the lines displays its values using only markers These objects represent Sidelobe intercept, \delta = 1.0, Mainlobe intercept, \delta = 1.0, Sidelobe intercept, \delta = 10.0, Mainlobe intercept, \delta = 10.0, Sidelobe intercept, \delta = 100.0, Mainlobe intercept, \delta = 100.0, LPI, \delta = 23653.54, \beta_i = 0.00 dB, \delta = 53.73, \beta_i = 26.44 dB.$

Based on the calculation, a processing gain of 440.3 (26.44 dB) is sufficient for the intercept receiver to detect radar signals beyond the radar's maximum detection range, indicating that the radar is no longer operating as an LPI system. It's important to note that this required gain is significantly smaller than the radar's processing gain. If the interceptor could perform coherent integration, it would need to integrate over less than 1% (440/50,000 < 0.01) of the radar's sweep time. However, since the interceptor lacks knowledge of the radar waveform, it can only enhance processing gain through noncoherent integration, which necessitates a much longer integration window.

To determine the appropriate duration for noncoherent integration by the intercept receiver, estimate the integration gain achieved by combining samples noncoherently over various time window lengths. Additionally, calculate the integration gain using an M-of-N integration scheme, where M = 0.5N, and N represents the total number of samples collected by the intercept receiver within the specified time interval.

% Integration window duration
taui = linspace(1e-4*Radar.SweepTime, Radar.SweepTime, 100);

% Number of samples collected by the interceptor
n = ceil(taui*Interceptor.Bandwidth);

Gnc = zeros(size(taui));
Gmn = zeros(size(taui));

for it = 1:numel(n)
    % Noncoherent integration gain
    Gnc(it) = SNRi - detectability(Interceptor.Pd, Interceptor.Pfa, n(it));

    % M-of-N integration gain
    Gmn(it) = Gnc(it) - binaryintloss(Interceptor.Pd, Interceptor.Pfa, n(it), ceil(0.5*n(it)));
end

Compare the computed noncoherent and M-of-N integration gains with the maximum achievable coherent integration gain.

% Corresponding coherent integration gain achievable only if the
% interceptor knows the waveform
Gc = pow2db(taui*Radar.Bandwidth);

figure;
semilogx(taui, Gc, LineWidth=2, DisplayName='Coherent integration');
hold on;
semilogx(taui, Gnc, LineWidth=2, DisplayName='Noncoherent integration');
semilogx(taui, Gmn, LineWidth=2, DisplayName='M-of-N integration');
grid on;
xlim([1e-5 1e-3]);
ylim([0 pow2db(Radar.Bandwidth*1e-3)]);

xlabel('Integration time (s)');
ylabel('Gain (dB)');
yline(pow2db(betai), '--', sprintf('%.1f dB',pow2db(betai)), LineWidth=2,...
    LabelHorizontalAlignment='right', LabelVerticalAlignment='bottom');

yyaxis(gca(),'right');
set(gca(),'ycolor','k')
ylabel('Integration efficiency')
legend('Coherent integration', 'Noncoherent integration', 'M-of-N integration', 'Gain required for R_i/R_{max}>1',...
    Location='southwest');
title('Intercept Receiver Integration Gain vs. Integration Time')

Figure contains an axes object. The axes object with title Intercept Receiver Integration Gain vs. Integration Time, xlabel Integration time (s), ylabel Integration efficiency contains 4 objects of type line, constantline.

Note that the computed gain accounts for the additional gain that results from the interceptor sampling the received signal at a rate higher than the radar waveform bandwidth. This gain arises because the noise is uncorrelated from sample to sample, allowing it to average out during integration, while the signal samples remain correlated and add constructively.

Calculate the duration of the integration window needed to achieve the desired gain that would result in an interceptor range advantage that is greater than one.

tauinc = taui(Gnc >= pow2db(betai));
tauinc = tauinc(1);
fprintf('Integration window duration for noncoherent integration: %.2f ms \n', tauinc*1e3);

Integration window duration for noncoherent integration: 0.11 ms

tauimn = taui(Gmn >= pow2db(betai));
tauimn = tauimn(1);
fprintf('Integration window duration for M-of-N integration: %.2f ms \n', tauimn*1e3);

Integration window duration for M-of-N integration: 0.24 ms

The M-of-N integration requires a window that is approximately twice as long as that of noncoherent integration. However, this window remains significantly shorter than the sweep time, allowing the intercept receiver to potentially perform multiple M-of-N integrations within a single sweep. In addition, the M-of-N scheme may offer greater robustness, as it does not rely on 100% of the received signal within a chosen integration window containing the radar waveform.

Radar and Intercept Receiver Simulation

The rest of the example performs waveform-level I/Q simulation of the radar and the intercept receiver to verify that the selected integration window at the intercept receiver is sufficient.

Introduce a target into the simulation to ensure that radar detection operates as expected. Position the target along the x-axis at a distance equal to the radar's maximum range.

targetRange = Rmax;
targetAngle = [0; 0];
targetPosition = targetRange*[cosd(targetAngle(1)); sind(targetAngle(1)); 0];
targetVelocity = [0; 0; 0];
helperPlotLPIScenario(Radar, Rmax, interceptorRange, interceptorAngle, targetRange, targetAngle);
xlim([-0.1 1.5]*Rmax*1e-3);
ylim([-0.4 1.2]*Rmax*1e-3);

$Figure contains an axes object. The axes object with title LPI Scenario, xlabel X (km), ylabel Y (km) contains 6 objects of type line. One or more of the lines displays its values using only markers These objects represent Bearing towards intercept receiver, Intercept receiver positions, R_{max}, Radar antenna pattern, Radar, Target.$

Received Signal Simulation

Create radar transmitter, radiator, collector, and receiver objects that model the radar system.

% Radar transmitter
Radar.Transmitter = phased.Transmitter(Gain=0, PeakPower=Radar.PeakPower);
Radar.Radiator = phased.Radiator(Sensor=Radar.Antenna, OperatingFrequency=Radar.OperatingFrequency,...
    SensorGainMeasure="dBi");

% Radar receiver
Radar.Collector = phased.Collector(Sensor=Radar.Antenna, OperatingFrequency=Radar.OperatingFrequency);
Radar.Receiver = phased.Receiver(SampleRate=Radar.Bandwidth, Gain=0, AddInputNoise=true, ...
    NoiseMethod="Noise figure", NoiseFigure=Radar.NoiseFigure, ...
    InputNoiseTemperature=Radar.ReferenceTemperature, ReferenceTemperature=Radar.ReferenceTemperature);

Create collector and receiver objects to model the intercept receiver.

Interceptor.Collector = phased.Collector(Sensor=phased.IsotropicAntennaElement,...
    OperatingFrequency=Interceptor.OperatingFrequency);

Interceptor.Receiver = phased.Receiver(SampleRate = Interceptor.Bandwidth, Gain=0, AddInputNoise=true,...
    NoiseMethod="Noise figure", NoiseFigure=Interceptor.NoiseFigure,...
    InputNoiseTemperature=Interceptor.ReferenceTemperature, ReferenceTemperature=Interceptor.ReferenceTemperature);

Note that the radar receiver sample rate is set to the radar's bandwidth, but the intercept receiver sample rate is set to the intercept receiver bandwidth. This is important for correctly modeling the amount of thermal noise added to the received signal. Given that the intercept receiver has a significantly larger bandwidth than the radar, it will also experience more noise.

Since the radar and the intercept receivers have different sample rates, the initial transmitted waveform is sampled at a much higher rate. A higher sample rate also allows for modeling radar frequency agility.

sampleRate = 600e6;
radarDownsampleFactor = sampleRate/Radar.Bandwidth;
interceptorDownsampleFactor = sampleRate/Interceptor.Bandwidth;

Create a transmitted waveform.

% Switch to PMCW waveform to verify that the results are independent of the
% transmitted waveform
waveformType = "FMCW";

if waveformType == "FMCW"
    waveform = phased.FMCWWaveform(SweepBandwidth=Radar.Bandwidth, SampleRate=sampleRate, ...
        SweepInterval="Symmetric", SweepTime=Radar.SweepTime, NumSweeps=1, SweepDirection="Up");
else
    % Phase coded waveform based on a quadratic residue sequence of length
    % 49999 (time-bandwidth product is approximately 50000)
    numChips = 49999;
    prf = Radar.Bandwidth/numChips;
    waveform = phased.PhaseCodedWaveform(ChipWidth=1/Radar.Bandwidth, SampleRate=sampleRate, Code="Quadratic Residue Sequence",...
        NumChips=49999, NumPulses=1, PRF=prf);
end

Generate waveform samples.

wav = waveform();
t = (0:numel(wav)-1).'/sampleRate;

Create a free space channel object to model a two-way radar signal propagation from the radar transmitter to the target and back to the radar receiver.

radarTargetChannel = phased.FreeSpace(SampleRate=sampleRate, TwoWayPropagation=true,...
    OperatingFrequency=Radar.OperatingFrequency);

Create a free space channel object to model a one-way radar signal propagation from the radar transmitter to the intercept receiver.

radarInterceptorChannel = phased.FreeSpace(SampleRate=sampleRate, TwoWayPropagation=false,...
    OperatingFrequency=Radar.OperatingFrequency);

Simulate radar waveform propagation to the target and back to the radar receiver.

% Generate radar transmit waveform
transmittedSignal = Radar.Transmitter(wav);

% Radiate the waveform by the radar transmit antenna towards the target
radiatedToTargetSignal = Radar.Radiator(transmittedSignal, targetAngle);

% Propagate the waveform through the two-way radar propagation channel
propagatedToTargetSignal = radarTargetChannel(radiatedToTargetSignal, Radar.Position, targetPosition, Radar.Velocity, targetVelocity);

% To ensure that there is a detection simulate multiple instances of the
% received radar signal
numRadarTrials = 5;
receivedRadarSignal = zeros(numel(wav)/radarDownsampleFactor, numRadarTrials);

for it = 1:numRadarTrials
    % Reflect the propagated signal off the target
    reflectedToRadarSignal = Target(propagatedToTargetSignal, false);

    % Collect the return signal by the radar receive antenna
    collectedAtRadarSignal = Radar.Collector(reflectedToRadarSignal, targetAngle);

    % Change the sampling rate of the signal to be equal to the radar bandwidth
    % before passing it through the receiver
    collectedAtRadarSignal = downsample(collectedAtRadarSignal, radarDownsampleFactor);

    % Pass the collected signal through the receiver
    receivedRadarSignal(:, it) = Radar.Receiver(collectedAtRadarSignal);
end

Loop over the set of intercept receiver ranges and simulate radar waveform propagation to the intercept receiver.

% Radiate the waveform towards the intercept receiver
radiatedToInterceptorSignal = Radar.Radiator(transmittedSignal, interceptorAngle);

% Preallocate space for radar signals collected at different intercept
% receiver positions
collectedAtInterceptorSignal = zeros(numel(wav), numInterceptorRanges);

for ir = 1:numInterceptorRanges
    % Reset the channel each time
    reset(radarInterceptorChannel);
    Interceptor.Position = interceptorRange(ir)*[cosd(interceptorAngle(1)); sind(interceptorAngle(1)); 0];

    % Propagate the waveform through the one-way propagation channel
    % towards the interceptor
    propagatedToInterceptorSignal = radarInterceptorChannel(radiatedToInterceptorSignal, Radar.Position, Interceptor.Position, Radar.Velocity, Interceptor.Velocity);
    collectedAtInterceptorSignal(:, ir) = Interceptor.Collector(propagatedToInterceptorSignal, interceptorAngle);
end

Apply frequency offset due to the difference between the radar and the intercept receiver operating frequencies.

% Frequency offset between the radar and the intercept receiver
frequencyOffset = Radar.OperatingFrequency - Interceptor.OperatingFrequency;
interceptedSignal = collectedAtInterceptorSignal.*exp(1i*2*pi*frequencyOffset.*t);

Bandpass filter the signal to fit the intercept receiver bandwidth. Then resample it to the sample rate of the intercept receiver.

interceptedSignal = lowpass(interceptedSignal, Interceptor.Bandwidth/2, sampleRate);
interceptedSignal = downsample(interceptedSignal, interceptorDownsampleFactor);

Finally, pass the signal through the receiver to add the intercept receiver noise. Perform this 50 times to collect data for computing the probability of detection statistics.

numInterceptorTrials = 50;
receivedInterceptedSignal = zeros(size(interceptedSignal, 1), numInterceptorTrials, numInterceptorRanges);
for it = 1:numInterceptorTrials
    receivedInterceptedSignal(:, it, :) = Interceptor.Receiver(interceptedSignal);
end

Radar Receiver Signal Processing

Check that the radar can see the target as expected.

if waveformType == "FMCW"
    % For FMCW waveform process the received radar signal by first
    % dechirping it and then performing an FFT to obtain the range response
    wavref = downsample(wav, radarDownsampleFactor);
    dechirpedRadarSignal = dechirp(receivedRadarSignal, wavref);
    nfft = 2^nextpow2(numel(dechirpedRadarSignal));

    resp = fft(dechirpedRadarSignal, nfft);

    sweepSlope = waveform.SweepBandwidth/waveform.SweepTime;
    fb = ((0:nfft-1).'/nfft)*Radar.Bandwidth;
    rnggrid = beat2range(fb, sweepSlope);

else
    % For PMCW waveform process the received radar signal by applying the
    % matched filter
    matchedFilter = phased.MatchedFilter; 
    matchedFilter.CoefficientsSource = 'input port'; 
    
    coef = getMatchedFilter(waveform);   
    coef = downsample(coef, radarDownsampleFactor);

    % Since the transmitted waveform has 100% duty cycle, the matched
    % filter has a group delay equal to the waveform duration. Thus, the
    % matched filter object must be called twice to see the target.
    matchedFilter(receivedRadarSignal, coef);
    resp = matchedFilter(zeros(size(receivedRadarSignal)), coef);

    rnggrid = time2range((1:size(resp, 1)).'/Radar.Bandwidth);
end

resp = abs(resp).^2;

Compute the detection threshold assuming that a single sweep interval is processed at a time.

detectionTreshold = Radar.NoisePower*timeBandwidthProduct*db2pow(npwgnthresh(Radar.Pfa, 1, 'noncoherent'));
detectionIdxs = resp > detectionTreshold;

Visualize the range response together with the computed detection threshold.

figure;
hold on;
inRangeIdxs = rnggrid <= 1.2*Rmax;
plot(rnggrid(inRangeIdxs)*1e-3, pow2db(resp(inRangeIdxs, :)));

for i = 1:numRadarTrials
    idxs = detectionIdxs(:, i);
    plot(rnggrid(idxs)*1e-3, pow2db(resp(idxs, i)), 'r*');
end

yline(pow2db(detectionTreshold), '--', 'Detection threshold', LabelHorizontalAlignment='left', LabelVerticalAlignment='middle');
colors = gca().ColorOrder;

xline(Rmax*1e-3, '--', 'Rmax', 'Color', 'k', 'LabelHorizontalAlignment', 'center', 'LabelVerticalAlignment', 'top');

xlim([0 1.2*Rmax*1e-3]);
ylim([-90 -50]);
xlabel('Range (km)');

ylabel('(dBW)');
title('Radar Range Response');
grid on;

Figure contains an axes object. The axes object with title Radar Range Response, xlabel Range (km), ylabel (dBW) contains 9 objects of type line, constantline. One or more of the lines displays its values using only markers

As expected, the target located at the radar's maximum detection range generates detection threshold crossings.

Intercept Receiver Signal Processing

The intercept receiver utilizes M-of-N (binary) integration with M = 0.5N, effectively implementing a majority vote rule. Within each integration window, the power of each sample is compared against a threshold. A radar waveform detection is declared if at least half of the samples within the window exceed this threshold.

The computed duration of the integration window that should result in the interceptor range being greater than the radar's maximum detection range is

fprintf('%.2f ms \n', tauimn*1e3);

0.24 ms

Compute N, the number of samples within the integration window, and M, the number of threshold crossings required to declare a detection.

% Number of samples in one integration window
N = ceil(tauimn*Interceptor.Bandwidth);

% Number of binary detection
M = ceil(N*0.5);

Use the binaryintloss function to compute the expected binary integration loss, as well as the detection and false alarm probabilities for a single sample.

[Lb, pd, pfa] = binaryintloss(Interceptor.Pd, Interceptor.Pfa, N, M);

Compute the single sample detection threshold.

interceptorDetectionThreshold = Interceptor.NoisePower * db2pow(npwgnthresh(pfa, 1, 'noncoherent'));

Use the helperEstimatePd helper function to estimate the probability of intercepting the radar waveform based on the generated received signals.

pdInBand = helperEstimateInterceptorPd(receivedInterceptedSignal, N, M, interceptorDetectionThreshold);

Plot the estimated probability of intercepting the radar waveform as a function of the range to the intercept receiver.

figure;
hold on;
plot(interceptorRange*1e-3, pdInBand, 's--', LineWidth=2, MarkerFaceColor='auto');
grid on;
xlabel('Range to intercept receiver (km)');
ylabel('P_d');
title('Probability of Intercepting Radar Waveform')
xline(Rmax*1e-3, '--', 'Rmax', 'LabelHorizontalAlignment', 'center', 'LabelVerticalAlignment', 'top');
ylim([0 1.1]);

% Mark the expected maximum interceptor range computed based on the interceptor range equation. 
SNRiMN = detectability(pd, pfa, 1);
RiMN = sqrt((Radar.PeakPower*db2pow(Gtxi + Gi)*Radar.Wavelength^2)./((4*pi)^2*Interceptor.NoisePower*db2pow(SNRiMN)));
xline(RiMN*1e-3, '--', 'Ri', 'LabelHorizontalAlignment', 'center', 'LabelVerticalAlignment', 'bottom');
yline(Interceptor.Pd, '--', 'Required P_d');

Figure contains an axes object. The axes object with title Probability of Intercepting Radar Waveform, xlabel Range to intercept receiver (km), ylabel P indexOf d baseline P_d contains 4 objects of type line, constantline.

The estimated probability of detecting the presence of the radar waveform at the range equal to the radar's maximum detection range is close to the desired value 0.95.

Frequency Agility

To avoid interception, the radar can alter its center frequency so that most or some of the waveform's power falls outside the interceptor's band. This strategy prevents the selected integration window from accumulating sufficient energy to reach the desired probability of detection.

Assume that the radar center frequency hopped from 9.375 GHz to 9.330 GHz.

Radar.NewOperatingFrequency = 9.320e9;

Compute the new frequency offset between the radar and the intercept receiver.

% Frequency offset between the radar and the intercept receiver
frequencyOffset = Radar.NewOperatingFrequency - Interceptor.OperatingFrequency;

pinband = helperWithinInterceptorBandwidth(Radar.NewOperatingFrequency, Radar.Bandwidth, Interceptor.OperatingFrequency, Interceptor.Bandwidth);
fprintf('%.2f %% of radar waveform bandwidth is within the intercept receiver band.\n',pinband);

60.00 % of radar waveform bandwidth is within the intercept receiver band.

Regenerate the signals received at the intercept receiver giving the new frequency offset.

interceptedSignal = collectedAtInterceptorSignal.*exp(1i*2*pi*frequencyOffset.*t);
interceptedSignal = lowpass(interceptedSignal, Interceptor.Bandwidth/2, sampleRate);
interceptedSignal = resample(interceptedSignal, Interceptor.Bandwidth, sampleRate);

receivedInterceptedSignal = zeros(size(interceptedSignal, 1), numInterceptorTrials, numInterceptorRanges);

for it = 1:numInterceptorTrials
    receivedInterceptedSignal(:, it, :) = Interceptor.Receiver(interceptedSignal);
end

Estimate the probability of intercepting the radar waveform and plot it as a function of the range to the interceptor.

estimatedPd2 = helperEstimateInterceptorPd(receivedInterceptedSignal, N, M, interceptorDetectionThreshold);

figure;
hold on;
plot(interceptorRange*1e-3, pdInBand, 's--', LineWidth=2, MarkerFaceColor='auto');
plot(interceptorRange*1e-3, estimatedPd2, '^--', LineWidth=2, MarkerFaceColor='auto');
grid on;
xlabel('Range to interceptor (km)');
ylabel('P_d');
title('Probability of Intercepting Radar Waveform')
legend("100% radar power intercepted", sprintf('%.0f%% radar power intercepted\n (frequency hopping)', pinband), Location="southoutside");
ylim([0 1.1]);

xline(Rmax*1e-3, '--', 'Rmax', 'LabelHorizontalAlignment', 'center', 'LabelVerticalAlignment', 'top', 'HandleVisibility','off');
xline(RiMN*1e-3, '--', 'Ri', 'LabelHorizontalAlignment', 'center', 'LabelVerticalAlignment', 'bottom', 'HandleVisibility','off');
yline(Interceptor.Pd, '--', 'Required P_d', 'HandleVisibility','off');

Figure contains an axes object. The axes object with title Probability of Intercepting Radar Waveform, xlabel Range to interceptor (km), ylabel P indexOf d baseline P_d contains 2 objects of type line. These objects represent 100% radar power intercepted, 60% radar power intercepted (frequency hopping).

When a significant portion of the waveform's power falls outside the interceptor's bandwidth, the probability of intercepting the radar waveform is only achieved at much shorter ranges within the radar's maximum range. As a result, the radar system effectively functions as an LPI system. This illustrates that the LPI status is contingent upon the parameters of both the radar and the intercept receiver.

Conclusion

This example demonstrates how to evaluate the performance of a radar system and a non-cooperative intercept receiver to determine the radar's LPI status. It starts with a power-level analysis that calculates the processing gain and corresponding noncoherent integration time needed for the intercept receiver to detect the radar waveform beyond the radar's maximum range. Additionally, the example simulates I/Q signals of both the radar and the intercept receiver to confirm that the desired probability of intercepting the radar waveform is achieved. This illustrates that if the intercept receiver employs energy detection and a sufficient portion of the radar waveform's power resides within the interceptor's bandwidth, interception is always possible regardless of the waveform type. A viable strategy for the radar to evade detection is to alter its center frequency so that at least some of its transmitted power falls outside the intercept receiver's bandwidth.

References

Pace, Phillip E. Detecting and classifying low probability of intercept radar. Artech house, 2009.
Schleher, D. C. "LPI radar: Fact or fiction." IEEE Aerospace and Electronic Systems Magazine 21, no. 5 (2006): 3-6.
Richards, Mark A. Fundamentals of radar signal processing. Vol. 1. New York: Mcgraw-hill, 2005.
Schrick, Gerd, and Richard G. Wiley. "Interception of LPI radar signals." In IEEE international conference on radar, pp. 108-111. IEEE, 1990.

Supporting Functions

function rsa = helperReceiverSensitivityRatio(SNRi, Dx, Bi, Br, betai, betar, Tsi, Tsr)
    rsa = Tsi*Bi*(db2pow(SNRi)/betai)/(Tsr*Br*(db2pow(Dx)/betar));
end

function pdest = helperEstimateInterceptorPd(X, N, M, threshold)
    [numSamples, numTrials, numRanges] = size(X);
    K = floor(numSamples/N);
    y = zeros(K-1, numTrials, numRanges);

    for ir = 1:numRanges
        for k = 2:K
            m = sum(abs(X((k-1)*N+1:k*N, :, ir)).^2 > threshold);
            y(k-1, :, ir) = m >= M;
        end
    end
    
    y = reshape(y, (K-1)*numTrials, numRanges);
    pdest = mean(y);
end

function p = helperWithinInterceptorBandwidth(fcr, Br, fci, Bi)
    fminr = fcr - Br/2;
    fmaxr = fcr + Br/2;
    
    fmini = fci - Bi/2;
    fmaxi = fci + Bi/2;

    if fmaxr < fmini || fminr > fmaxi
        p = 0;
    else
        fmax = min(fmaxr, fmaxi);
        fmin = max(fminr, fmini);
        p = (fmax-fmin)/Br * 100;
    end
end