Interference from Satellite Constellation on Communications Link

This example shows how to analyze interference on a downlink from a constellation of satellites. The downlink occurs from a geosynchronous satellite to a ground station located in the Pacific Ocean. The interfering constellation consists of 40 satellites in a low-Earth orbit. This example determines the times at which the downlink is closed, the carrier to noise plus interference ratio, and the link margin.

Create Satellite Scenario

Create a satellite scenario. Define the start time and stop time of the scenario. Set the sample time to 60 seconds.

startTime = datetime(2021,3,17);                       % 17 March 2021 12:00 AM UTC
stopTime = startTime + days(1);                        % 18 March 2021 12:00 AM UTC
sampleTime = 60;                                       % In s
sc = satelliteScenario(startTime,stopTime,sampleTime);

Add a satellite in a geosynchronous orbit by specifying its Keplerian orbital elements. This satellite is the satellite from which data is downlinked.

earthAngularVelocity = 0.0000729211585530;               % In rad/s
orbitalPeriod = 2*pi/earthAngularVelocity;               % In s
earthStandardGravitationalParameter = 398600.4418e9;     % In m^3/s^2
semiMajorAxis = (earthStandardGravitationalParameter ...
/(earthAngularVelocity^2))^(1/3);                    % In m
eccentricity = 0;
inclination = 8;                                         % In degrees
raan = 0;                                                % In degrees
argOfPeriapsis = 0;                                      % In degrees
trueAnomaly = 0;                                         % In degrees
geoSat = satellite(sc, semiMajorAxis, ...
eccentricity, ...
inclination, ...
raan, ...
argOfPeriapsis, ...
trueAnomaly, ...
"Name","Geo Sat");

Add the interfering satellite constellation from a two-line-element (TLE) file. Declutter the visualizations that you launch later in this example by setting the ShowLabel property of this satellite constellation to false.

interferingSat = satellite(sc,"leoSatelliteConstellation.tle");
interferingSat.ShowLabel = false;

Launch Satellite Scenario Viewer

Launch the Satellite Scenario Viewer. Hide the orbits of the interfering satellites to further declutter the visualizations.

v = satelliteScenarioViewer(sc);
hide(interferingSat.Orbit); Add Transmitter to Geosynchronous Orbit Satellite

Add a transmitter to the geosynchronous orbit satellite. This transmitter is used for the downlink. Define the antenna specifications and set the operating carrier frequency to 3 GHz.

txGeoFreq = 3e9;                   % In Hz
txGeoSat = transmitter(geoSat, ...
"Frequency",txGeoFreq, ...     % In Hz
"Power",25);                   % In dBW
gaussianAntenna(txGeoSat, ...
"DishDiameter",1);             % In m

Add a transmitter to each satellite in the interfering constellation, and then define the antenna specifications. These transmitters are the ones that interfere with the downlink from the geosynchronous orbit satellite. Set the operating carrier frequency of the interfering satellites to 2.99 GHz. The example assigns each interfering satellite a random power in the range from 10 to 20 dBW.

interferenceFreq = 2.99e9;                                  % In Hz
rng("default");
for idx = 1:numel(interferingSat)
txInterferingSat = transmitter(interferingSat(idx), ...
"Frequency",interferenceFreq, ...                   % In Hz
"Power",10+10*rand);                                % In dBW
gaussianAntenna(txInterferingSat, ...
"DishDiameter",0.2);                                % In m
end

Add a ground station to satellite scenario by specifying its latitude and longitude.

gs = groundStation(sc, ...
0, ...                    % Latitude in degrees
180, ...                  % Longitude in degrees
"Name","Ground station");

Add a gimbal to the ground station by specifying its mounting location and angles. This gimbal is used to make the receiver antenna of the ground station track the geosynchronous orbit satellite.

gim = gimbal(gs, ...
"MountingLocation",[0;0;-5], ... % In m
"MountingAngles",[0;180;0]);     % In m

Add a receiver to the ground station gimbal, specifying the receiver's mounting location with respect to the gimbal. Define the specifications of the receiver antenna.

"MountingLocation",[0;0;1]); % In m
gaussianAntenna(rxGs, ...
"DishDiameter",0.8);         % In m

Set Tracking Targets for Satellites and Gimbal

Set the satellites to track the ground station, and the ground station gimbal to track the geosynchronous orbit satellite. This ensures that the transmitter antennas on board each satellite track the ground station and that the ground station antenna tracks the geosynchronous orbit satellite. Setting the interfering satellite transmitters to track the ground station results in the worst-case interference on the downlink.

