Введение

Чтобы правильно оценить эффективность радиолокационных и беспроводных систем связи, важно понять окружение распространения. Принятая степень сигнала моностатического радара задается уравнением радиолокационной области значений:

где - переданная степень, - коэффициент усиления антенны, - целевое радиолокационное сечение (RCS), - длина волны и - расстояние распространения. Все потери распространения, кроме потерь при распространении в свободном пространстве, включены в термин. Остальная часть примера показывает, как оценить этот термин в различных сценариях.

Потери при распространении в свободном пространстве

Потери при распространении в свободном пространстве вычисляются как функция расстояния распространения и частоты. В свободном пространстве радиочастотные сигналы распространяются со скоростью света во всех направлениях. На достаточно дальнем расстоянии источник излучения выглядит как точка в пространстве, и волновой фронт образует сферу, радиус которой равен. Плотность степени в волновом фронте обратно пропорциональна:

где - степень передаваемого сигнала. Для моностатического радара, где сигнал должен перемещаться в обоих направлениях (от источника к цели и назад), зависимость фактически обратно пропорциональна, как показано ранее в уравнении радара. Потеря, связанная с этим механизмом распространения, упоминается как потери при распространении в свободном пространстве, иногда также называемая потерями расширения. Количественно потери при распространении в свободном пространстве также являются функцией частоты, заданной [5]:

Как условие, потери на распространение часто выражаются в дБ. Эта конвенция значительно облегчает вывод двухсторонних потерь при распространении в свободном пространстве путем простого удвоения односторонних потерь свободного пространства.

Использование fspl вычислить потери при распространении в свободном пространстве и построить график потерь для частот от 1 до 1000 ГГц для различных областей значений.

c = physconst('lightspeed');
R0 = [100 1e3 10e3];
freq = (1:1000).'*1e9;
apathloss = fspl(R0,c./freq);
loglog(freq/1e9,apathloss);
grid on;
ylim([90 200]);
legend('Range: 100 m', 'Range: 1 km', 'Range: 10 km','Location','northwest');
xlabel('Frequency (GHz)');
ylabel('Path Loss (dB)');
title('Free Space Path Loss');

Рисунок иллюстрирует, что потери распространения увеличиваются с областью значений и частотой.

Потери распространения из-за осадков и атмосферы

В действительности сигналы не всегда перемещаются в вакууме, поэтому потери при распространении в свободном пространстве описывают только часть ослабления сигнала. Сигналы взаимодействуют с частицами в воздухе и теряют энергию вдоль пути распространения. Потеря изменяется с различными факторами, такими как давление, температура и плотность воды.

Потеря из-за дождя и снега

Дождь может быть основным ограничивающим фактором для радиолокационных систем, особенно при работе выше 5 ГГц. В модели ITU в [2] дождь характеризуется скоростью дождя (в мм/ч). Согласно [6], скорость дождя может варьироваться от менее 0,25 мм/ч для очень небольшого дождя до более 50 мм/ч для экстремальных дождей. В сложение, из-за формы падения дождя и его относительного размера по сравнению с длиной волны радиочастотного сигнала, потери распространения из-за дождя также являются функцией поляризации сигнала. В целом горизонтальная поляризация представляет худший случай потерь распространения из-за дождя.

Функции rainpl и cranerainpl может использоваться для вычисления потерь из-за дождя по моделям ITU и Crane, соответственно. Обе модели действительны между 1 ГГц и 1 ТГц. Пусть поляризация будет горизонтальной, поэтому угол наклона равен 0, и пусть сигнал распространяется параллельно земле, поэтому угол возвышения равен 0. Постройте графики потерь, вычисленных обеими моделями и сравните.

R0 = 5e3;                % 5 km range
rainrate = [1 4 20];     % rain rate in mm/h
el = 0;                  % 0 degree elevation
tau = 0;                 % horizontal polarization

for m = 1:numel(rainrate)
    rainloss_itu(:,m) = rainpl(R0,freq,rainrate(m),el,tau)';
    rainloss_crane(:,m) = cranerainpl(R0,freq,rainrate(m),el,tau)';
