interference from satellite constellation on communications link -凯发k8网页登录
this example shows how to analyze interference on a downlink from a constellation of satellites in a medium-earth orbit (meo) to a ground station located in the pacific ocean. the interfering constellation consists of 40 satellites in a low-earth orbit (leo). this example determines the times at which the downlink is closed, the carrier to noise plus interference ratio, and the link margin.
this example requires satellite communications toolbox. if you have also have antenna toolbox™, you can use this example to learn how to import antennas from antenna toolbox into satellite scenario. if you also have phased array system toolbox™, you can use this example to learn how to import antennas from phased array system toolbox into satellite scenario and use beamforming to improve the carrier-to-noise-and-interference ratio (cnir).
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,22,52,0); % 17 march 2021 10:52 pm utc stoptime = starttime minutes(10); % 17 march 2021 11:02 pm utc sampletime = 60; % in s sc = satellitescenario(starttime,stoptime,sampletime);
add medium-earth orbit satellite
add a satellite in an meo by specifying its keplerian orbital elements. this satellite is the satellite from which data is downlinked.
semimajoraxis = 12000000; % in m eccentricity = 0; inclination = 8; % in degrees raan = 0; % right ascension of ascending node, in degrees argofperiapsis = 0; % in degrees trueanomaly = 343.9391; % in degrees meosat = satellite(sc, semimajoraxis, ... eccentricity, ... inclination, ... raan, ... argofperiapsis, ... trueanomaly, ... name = "meo satellite", ... orbitpropagator = "two-body-keplerian");
add interfering satellite constellation
add the interfering satellite constellation from a two-line-element (tle) file. these satellites are placed in leo.
interferingsat = satellite(sc,"leosatelliteconstellation.tle");
add transmitter to meo satellite
add a transmitter to the meo satellite. this transmitter is used for the downlink. define the antenna specifications and set the operating carrier frequency to 3 ghz.
txmeofreq = 3e9; % in hz txmeosat = transmitter(meosat, ... frequency = txmeofreq, ... % in hz power = 11); % in dbw gaussianantenna(txmeosat, ... dishdiameter = 1); % in m
add transmitter to leo satellites
add a transmitter to each satellite in the leo constellation, and then define the antenna specifications. these transmitters are the ones that interfere with the downlink from the meo 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"); txinterferingsat = transmitter(interferingsat, ... frequency = interferencefreq, ... % in hz power = 10 10*rand(1,numel(interferingsat))); % in dbw gaussianantenna(txinterferingsat, ... dishdiameter = 0.2); % in m
add ground station
add a ground station to the satellite scenario by specifying its latitude and longitude.
gs = groundstation(sc, ... 0, ... % latitude in degrees 180, ... % longitude in degrees name = "ground station");
specify ground station antenna type
for this example, you can choose from one of the following antennas:
gaussian antenna
parabolic reflector from antenna toolbox
uniform rectangular array from phased array system toolbox
% select the desired ground station antenna. groundstationantennatype = "gaussian antenna";
add receiver to ground station
add a receiver to the ground station. if you selected gaussian antenna or parabolic reflector, attach the receiver to a gimbal, which is in turn attached to the ground station. configure the gimbal to track the meo satellite, so that the antenna also tracks the meo satellite. if you selected uniform rectangular array, attach the receiver directly to the ground station. specify the mounting location and mounting angles of the gimbal and receiver, and the antenna specifications appropriately.
switch groundstationantennatype case {"gaussian antenna","parabolic reflector"} % when gaussian antenna or parabolic reflector is selected, attach % a gimbal to the ground station. gim = gimbal(gs, ... mountinglocation = [0;0;-5], ... % in m mountingangles = [0;180;0]); % in degrees % set the gimbal to track the meo satellite. pointat(gim,meosat); if groundstationantennatype == "gaussian antenna" % when gaussian antenna is selected % create the receiver object and add it to the gimbal rxgs = receiver(gim, ... mountinglocation = [0;0;1]); % in m % provide the gaussian antenna specifications gaussianantenna(rxgs, ... dishdiameter = 0.8); % in m else % when parabolic reflector is selected % size the antenna based on the frequency % requires antenna toolbox(tm) ant = design(reflectorparabolic,txmeofreq); % create the receiver object and add it to the gimbal rxgs = receiver(gim, ... antenna = ant, ... mountinglocation = [0;0;1]); % in m end case "uniform rectangular array" % when uniform rectangular array is selected % determine the wavelength of the downlink signal c = physconst('lightspeed'); lambda = c/txmeofreq; % define array size nrow = 8; ncol = 8; % define element spacing drow = lambda/2; dcol = lambda/2; % create a back-baffled 6-by-6 antenna array % requires phased array system toolbox(tm) ant = phased.ura(size = [nrow ncol], ... elementspacing = [drow dcol]); ant.element.backbaffled = true; % create the receiver object and add it to the ground station rxgs = receiver(gs, ... antenna = ant, ... mountingangles = [0;90;0]); % in degrees 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 in the satellite scenario viewer that will be launched later. any time a satellite in the constellation is visible to the ground station, there is some level of interference from that visible satellite.
ac = access(interferingsat,gs);
ac.linecolor = [1 1 0]; % yellow
set tracking targets for satellites
set the satellites to track the ground station. this ensures that the transmitter antennas on board each satellite track the ground station. setting the interfering satellite transmitters to track the ground station results in the worst-case interference on the downlink.
pointat([meosat interferingsat],gs);
calculate weights of uniform rectangular array
if you selected uniform rectangular array as the ground station antenna, compute the weights that are required to point the main lobe towards the meo satellite, and the nulls towards the interfering satellites, thereby cancelling the interference. assign the computed weights using pointat
.
if groundstationantennatype == "uniform rectangular array" % find the leo satellites that are in the line of sight of the ground % station. these satellites are the potential interferers. currentinterferingsat = interferingsat(accessstatus(ac,sc.starttime) == true); % calculate the direction of the meo satellite with respect to the % array. this is the lookout direction. [azd,eld] = aer(rxgs,meosat,sc.starttime,coordinateframe='body'); % calculate the directions of the potentially interfering satellites % with respect to the array. these are the null directions. [azn,eln] = aer(rxgs,currentinterferingsat,sc.starttime,coordinateframe='body'); % calculate the steering vectors for the lookout direction. % requires phased array system toolbox. wd = steervec(getelementposition(ant)/lambda,[wrapto180(azd);-eld]); % calculate the steering vector for null directions. % requires phased array system toolbox. wn = steervec(getelementposition(ant)/lambda,[wrapto180(azn)';-eln']); % compute the response of the desired steering at null directions. rn = (wn'*wn)\(wn'*wd); % sidelobe canceler - remove the response at null directions. w = wd-wn*rn; % assign the weights to the phased array. pointat(rxgs,weights=w); end
create desired downlink
create a downlink from the transmitter on board the meo satellite to the receiver on board the ground station. this link is the downlink which encounters interference from the leo constellation.
downlink = link(txmeosat,rxgs);
create interfering links
create a link between the transmitter on board each satellite in the leo constellation and the receiver on board the ground station. these links are the interferer links with the desired downlink.
lnkinterference = link(txinterferingsat,rxgs);
launch satellite scenario viewer
launch the satellite scenario viewer with showdetails
set to false. when the showdetails
property is set to false
, only the satellites, the ground station, accesses, and links will be shown. the labels and orbits will be hidden. mouse over the satellites and the ground stations to show their labels. click on the meo satellite so that its orbit projected up to the scenario stoptime
and its label are visible without mousing over. click on the ground station so that its label is visible without mousing over. the presence of the green line between the transmitter on board the meo 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 they have access to one another, and as a result, interference from that satellite exists.
v = satellitescenarioviewer(sc,showdetails=false);
visualize radiation pattern of antennas involved in downlink
visualize the radiation pattern of the transmitter antenna on board the meo satellite and the receiver on board the ground station.
pattern(txmeosat, ... size = 1000000); % in m pattern(rxgs,txmeofreq, ... size = 1000000); % in m
radiation pattern of meo satellite antenna
set camera to view ground station antenna radiation pattern
% set camera position and orientation to view the ground station antenna % radiation pattern. campos(v,-8,172,2500000); camheading(v,40); campitch(v,-60);
gaussian antenna
parabolic reflector
uniform rectangular array
simulate scenario and visualize
with gaussian antenna or parabolic reflector
if you selected gaussian antenna or parabolic reflector (requires antenna toolbox), use play
to visualize the scenario from starttime
to stoptime
. this will automatically simulate the scenario before playing back the visualization. note how the antenna pointing changes as the gimbal tracks the meo satellite.
if groundstationantennatype == "gaussian antenna" || groundstationantennatype == "parabolic reflector" play(sc); campos(v,-8,172,2500000); camheading(v,40); campitch(v,-60); end
with gaussian antenna
with parabolic reflector
with uniform rectangular array
if you selected uniform rectangular array (requires phased array system toolbox), you must manually step through the simulation so that you can recompute the weights at each time step based on the new position of the meo satellite and the interfering leo satellites. to manually step through the simulation, first set autosimulate
to false. following this, you can call advance
to move the simulation by one time step. the first call to advance
will compute the simulation states at starttime
. subsequent calls will advance the time step by one sampletime
and compute the states accordingly.
if groundstationantennatype == "uniform rectangular array" % set autosimulate to false. sc.autosimulate = false; % manually step through the simulation. while advance(sc) % determine the access status history for each leo satellite % corresponding to the current simulationtime. acstatushistory = accessstatus(ac); acstatus = acstatushistory(:,end); % determine the leo satellites that are visible to the ground % station. these are the satellites that will potentially % interfere with the ground station at the current simulation % time. currentinterferingsat = interferingsat(acstatus == true); % determine the direction of the meo satellite in the body frame of % the uniform rectangular array. this is the lookout direction of % the array. [azdhistory,eldhistory] = aer(rxgs,meosat,coordinateframe='body'); azd = azdhistory(:,end); eld = eldhistory(:,end); % determine the direction of these interfering satellites in % the body frame of the uniform rectangular array. these are % the directions in which the array must point a null. [aznhistory,elnhistory] = aer(rxgs,currentinterferingsat,coordinateframe='body'); azn = aznhistory(:,end); eln = elnhistory(:,end); % calculate the steering vectors for lookout direction. % requires phased array system toolbox. wd = steervec(getelementposition(ant)/lambda,[wrapto180(azd);-eld]); % calculate the steering vector for null directions. % requires phased array system toolbox. wn = steervec(getelementposition(ant)/lambda,[wrapto180(azn)';-eln']); % compute the response of desired steering at null direction. rn = (wn'*wn)\(wn'*wd); % sidelobe canceler - remove the response at null direction. w = wd-wn*rn; % assign the weights to the phased array. pointat(rxgs,weights=w); end end
plot downlink closure status neglecting interference
determine the closure status of the desired downlink from the meo satellite. the linkstatus
function neglects interference from other transmitters. any time the downlink is closed, the status is true. otherwise, the status is false. the status is indicated by 1 and 0, respectively in the plot.
[downlinkstatus,t] = linkstatus(downlink); plot(t,downlinkstatus,"-g",linewidth=2); xlabel("time"); ylabel("downlink closure status"); title("link status as a function of time"); grid on;
calculate downlink closure status with interference
calculate the downlink closure status with interference by first calculating the meo downlink and interference signal power levels at the ground station receiver input using sigstrength.
the locations of received power measurements and losses are illustrated in the ground station receiver diagram below.
% calculate the power at receiver input corresponding to the downlink from % the meo satellite. [~,downlinkpowerrxinput] = sigstrength(downlink); % in dbw % calculate the interference power at receiver input corresponding to each % leo satellite. [~,interferencepowerrxinput] = sigstrength(lnkinterference); % in dbw
calculate total interfering signal power at the receiver input. get this quantity by summing the individual power levels from the interfering leo satellites in watts.
interferencepowerrxinputw = 10.^(interferencepowerrxinput/10); % w interferencepowerrxinputsumw = sum(interferencepowerrxinputw); % w
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 the meo satellite as constant across the whole bandwidth of respective meo 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 meo 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 the have same interference power across % the whole bandwidth. overlapfactor = getoverlapfactor(txmeofreq,txbandwidth, ... interferencefreq,interferencebandwidth); % get the interference power that contributes as interference to the signal % of interest from the total interference power interferencepowerrxinputactual = interferencepowerrxinputsumw*overlapfactor; % in w
interference is modeled by treating the contribution of interfering signal power in the overlapped bandwidth as noise. accordingly, add this quantity to the thermal noise at the ground station receiver input. note that the interference and noise power levels must be added in watts.
% calculate the thermal noise at the ground station receiver input. t = helpergetnoisetemperature(txmeofreq,rxgs); % in k kb = physconst("boltzmann"); thermalnoise = kb*t*txbandwidth; % in w % calculate the noise plus interference power at the receiver input. noiseplusinterferencepowerw = thermalnoise interferencepowerrxinputactual; % in w noiseplusinterferencepower = 10*log10(noiseplusinterferencepowerw); % in dbw
calculate the carrier to noise plus interference power spectral density ratio at the demodulator input as follows:
,
where:
is the carrier to noise plus interference power density ratio at the demodulator input (in db).
is the received downlink power from the meo satellite measured at the ground station receiver input (in dbw).
is the sum of receiver system thermal noise and the contribution of interfering signal power in the overlapped bandwidth measured at the receiver input (in dbw).
is the downlink transmission bandwidth from meo satellite (in hz).
is the loss that occurs between the receiver input and the demodulator input (in db).
% calculate loss that occurs between the receiver input and the demodulator % input. rxgsloss = rxgs.systemloss - rxgs.prereceiverloss; % calculate c/(n0 i0) at the demodulator input. cnoplusinterference = downlinkpowerrxinput - ... noiseplusinterferencepower 10*log10(txbandwidth) - rxgsloss;
calculate the energy per bit to noise plus interference power spectral density ratio at the demodulator input as follows:
,
where:
is the energy per bit to noise plus interference power spectral density ratio at the demodulator input (in db).
is the bit rate of the downlink from the meo satellite (in mbps).
bitrate = txmeosat.bitrate; ebnoplusinterference = cnoplusinterference - 10*log10(bitrate) - 60;
calculate the link margin as follows:
where:
is the link margin (in db).
is the minimum received energy per bit to noise power spectral density ratio at the demodulator input that is required to close the link (in db).
marginwithinterference = ebnoplusinterference - rxgs.requiredebno;
calculate the downlink closure status with interference. the status is true whenever the link margin is greater than or equal to 0 db.
downlinkstatuswithinterference = marginwithinterference >= 0;
calculate energy per bit to noise power spectral density ratio
calculate the energy per bit to noise power spectral density ratio (eb/n0) of the downlink and interfering links at the demodulator input for analysis later.
ebnodownlink = ebno(downlink); % in db ebnointerference = ebno(lnkinterference); % in db
plot downlink closure status with interference
plot the new downlink closure status that accounts for interference. compare the new link status with the previous case when interference was neglected.
plot(t,downlinkstatuswithinterference,"-r",t,downlinkstatus,"--g",linewidth=2); legend("interference accounted","interference neglected"); xlabel("time"); ylabel("downlink closure status"); title("link status as a function of time"); ylim([0 1.2]); grid on
when gaussian antenna or parabolic reflector are chosen
the plot shows that at 10:54 pm, the downlink cannot be closed because of excessive interference. this occurs because satellite 10
of the leo constellation flies overhead, and its transmission is picked up by its main lobe. this can also be visually confirmed by setting the current time of the viewer to 10:54 pm and clicking on the satellite near the main lobe of the antenna. note that you require antenna toolbox to select parabolic reflector.
if groundstationantennatype == "gaussian antenna" || groundstationantennatype == "parabolic reflector" v.currenttime = datetime(2021,3,17,22,54,0); end
gaussian antenna
parabolic reflector
when uniform rectangular array with interference cancellation is chosen
if you selected uniform rectangular array (requires phased array system toolbox), the plot shows that the downlink can be closed for the duration of the scenario because the array is pointing nulls towards the direction of the interfering leo satellites. this can also be visually confirmed by setting the current time of the viewer to 10:54 pm and 10:55 pm. to be able to manually set the viewer currenttime
, you must change autosimulate
to true. note that this will clear the simulation data. also, you will be required to re-compute the weights for these times and assign them to the array using pointat
. the satellite that is overflying the ground station is satellite 10
. click on it to see its name and orbit. drag the mouse while holding down on the left mouse button or scroll button to bring the camera to the desired position and orientation. rotate the scroll wheel to control camera zoom. additionally, make the radiation pattern opaque to clearly visualize the position of satellite 10
with respect to the lobes. you can see that at both times, satellite 10
is in between lobes. this is because the array is pointing a null towards the satellite, thereby cancelling interference from it.
if groundstationantennatype == "uniform rectangular array" % set autosimulate to true. sc.autosimulate = true; % set viewer currenttime to 10:54 pm. time = datetime(2021,3,17,22,54,0); v.currenttime = time; % calculate the weights and assign them to the array. currentinterferingsat = interferingsat(accessstatus(ac,time) == true); [azd,eld] = aer(rxgs,meosat,time,coordinateframe='body'); [azn,eln] = aer(rxgs,currentinterferingsat,time,coordinateframe='body'); % requires phased array system toolbox. wd = steervec(getelementposition(ant)/lambda,[wrapto180(azd);-eld]); wn = steervec(getelementposition(ant)/lambda,[wrapto180(azn)';-eln']); rn = (wn'*wn)\(wn'*wd); w = wd-wn*rn; pointat(rxgs,weights=w); % make the radiation pattern opaque. pattern(rxgs,txmeofreq, ... size = 1000000, ... transparency = 1); end
at 10:54 pm utc
you can run the above code with time set to 10:55 pm and observe the nulls pointing towards the new positions of the interfering satellites.
at 10:55 pm utc
calculate and plot carrier to noise ratio and carrier to noise plus interference ratio
calculate the carrier to noise ratio (cnr) at the demodulator input as follows:
,
,
where:
is the carrier to noise ratio at the demodulator input (in db).
is the carrier to noise power spectral density ratio at the demodulator input (in db).
% calculate the carrier to noise power spectral density ratio. cnodownlink = ebnodownlink 10*log10(bitrate) 60; % calculate the carrier to noise ratio. cbyn = cnodownlink - 10*log10(txbandwidth);
calculate the carrier to noise plus interference ratio (cnir) at the demodulator input as follows:
cbynplusi = cnoplusinterference - 10*log10(txbandwidth);
plot cnr and cnir
plot(t,cbynplusi,"-r",t,cbyn,"--g",linewidth=2); legend("cnir", "cnr",location="south"); xlabel("time"); ylabel("cnr or cnir (db)"); title("cnr and cnir vs. time for " groundstationantennatype); grid on
when uniform rectangular array (requires phased array system toolbox) with meo satellite tracking and interference cancellation is chosen, both cnr and cnir overlap because there is no interference from the leo satellites. this is because the array is pointing nulls towards the leo satellites. this can be confirmed by noting that the maximum power at ground station receiver input from the interfering satellites is about -167.3 dbw, which is very low. for all other antennas used in this example, the maximum power at ground station receiver input is much higher (-125.9 dbw for gaussian antenna, and -127.3 dbw for parabolic reflector).
maxinterferencepowerrxinput = max(interferencepowerrxinput,[],'all'); disp("the maximum power at ground station receiver input from the interfering leo satellites over the entire scenario duration is " maxinterferencepowerrxinput " dbw.");
the maximum power at ground station receiver input from the interfering leo satellites over the entire scenario duration is -125.9301 dbw.
compare link margins with and without interference
calculate the link margin without interference as follows:
marginwithoutinterference = ebnodownlink - rxgs.requiredebno;
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 vs. time for " groundstationantennatype); grid on
any time the link margin is greater than or equal to 0 db, the downlink is closed. with gaussian antenna, parabolic reflector (requires antenna toolbox), and uniform rectangular array (requires phased array system toolbox) without interference cancellation, there exist times when the link margin dips below 0 db because of interference. at these times, the downlink is broken.
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 specifications of the transmitter and receiver antennas
weights if using a uniform rectangular array
the signal and interference bandwidth
modify these parameters to observe their influence on the level of interference on the link. you can also choose the different antennas from antenna toolbox and phased array system toolbox for transmitters and receivers and observe the link performance. when using phased arrays and if you are only interested in making the main lobe track a single target and not deal with pointing nulls, you can use pointat
to automatically track other satellites, ground stations, and geographic locations without having to manually simulate by setting autosimulate
to false. the limitation of calling play
when using dynamically steered phased arrays is that you cannot visualize the variation of their radiation pattern over the course of the simulation.
helper functions
the example uses the helper function 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
see also
objects
- | | | | | | | |
functions
- | |