pointAt(geoSat,gs);
pointAt(gim,geoSat);
for idx = 1:numel(interferingSat)
pointAt(interferingSat(idx),gs);
end

Visualize the radiation pattern of the transmitter antenna on board the geosynchronous orbit satellite and the receiver on board the ground station.

pattern(geoSat.Transmitters, ...
'Size',10000000);                           % In m
pattern(rxGs,geoSat.Transmitters.Frequency, ...
'Size',1000000);                            % In m Create a downlink from the transmitter on board the geosynchronous orbit satellite to the receiver on board the ground station. This link is the downlink which encounters interference from the satellite constellation.

Create a link between the transmitter on board each interfering satellite in the constellation and the receiver on board the ground station. These links are the interferer links with the desired downlink.

for idx = 1:numel(interferingSat)
end

Create Access Analysis Between Interfering Satellite Constellation and Ground Station

Add an access analysis between each satellite in the interfering constellation and the ground station. This analysis enables the visualization of interference. Any time a satellite in the constellation is visible to the ground station, there is some level of interference from that visible satellite.

for idx = 1:numel(interferingSat)
ac = access(interferingSat(idx),gs);
ac.LineColor = [1 1 0];              % Yellow
end

The presence of the green line between the transmitter on board the geosynchronous satellite and the receiver on board the ground station signifies that the downlink can be closed successfully assuming no interference from the satellite constellation exists. The presence of yellow lines between a given satellite in the constellation and the ground station signifies that an interference from that satellite exists. Play Scenario

Play the scenario in the Satellite Scenario Viewer. Set the target of the camera to the ground station. Place the camera at a -10 degree latitude, 170 degree longitude, and 4,000 km altitude. Additionally, set the camera heading and pitch angles to 40 and -60 degrees, respectively.

play(sc);
camtarget(v,gs);
campos(v,-10,170,4000000);
campitch(v,-60); Plot Downlink Closure Status Neglecting Interference

Determine the closure status of the desired downlink from the geosynchronous orbit satellite. The linkStatus function neglects interference from other transmitters. Any time the downlink is close, the status is true. Otherwise, the status is false. The status is indicated by 1 and 0, respectively in the plot. Because the satellite involved in the downlink is in a geosynchronous orbit above the ground station, the downlink is closed for the duration of the scenario.

xlabel("Time");
title("Link Status as a Function of Time");
grid on; Although the linkStatus function neglects interference when determining the times when the downlink is closed, manually calculate the interference based on data that can be extracted from the scenario is possible. To do this manual calculation, understanding how linkStatus calculates the link closure, as shown in these steps, is necessary.

1) Calculate the effective isotropic radiated power (EIRP):

$EIR{P}_{Tx}={P}_{Tx}-LOS{S}_{Tx}+{G}_{TxAntenna}$,

where:

• $EIR{P}_{Tx}$ is the EIRP of the transmitter antenna (in dBW).

• ${P}_{Tx}$ is the transmitter power (in dBW).

• $LOS{S}_{Tx}$ is the total system loss of the transmitter (in dB).

• ${G}_{TxAntenna}$ is the gain of the transmitter antenna in the direction of the receiver antenna (in dB).

2) In the satellite scenario, all path loss computations assume a free space propagation model. Calculate the free space path loss from the transmitter to the receiver antenna as

$FSPL=10{\mathrm{log}}_{10}\left({\left(\frac{4\pi df}{c}\right)}^{2}\right)$,

where:

• $FSPL$ is the free space path loss from the transmitter to the receiver antenna (in dB).

• $d$ is the distance between the transmitter and the receiver antenna (in m).

• $f$ is the transmitter frequency (in Hz).

• $c$ is the speed of light in vacuum (in m/s).

3) Calculate the received isotropic power as

$RI{P}_{Rx}=EIR{P}_{Tx}-FSPL$,

where $RI{P}_{Rx}$ is the received isotropic power at the received antenna (in dBW).

4) Calculate the carrier to noise density ratio as

$C/{N}_{0}=RI{P}_{Rx}+{\left(G/T\right)}_{Rx}-10{\mathrm{log}}_{10}{k}_{B}-LOS{S}_{RxSystem}$,

where:

• $C/{N}_{0}$ is the carrier to noise density ratio (in dB).

• ${\left(G/T\right)}_{Rx}$ is the gain to noise temperature ratio of the receiver antenna (in dB/K).

• ${k}_{B}$ is the Boltzmann constant (in J/K).