end
loglog(freq/1e9,rainloss_itu);
hold on;
set(gca,'ColorOrderIndex',1); % reset color index for better comparison
loglog(freq/1e9,rainloss_crane,'--');
hold off;
grid on;
legend('Light Rain (ITU)','Moderate Rain (ITU)','Heavy Rain (ITU)',...
    'Light Rain (Crane)','Moderate Rain (Crane)','Heavy Rain (Crane)', ...
    'Location','SouthEast');
xlabel('Frequency (GHz)');
ylabel('Attenuation at 1 km (dB)')
title('Rain Attenuation for Horizontal Polarization');

Потери, вычисленные с помощью модели Крана, в основном больше, чем потери, вычисленные с помощью модели ITU в этой области значений распространения. При меньших областях значений распространения и более низких частотах модель ITU может выдать меньшее значение ослабления, чем у Крана. Обратите внимание, что модели достаточно сильно различаются, что на более высоких частотах «легкое» количество осадков для одной модели может иметь такое же ослабление, как «умеренное» количество осадков для другой модели.

Подобно осадкам, снег также может оказать значительное влияние на распространение сигналов РФ. Распространенной практикой является обработка снега как осадков и вычисление потерь распространения на основе модели дождя, хотя этот подход имеет тенденцию немного переоценить потери. Ослабление из-за распространения через снег не рассматривается зависимым от поляризации, но сильно зависит от частоты. Модель потерь из-за снега параметризована эквивалентным содержимым жидкости вместо объема. Для заданного содержимого воды для снега требуется примерно в 10 раз больше объема, чем для дождя.

Использование snowpl вычислить потери из-за снега и построить график потерь по частоте. По умолчанию эта функция использует модель ослабления Gunn-East, которая обычно считается действительной примерно до 20 ГГц.

freq = (1:20)*1e9;
R0 = 1e3;               % 1 km range
snowrate = [0.1 1.5 4]; % equivalent liquid water content in mm/h

for m = 1:numel(snowrate)
    snowloss(:,m) = snowpl(R0,freq,snowrate(m));
end
loglog(freq/1e9,snowloss);
grid on;
legend('Light Snow','Moderate Snow','Heavy Snow', ...
    'Location','SouthEast');
xlabel('Frequency (GHz)');
ylabel('Attenuation at 1 km (dB)')
title('Snow Attenuation');

Потеря из-за тумана и облака

Туман и облако образуются и с каплями воды, хотя и гораздо меньше по сравнению с каплями дождя. Размер капель тумана обычно составляет менее 0,01 см. Туман часто характеризуется плотностью воды. Средний туман с видимостью около 300 метров, имеет плотность жидкой воды 0,05 г/м ^ 3. Для сильного тумана, где видимость падает до 50 метров, плотность жидкой воды составляет около 0,5 г/м ^ 3. Температура атмосферы (в Цельсии) также присутствует в модели ITU для потерь распространения из-за тумана и облака [3].

Использование fogpl вычислить потери из-за тумана и построить график потерь по частоте. Модель ITU для ослабления из-за тумана действительна между 10 ГГц и 1 ТГц.

freq = (10:1000)*1e9;
T = 15;                         % 15 degree Celsius
waterdensity = [0.01 0.05 0.5]; % liquid water density in g/m^3
for m = 1: numel(waterdensity)
    fogloss(:,m) = fogpl(R0,freq,T,waterdensity(m))';
end
loglog(freq/1e9,fogloss);
grid on;
legend('Light Fog','Medium Fog','Heavy Fog','Location','southeast');
xlabel('Frequency (GHz)');
ylabel('Attenuation at 1 km (dB)')
title('Fog Attenuation');

Обратите внимание, что в целом тумана нет, когда идет дождь.

Потеря из-за атмосферной абсорбции и линзирования

Даже когда нет тумана или дождя, атмосфера полна газов, которые все еще влияют на распространение сигнала. Модель ITU [4] описывает ослабление атмосферного газа как функцию от давления сухого воздуха, такого как кислород, измеренный в hPa, и плотности водяного пара, измеренной в г/м ^ 3.

Использование tropopl вычислить потери из-за атмосферного поглощения и построить график потерь по частоте. По умолчанию эта функция использует модель Mean Annual Global Reference Atmosphere (MAGRA), чтобы получить типичные значения температуры, давления и плотности водяного пара для заданной высоты. Мы также можем задать модель широты, чтобы использовать модель, адаптированную для определенной области значений широт. Некоторые модели широты также допускают спецификацию сезона. Пусть наша высота будет 2 км (обратите внимание, что тропосфера, для которой эта модель действительна, простирается до 10 км) и наш путь распространения будет подавлен на 5 степени. Эта функция вернёт общие потери из-за атмосферного поглощения по наклонному пути распространения, но не включает рассеивание из-за преломления (линзирования). Мы сравним потери между моделями низкой, средней и высокой широты.

height = 2e3;
el = -5; % elevation angle
atmloss_low = tropopl(R0,freq,height,el,'LatitudeModel','Low');
atmloss_mid = tropopl(R0,freq,height,el,'LatitudeModel','Mid');
atmloss_high = tropopl(R0,freq,height,el,'LatitudeModel','High');
loglog(freq/1e9,atmloss_low);
hold on;
loglog(freq/1e9,atmloss_mid);
loglog(freq/1e9,atmloss_high);
hold off;
grid on;
legend('Low Latitudes','Mid Latitudes','High Latitudes','Location','northwest');
xlabel('Frequency (GHz)');
ylabel('Attenuation at 1 km (dB)')
title('Atmospheric Gas Attenuation');

График предполагает, что происходит сильное поглощение из-за атмосферных газов на частоте около 60 ГГц.

Другой источник потерь из-за атмосферы - от атмосферного линзирования. Это является явлением, посредством которого угловая протяженность передачи увеличивается с областью значений из-за градиента преломления. Это расширение энергии уменьшает плотность энергии вдоль номинального (прямого) пути распространения, независимо от частоты.

Атмосферное давление и, таким образом, рефрактивность изменяются с высотой. Так что для заданной высоты угла возвышения пути распространения достаточно, чтобы определить потери из-за этого эффекта.

Использование lenspl вычислить эти потери и построить график по частоте. Поскольку эта потеря не зависит от частоты, мы построим график потерь от области значений распространения для набора высот. Используйте угол возвышения 0,05 степеней для наклонного пути распространения.

R = 1e3:1e3:100e3;      % propagation range
el = 0.05;              % elevation angle
heights = [10 100 200]; % radar platform heights
for m = 1:numel(heights)
    lenloss(:,m) = lenspl(R,heights(m),el);
end
semilogy(R/1e3,lenloss);
grid on;
legend('Height: 10 m','Height: 100 m','Height: 200 m','Location',...
    'southeast');
xlabel('Propagation Range (km)');
ylabel('Attenuation (dB)')
title('Atmospheric Lensing Attenuation');

Ослабление из-за линзирования уменьшается с увеличением высоты. Для удобства ослабление из-за линзирования также обеспечивается как вторичный выход от tropopl.

Потеря из-за несоответствия поляризации

Некоторые типы потерь распространения зависят от поляризации передаваемого излучения, например, с потерями дождя. Это является результатом химических и структурных свойств среды. Однако даже в свободном пространстве могут быть потери из-за несоответствия распространенного вектора поляризации и поляризации приемной антенны. Для примера, если распространенный вектор поляризации ортогональен поляризации приемной антенны, энергия прямого сигнала не будет принята. Обратите внимание, что «распространенный вектор поляризации» в целом не такой, как переданный вектор поляризации, так как направление распространения должно быть учтено. Обратите внимание, что другие функции потерь, которые берут поляризацию как вход, не вычисляют потери из-за этого несоответствия. Поляризационно-зависимые потери из-за свойств среды распространения могут обрабатываться отдельно от потерь из-за несоответствия поляризации, поскольку последняя сильно зависит от ориентации передатчика/приемника.