• $LOS{S}_{Rx}$ is the total system loss of the transmitter.

The receiver antenna gain is computed in the direction of the transmitter antenna. The noise temperature is assumed to be a constant and is derived from the GainToNoiseTemperatureRatio property of the receiver, which in turn is based on the receiver antenna gain along the $z$-axis of the receiver. In essence:

${\left(G/T\right)}_{z-axis}={\left({G}_{Rx}\right)}_{z-axis}-10{\mathrm{log}}_{10}T$,

where:

• ${\left(G/T\right)}_{z-axis}$ is the receiver antenna gain to noise temperature ratio along the $z$-axis of the receiver (in dB).

• ${\left({G}_{Rx}\right)}_{z-axis}$ is the receiver antenna gain in the direction of the $z$-axis of the receiver (in dB).

• $T$ is the noise temperature (in K).

Rearranging the above equation results in

$10{\mathrm{log}}_{10}T={\left({G}_{Rx}\right)}_{z-axis}-{\left(G/T\right)}_{z-axis}$.

Consequently,

${\left(G/T\right)}_{Rx}=\left({G}_{Rx}\right)-10{\mathrm{log}}_{10}T$.

5) Calculate the receiver energy per bit to noise power spectral density ratio as

${E}_{b}/{N}_{0}=C/{N}_{0}-10{\mathrm{log}}_{10}\left(BITRATE\right)-60$,

where:

• ${E}_{b}/{N}_{0}$ is the received energy per bit to noise power density ratio (in dB).

• $BITRATE$ is the bit rate of the link (in Mbps).

The 60 value appears in the equation because the bit rate is in Mbps.

6) Calculate the link margin as

$MARGIN={E}_{b}/{N}_{0}-{\left({E}_{b}/{N}_{0}\right)}_{Required}$,

where:

• $MARGIN$ is the link margin (in dB).

• ${\left({E}_{b}/{N}_{0}\right)}_{Required}$ is the minimum received energy per bit to noise power spectral density ratio that is required to close the link (in dB).

The link is closed when the link margin is greater than or equal to 0 dB.

To account for interference, the signal received at the ground station antenna from each interfering transmitter must be treated as noise and must be added to the noise power. To compute the new received Eb/No and link margin, the received signal power from each transmitter in the scenario must be calculated.

The received power after antenna is the quantity $RI{P}_{Rx}+{G}_{RxAntenna}$. This received power can be obtained by starting from ${E}_{b}/{N}_{0}$ and using the equations in the previous section backward. First, the received signal power corresponding to the transmitter on board the geosynchronous orbit satellite must be computed. Retrieve the received ${E}_{b}/{N}_{0}$ corresponding to this transmitter.

Use these equations to calculate the corresponding received signal power after antenna.

$C/{N}_{0}={E}_{b}/{N}_{0}+10{\mathrm{log}}_{10}\left(BITRATE\right)+60$

$RI{P}_{Rx}+{G}_{RxAntenna}=C/{N}_{0}+10{\mathrm{log}}_{10}{k}_{B}+10{\mathrm{log}}_{10}T+LOS{S}_{Rx}$

T = HelperGetNoiseTemperature(txGeoSat.Frequency,rxGs);                        % In K
kb = physconst("Boltzmann");
downlinkSigPower = CNoDownlink + 10*log10(kb) + 10*log10(T) + rxGs.SystemLoss; % In dBW

Calculate Total Interfering Signal Power After Receiver Antenna from All Interfering Transmitters

The total interfering signal power after the receiver antenna from all interfering transmitters is calculated using these steps.

1) Calculate $RI{P}_{Rx}+{G}_{Rx}$, which corresponds to each interfering transmitter, using the same equations used for the transmitter from the geosynchronous orbit satellite.

2) Convert the quantities to Watts.

3) Add the quantities in Watts and the final sum is the total interfering power, ${I}_{W}$, in Watts.

% Initialize an array to store the received Eb/No history corresponding to
% each interfering transmitter
ebnoInterference = zeros(numel(lnkInterference), numel(t));

% Retrieve the Eb/No history corresponding to each interfering transmitter
% and store it in the array
for idx = 1:numel(lnkInterference)
ebnoInterference(idx,:) = ebno(lnkInterference(idx));
end

% Calculate the interfering signal power corresponding to each interfering
% transmitter in Watts (in this example, the bit rate of each interfering
% transmitter is the same as that of that of the geosynchronous orbit
% satellite)
bitRate = txGeoSat.BitRate;
CNoInterference = ebnoInterference + 10*log10(bitRate) + 60;
interferenceSigPower = CNoInterference + 10*log10(kb) ...
+ 10*log10(T) + rxGs.SystemLoss;
interferenceSigPowerW = 10.^(interferenceSigPower/10);       % In W