Использование polloss вычислить потери из-за несоответствия поляризации для заданной поляризации передачи/приема, положения платформы и ориентации платформы. Поместите передающую платформу в источник без вращения от инерционной. Поместите приемную платформу вдоль оси X и вычислите поляризационные потери для области значений углов крена. Пусть обе поляризации антенны будут вертикальными.

poltx = [0;1];  % [H;V] polarization
polrx = [0;1];
postx = [0;0;0];
posrx = [100;0;0];
frmtx = eye(3); % transmit frame aligned with inertial
rolls = 0:180;

for m = 1:numel(rolls)
    frm_r = rotx(rolls(m));
    rho(m) = polloss(poltx,polrx,posrx,frm_r,postx,frmtx);
end

semilogy(rolls,rho);
grid on;
xlabel('Roll Angle (deg)');
ylabel('Attenuation (dB)');
title('Attenuation Due to Polarization Mismatch');

Затухание приближается к бесконечности под углом крена 90 степеней.

Коэффициент распространения радара и схема вертикального покрытия

При передаче на широкий угол или от антенны, близкой к земле, многолучевой канал от отскока земли, наряду с преломлением от атмосферы, приводит к диаграмме направленности излучения в заданной области значений, которая может сильно отличаться от номинального шаблона. Это фиксируется коэффициентом распространения радара, который является отношением фактической напряженности поля относительно того, какая напряженность поля была бы в свободном пространстве. Коэффициент распространения может сильно изменяться, когда изменяется относительная фаза между сигналами прямого и косвенного пути.

Вертикальная схема покрытия (график Блейка) является компактным способом отображения контуров энергии фиксированного сигнала (такая как минимальная степень сигнала для обнаружения) как функции области значений и угла возвышения распространения. Необходимо учитывать только вертикальную плоскость, в которой распространяются как сигналы прямого, так и косвенного пути.

Функция radarvcd принимает ссылку области значений как вход и возвращает область значений, в которой мощности приемника в многолучевой среде равны тому, что было бы в свободном пространстве. Этот «эффективный» диапазон нанесен на график Угол высоты области значений. Это может быстро дать, например, фактический диапазон обнаружения, заданный диапазон обнаружения свободного пространства, как функцию области значений, высоты или угла возвышения.

Мы будем использовать область значений обнаружения свободного пространства 100 км, частоты передачи в L-диапазоне и C-диапазоне и высоту антенны 12 м. По умолчанию используется шаблон передачи sinc.

freq = [1.06 5.7]*1e9; % L-Band and C-Band transmit frequencies (Hz)
antht = 12;            % height of antenna (m)
rngfs = 100;           % free-space detection range (km)
for m = 1:numel(freq)
    [vcp{m}, vcpang{m}] = radarvcd(freq(m),rngfs,antht);
end

blakechart принимает эти области значений и углы обнаружения вместе с дополнительными атмосферными свойствами для создания графика Блейка. Использование refractiveidx и вычислите соответствующий показатель преломления для входов blakechart.

[~,N] = refractiveidx(0); % atmospheric refractivity at the surface
helperPlotBlakeChart(vcp,vcpang,N)

Помеха при отказе от земли доминирует над коэффициентом распространения для более коротких областей значений, «интерференционной области», но при более длинных областях значений и низких углах возвышения коэффициент распространения становится доминирующим за счет дифракции над горизонтом, «дифракционной областью». Мы можем использовать radarpropfactor вычислить коэффициент распространения для интервала областей значений и наблюдать различие между этими двумя областями.

Вычислите коэффициент распространения для фиксированной высоты над поверхностью 1 км и областей значений распространения между 50 и 200 км. Установите значение стандартного отклонения наклона и высоты поверхности равным 0, чтобы представлять собой гладкую поверхность. Выполните анализ для двух полос частот.