% Add the interfering signal power from each transmitter
interferenceSigPowerSumW = sum(interferenceSigPowerW);       % In W

Calculate Contribution of Interfering Signal Power in Overlapped Bandwidth

Calculate the amount of total interfering signal power that contributes to interference in the signal bandwidth by following these steps.

1) Calculate the overlapping portion of the signal bandwidth with the bandwidth of the interferers. This example considers the transmission power of interfering satellites and geosynchronous satellite as constant across the whole bandwidth of respective geosynchronous satellite and interfering satellites.

2) Calculate the amount of interference power that acts as interference to signal bandwidth.

This diagram shows the power spectral density (PSD) plot, which shows the actual interference power and modeled interference power when the transmission bandwidth and interfering bandwidth overlap. The actual interference power is the area occupied by the interference power density in the overlapped bandwidth region. This actual interference power is then spread across the entire transmission bandwidth and assumed to be noise-like. This example assumes that the transmission (or signal) bandwidth of the geosynchronous satellite is 30 MHz and that the bandwidth of the interfering signal is 20 MHz.

txBandwidth = 30e6;                                               % In Hz
interferenceBandWidth = 20e6;                                     % In Hz

% Get the overlap portion of the interfering bandwidth and the bandwidth of
% interest. The assumption is to have same power across the whole
% bandwidth.
overlapFactor = getOverlapFactor(txGeoFreq,txBandwidth, ...
interferenceFreq,interferenceBandWidth);

% Get the interference power that contributes as interference to the signal
% of interest from the total interference power
interferenceSigPowActual = interferenceSigPowerSumW*overlapFactor; % In W

Calculate Downlink Closure Status That Accounts for Interference

Calculate the downlink closure status that accounts for interference by following these steps.

1) Calculate the carrier to noise plus interference power density ratio as

$C/\left({N}_{0}+{I}_{0}\right)=RI{P}_{Rx}+{G}_{Rx}-10{\mathrm{log}}_{10}\left(\frac{{I}_{W}}{TxBandwidth}+{k}_{B}T\right)-LOS{S}_{Rx}$,

where:

• $C/\left({N}_{0}+{I}_{0}\right)$ is the carrier to noise plus interference power density ratio (in dB).

• ${I}_{W}$ is the interference signal power after receiver antenna that interferes with the signal bandwidth (in W).

• $TxBandwidth$ is the downlink transmission bandwidth from geosynchronous satellite (in Hz).

2) Calculate the energy per bit to noise plus interference power spectral density ratio as

${E}_{b}/\left({N}_{0}+{I}_{0}\right)=C/\left({N}_{0}+{I}_{0}\right)-10{\mathrm{log}}_{10}\left(BITRATE\right)-60$,

where ${E}_{b}/\left({N}_{0}+{I}_{0}\right)$ is the energy per bit to noise plus interference power spectral density ratio (in dB).

3) Calculate the link margin accounting for interference:

$MARGIN={E}_{b}/\left({N}_{0}+{I}_{0}\right)-{\left({E}_{b}/{N}_{0}\right)}_{Required}$.

% Calculate link status in the presence of interference
10*log10((interferenceSigPowActual/txBandwidth) + T*kb) - rxGs.SystemLoss;
ebNoPlusInterference = CNoPlusInterference - 10*log10(bitRate) - 60;
marginWithInterference = ebNoPlusInterference - rxGs.RequiredEbNo;

Plot Downlink Closure Status Accounting for Interference

Plot the new downlink closure status that accounts for interference. Compare the new link status with the previous case when interference was neglected.

legend("Interference accounted","Interference neglected");
xlabel("Time");
title("Link Status as a Function of Time");
ylim([0 1.2]);
grid on The plot shows times exist when the downlink cannot be closed because of interference. The interference is particularly severe when at least one satellite flies overhead. For instance, the plot shows that the downlink is broken at 3:19 AM, 6:49 AM, and 10:54 PM. Set the CurrentTime property of the Satellite Scenario Viewer to these times to confirm that an interfering satellite overflying the ground station exists.

v.CurrentTime = datetime(2021,3,17,3,19,0); v.CurrentTime = datetime(2021,3,17,6,49,0); v.CurrentTime = datetime(2021,3,17,22,54,0); Calculate Carrier to Noise Ratio and Carrier to Noise Plus Interference Ratio