tgtht = 1e3;           % target height (m)
R = (50:200)*1e3;      % propagation range (m)
Re = effearthradius;   % effective Earth radius (m)
Rd = sqrt(2*Re)*(sqrt(antht) + sqrt(tgtht)); % diffraction range
F = zeros(numel(freq),numel(R));
for m = 1:numel(freq)
    F(m,:) = radarpropfactor(R,freq(m),antht,tgtht,'SurfaceHeightStandardDeviation',0,'SurfaceSlope',0);
end
helperPlotPropagationFactor(R,F,Rd)

Коэффициент распространения колеблется в интерференционной области, затем быстро уменьшается в дифракционной области.

Объедините помехи от отскока земли и потери атмосферного поглощения. Предположим в этом вычислении, что радар надводного корабля диапазона S 3,3 ГГц на 20 м выше воды и имеет повышение ширину луча 30 o.

freq = 3.3e9;                   % Frequency (Hz)
elbw = 30;                      % Elevation beamwidth (deg)
Rkm = 1:0.1:120;                % Range (km)
R = Rkm.*1e3;                   % Range (m)
[htsd,beta0] = searoughness(1); % Sea surface
anht = 20 + 2*htsd;             % Radar height (m)
tgtht = (anht+1):1:300;         % Target height (m)

% Calculate combined environment losses for different heights and ranges
[PLdB, PLdBNorm] = helperCombineEnvLosses(R,freq,anht,tgtht,htsd,beta0,elbw);

% Plot combined losses for different heights and ranges
helperPlotCombinedEnvLosses(Rkm,freq,anht,tgtht,PLdBNorm)

Многолучевое распространение, задержки и допплеровские сдвиги

Как мы видели, сигналы могут не всегда распространяться вдоль линии зрения, но поступать в пункт назначения через различные пути и могут складываться или конструктивно, или разрушительно. Этот многолучевой эффект может вызвать значительные колебания в степень принимаемого сигнала.

Функции, упомянутые выше для вычисления потерь распространения, полезны для установления бюджетных ссылок, но чтобы симулировать распространение произвольных сигналов, нам также нужно применить зависящие от области значений задержки времени, усиления и сдвиги фазы. Различные объекты канала доступны для моделирования распространения многолучевого распространения. Для простого пути видимости, phased.LOSChannel может использоваться, чтобы смоделировать распространение, подверженное любому из типов потерь, описанных выше.

Отражение Земли является общим явлением для многих радиолокационных или беспроводных систем связи. Например, когда радар наземного или морского базирования освещает цель, сигнал не только распространяется вдоль прямой линии зрения, но также отражается от земли. twoRayChannel может использоваться, чтобы смоделировать комбинацию прямого пути и пути с одним отскоком, например, для отражения земли.

Задержки и допплеровские сдвиги

Сначала задайте переданный сигнал. В этом случае будет использоваться прямоугольная форма волны.

waveform = phased.RectangularWaveform('PRF',250);
wav = waveform();

Предположим, что рабочая частота L-диапазона составляет 1,9 ГГц. Моделируйте канал с

fc = 1.9e9;
channel = twoRayChannel('PropagationSpeed',c,'OperatingFrequency',fc);

Предположим, что целевой модуль находится на высоте 1,65 км над землей, антенна радара на расстоянии 50 км находится на высоте 12 метров над землей. Симулируйте сигнал, когда он достигает цели.

pos_radar = [0;0;12];
pos_target = [50e3;0;1.65e3];
vel_radar = [0;0;0];
vel_target = [-200;0;0];
y2ray = channel(wav,pos_radar,pos_target,vel_radar,vel_target);

Визуализируйте переданный и распространенный импульс и их нормированный спектр. Канал ввел задержку 167, которая соответствует 50-километровой области значений цели, разделенной на скорость света.

[delay, dop] = helperPlotDelayAndDopplerShift(wav,y2ray,waveform.SampleRate);

estRange = delay*c*1e-3 % km

estRange =

   49.9954

Канал также применил доплеровский сдвиг, который соответствует скорости области значений цели. Сравните расчетное значение с основной истиной -200 м/с с помощью dop2speed и freq2wavelen.

estRangeRate = -dop2speed(dop,freq2wavelen(fc)) % m/s

estRangeRate =

 -201.9038

Многолучевое замирание

Потеря сигнала, понесенная в этом канале, может быть вычислена как

L_2ray = pow2db(bandpower(wav))-pow2db(bandpower(y2ray))

L_2ray =

  151.5888

Это потери при распространении в свободном пространстве задаются как

L_ref = fspl(norm(pos_target-pos_radar),c/fc)

L_ref =

  132.0069

Результат предполагает, что в этом строении канал вводит дополнительные потери на 19,6 дБ к полученному сигналу по сравнению со случаем свободного пространства. Теперь предположим, что цель летит немного выше на 1,8 км над землей. Повторение симуляции выше предполагает, что на этот раз отражение земли фактически обеспечивает коэффициент усиления 6 дБ! Хотя потери при распространении в свободном пространстве по существу одинаковы в двух сценариях, перемещение на 150 м вызвало колебание степени сигнала на 25,6 дБ.

pos_target = [50e3;0;1.8e3];
y2ray  = channel(wav,pos_radar,pos_target,vel_radar,vel_target);
L_2ray = pow2db(bandpower(wav))-pow2db(bandpower(y2ray))
L_ref  = fspl(norm(pos_target-pos_radar),c/fc)

L_2ray =

  126.0374


L_ref =

  132.0078

Увеличение полосы пропускания системы увеличивает пропускную способность ее канала. Это позволяет повысить скорость передачи данных в системах связи и более мелкие разрешения области значений для радиолокационных систем. Увеличенная полоса пропускания может также улучшить робастность к многолучевому замиранию для обеих систем.

Как правило, широкополосные системы работают с шириной полосы более 5% от их центральной частоты. Напротив, узкополосные системы работают с пропускной способностью 1% или менее от центральной частоты системы.

Было показано, что узкополосный канал в предыдущей секции очень чувствителен к многолучевому замиранию. Незначительные изменения высоты цели привели к значительным потерям сигнала. Характеристики замирания канала могут быть нанесены путем изменения высоты цели в диапазоне рабочих высот для этой радиолокационной системы. Выбран размах высот от 1 км до 3 км.

% Simulate the signal fading at the target for heights from 1 km to 3 km
hTarget = linspace(1e3,3e3);
pos_target = repmat([50e3;0;1.6e3],[1 numel(hTarget)]);
pos_target(3,:) = hTarget;
vel_target = repmat(vel_target,[1 numel(hTarget)]);

release(channel);
y2ray = channel(repmat(wav,[1 numel(hTarget)]),pos_radar,pos_target,vel_radar,vel_target);

Потеря сигнала, наблюдаемая в цели для узкополосной системы, теперь может быть нанесена.

L2ray = pow2db(bandpower(wav))-pow2db(bandpower(y2ray));

clf;
plot(hTarget,L2ray);
xlabel('Target Height (m)');
ylabel('One-Way Propagation Loss (dB)');
title('Multipath Fading Observed at the Target');
grid on;

Чувствительность потерь канала к высоте цели для этой узкополосной системы четко. Глубокие затухания сигнала происходят на высотах, которые, вероятно, находятся в пределах зоны наблюдения радара.

Увеличение пропускной способности канала может улучшить робастность к этим многолучевым исчезновениям. Для этого задается широкополосная форма волны с пропускной способностью 8% от центральной частоты ссылки.

bw = 0.08*fc;
pulse_width = 1/bw;
fs = 2*bw;

waveform = phased.RectangularWaveform('SampleRate',fs,'PRF',2000,'PulseWidth',pulse_width);
wav = waveform();

Широкополосная версия этой модели канала, widebandTwoRayChannelможет использоваться, чтобы симулировать многолучевые отражения этого широкополосного сигнала вне земли между радаром и целью и вычислить соответствующие потери канала.