Calculate the carrier to noise ratio (CNR) and carrier to noise plus interference ratio (CNIR) from the carrier to noise density ratio and carrier to noise plus interference power density as:

$C/N=C/{N}_{0}-10{\mathrm{log}}_{10}\left(TxBandwidth\right)$ and

$C/\left(N+I\right)=C/\left({N}_{0}+{I}_{0}\right)-10{\mathrm{log}}_{10}\left(TxBandwidth\right)$,

where:

• $C/N$ is the carrier to noise ratio

• $C/\left(N+I\right)$ is the carrier to noise plus interference ratio.

• $TxBandwidth$ is the downlink transmission bandwidth from geosynchronous satellite (in Hz).

cByNPlusI = CNoPlusInterference - 10*log10(txBandwidth);

Plot $C/N$ and $C/\left(N+I\right)$.

plot(t,cByNPlusI,"-r",t,cByN,"--g","LineWidth",2);
legend("C/(N+I)", "C/N","Location","south");
xlabel("Time");
ylabel("C/N or C/(N+I) (dB)");
title("CNR and CNIR as a Function of Time");
grid on The downward spikes indicate interference. The more severe spikes are caused when simultaneous interference from multiple satellites exists or when an interfering satellite overflies the ground station.

Compare Link Margins with and without Interference

Calculate the link margin without interference.

Plot the link margins with and without interference.

plot(t,marginWithInterference,"-r",t,marginWithoutInterference,"--g","LineWidth",2);
legend("With interference","Without interference","Location","south");
xlabel("Time");
ylabel("Margin (dB)");
title("Link Margin as a Function of Time");
grid on Any time the link margin is greater than or equal to 0 dB, the downlink is closed. With interference, times exist when the link margin dips below 0 dB and at these times the downlink is broken because of interference.

Further Exploration

This example demonstrates how to analyze interference on a satellite communication link. The link closure times are a function of these parameters:

• The orbit of the satellites

• The position of the ground station

• The specifications of the transmitters and the receiver

• The signal and interference bandwidth

Modify these parameters to observe their influence on the level of interference on the link.

Helper Functions

The example uses the helper function HelperGetNoiseTemperature to obtain the noise temperature of the receiver antenna.

The example also uses this local function to compute the amount of overlap between the transmission bandwidth and the interfering bandwidth.

function overlapFactor = getOverlapFactor(txFreq,txBW,interferenceFreq,interferenceBW)
% getOverlapFactor provides the amount of interference bandwidth overlapped
% with transmission bandwidth

txFreq_Limits = [txFreq-(txBW/2) txFreq+(txBW/2)];
interferenceFreq_Limits = [interferenceFreq-(interferenceBW/2) ...
interferenceFreq+(interferenceBW/2)];
if (interferenceFreq_Limits(2) < txFreq_Limits(1)) || ...
(interferenceFreq_Limits(1) > txFreq_Limits(2))
% If no overlap exists between transmission bandwidth and
% interfering bandwidth, then overlap factor is 0
overlapFactor = 0;
elseif (interferenceFreq_Limits(2) <= txFreq_Limits(2)) && ...
(interferenceFreq_Limits(1) >= txFreq_Limits(1))
% If interfering bandwidth lies completely within transmission
% bandwidth, then overlap factor is 1
overlapFactor = 1;
elseif (interferenceFreq_Limits(2) > txFreq_Limits(2)) && ...
(interferenceFreq_Limits(1) < txFreq_Limits(1))
% If transmission bandwidth lies completely within interfering
% bandwidth, then overlap factor is the ratio of transmission
% bandwidth with that of interference bandwidth
overlapFactor = txBW/interferenceBW;
elseif (interferenceFreq_Limits(2) <= txFreq_Limits(2)) && ...
(interferenceFreq_Limits(1) <= txFreq_Limits(1))
% If the start edge of transmission bandwidth lies within
% interfering bandwidth, then overlap factor is the ratio of
% difference from last edge of interfering bandwidth and first edge
% of signal bandwidth, with that of interference bandwidth
overlapFactor = (interferenceFreq_Limits(2)-txFreq_Limits(1))/interferenceBW;
else
% If the last edge of transmission bandwidth lies within
% interfering bandwidth, then overlap factor is the ratio of difference
% from last edge of signal bandwidth and first edge of interfering
% bandwidth, with that of interference bandwidth
overlapFactor = (-interferenceFreq_Limits(1)+txFreq_Limits(2))/interferenceBW;
end

end

Related Examples 