channel = widebandTwoRayChannel('PropagationSpeed',c,'OperatingFrequency',fc,'SampleRate',fs);

Сигнал в цели для различных рабочих высот теперь может быть моделирован для этой широкополосной системы.

y2ray_wb = channel(repmat(wav,[1 numel(hTarget)]),pos_radar,pos_target,vel_radar,vel_target);
L2ray_wb = pow2db(bandpower(wav))-pow2db(bandpower(y2ray_wb));

hold on;
plot(hTarget,L2ray_wb);
hold off;
legend('Narrowband','Wideband');

Как ожидалось, широкополосный канал обеспечивает гораздо лучшую эффективность на широкой области значений высот для цели. На самом деле, когда высота цели увеличивается, влияние многолучевого замирания почти полностью исчезает. Это связано с тем, что различие в задержке распространения между сигналами прямого и обратного путей увеличивается, уменьшая величину когерентности между этими двумя сигналами при приеме в цель.

Заключение

В этом примере представлен обзор потерь распространения РФ из-за атмосферных и погодных эффектов. Это также вводит многолучевые колебания сигнала из-за отскока земли. Он осветил функции и объекты, чтобы симулировать потери ослабления для узкополосных и широкополосных каналов с одним отказом.

Ссылки

[1] John Seybold, Introduction to RF Propagation, Wiley, 2005

[2] Recommendation ITU-R P.838-3, 2005

[3] Recommendation ITU-R P.840-3, 2013

[4] Recommendation ITU-R P.676-10, 2013

[5] Recommendation ITU-R P.525-2, 1994

[6] Rain, A Water Resource (Pamphlet), U.S. Geological Survey, 1988

Вспомогательные функции

helperPlotPropagationFactor

function helperPlotPropagationFactor(R,F,Rd)

% Plot interference and diffraction region patches
[minF, maxF] = bounds(F(:));
maxF = ceil((maxF+10)/10)*10;
minF = floor((minF-10)/10)*10;
yPatch = [minF minF maxF maxF];
c1 = [0.3010 0.7450 0.9330];
c2 = [0 0.4470 0.7410];
clf % clear current figure
fill([R(1) Rd Rd R(1)]/1e3,yPatch,c1,'EdgeColor','none','FaceAlpha',0.25)
hold on
fill([Rd R(end) R(end) Rd]/1e3,yPatch,c2,'EdgeColor','none','FaceAlpha',0.25)

% Plot one-way propagation factor
set(gca,'ColorOrderIndex',1); % reset color index
plot(R/1e3,F);
ylim([minF maxF])
grid on;
xlabel('Range (km)');
ylabel('Propagation Factor (dB)');
title('One-Way Propagation Factor at 1 km above Surface');
legend('Interference Region', 'Diffraction Region',...
    'L-Band (1.06 GHz)', 'C-Band (5.7 GHZ)',...
    'Location','SouthWest')
hold off
end

helperPlotBlakeChart

function helperPlotBlakeChart(vcp,vcpang,N)
% Calculate refraction exponent
DelN = -7.32*exp(0.005577*N);
rexp = log(N./(N + DelN));

subplot(211)
blakechart(vcp{1},vcpang{1},'SurfaceRefractivity',N,'RefractionExponent',rexp);
legend('L-Band (1.06 GHz)')
xlabel('')
title ('Blake Chart - Antenna Height: 12 m')
subplot(212)
blakechart(vcp{2},vcpang{2},'SurfaceRefractivity',N,'RefractionExponent',rexp);
allc = get(gca,'Children');
set(allc(11),'Color',[0.8500 0.3250 0.0980]) % Change line color
title('')
legend('C-Band (5.7 GHz)')
end

helperPlotDelayAndDopplerShift

function [delay, dop] = helperPlotDelayAndDopplerShift(wav,y2ray,Fs)
% Plot transmitted and propagated pulse
t = 1e6*(0:numel(wav)-1)'/Fs;
subplot(211)
yyaxis left
plot(t,abs(wav))
ylabel('Magnitude')
yyaxis right
plot(t,abs(y2ray))
grid on
axis padded
xlim([0 300])
xlabel(['Time (' char(0x00B5) 's)'])
ylabel('Magnitude')
title('Transmitted and Propagated Pulse')

% Annotation
delay = midcross(abs(y2ray),t/1e6,'MidPercentReferenceLevel',80); % seconds
delay = delay(1);
xl = xline(1e6*delay,'-.',... % Annotation
    {[num2str(round(1e6*delay)),' ',char(0x00B5) 's delay']},'Color',[0.8500 0.3250 0.0980]);
xl.LabelVerticalAlignment = 'middle';
xl.LabelHorizontalAlignment = 'left';
xl.LineWidth = 2;

% Plot power spectrum
subplot(212)
[p,f] = pspectrum([wav y2ray],Fs,'FrequencyLimits',[-20e3 20e3]);
p = abs(p);
plot(1e-3*f,rescale(p,'InputMin',min(p),'InputMax',max(p)));
axis padded
grid on
[~,idx]=max(p);
dop = f(idx(2))-f(idx(1));   % Hz
xlabel('Frequency (kHz)')
ylabel('Magnitude')
title('Normalized Spectrum')

xl = xline(1e-3*dop,'-.',... % Annotation
    {'Doppler shift',[num2str(round(dop)*1e-3) ' kHz']},'Color',[0.8500 0.3250 0.0980]);
xl.LabelVerticalAlignment = 'bottom';
xl.LineWidth = 2;
legend('Transmitted','Propagated')
end

helperCombineEnvLosses

function [PLdB, PLdBNorm] = helperCombineEnvLosses(R,freq,anht,tgtht,htsd,beta0,elbw)
% Calculate the combined environment losses
numHt = numel(tgtht);
numR = numel(R);
F = zeros(numHt,numR);
for ih = 1:numHt
    F(ih,:) = radarpropfactor(R, freq, anht, tgtht(ih),...
    'SurfaceHeightStandardDeviation',htsd,'SurfaceSlope',beta0,...
    'ElevationBeamwidth', elbw);
end

% Free space spreading loss
Lspl_dB = 2*fspl(R,freq2wavelen(freq));  % Factor of 2 for two-way

% Perform tropospheric losses calculation for a subset of elevation angles,
% since the ray refracting can take a long time.
numEl = 10;
minEl = height2el(tgtht(1),anht,R(end)); % Min elevation angle (deg)
maxEl = height2el(tgtht(end),anht,R(1)); % Max elevation angle (deg)
elSubset = linspace(minEl,maxEl,numEl);
LtropoSubset = zeros(numEl,numR);
for ie = 1:numEl
    LtropoSubset(ie,:) = tropopl(R,freq,anht,elSubset(ie));
end
% Interpolate tropospheric losses for all elevation angles of interest
Ltropo = zeros(numHt,numR);
for ir = 1:numR
    el = height2el(tgtht,anht,R(ir));
    Ltropo(:,ir) = interp1(elSubset,LtropoSubset(:,ir),el);
end

PLdB = 2*F - Lspl_dB - Ltropo;           % Factor of 2 for two-way
PLdBNorm = PLdB - max(PLdB(:));
end

helperPlotCombinedEnvLosses

function helperPlotCombinedEnvLosses(Rkm,freq,anht,tgtht,PLdBNorm)
% Plot combined losses for different heights and ranges
hP = pcolor(Rkm,tgtht,PLdBNorm);
set(hP, 'EdgeColor', 'none');
title([num2str(freq/1e9) ' GHz S-Band Radar'])
subtitle([num2str(round(anht)) ' m above water'])
xlabel('Range (km)')
ylabel('Height (m)')
colormap('jet');
caxis([-150 0])
hC = colorbar;
hC.Label.String = 'Normalized Two-Way Propagation Loss (dB)';
end