Задайте требования к радару

Установите желаемую вероятность обнаружения равной 0,9 и вероятность ложного предупреждения равной $$10^{-6}$$. Установите максимальную область значений радара равную 4000 метрам, а разрешение области значений - 50 метрам. Установите фактическую целевую область значений 3000 метров. Установите поперечное сечение радара цели на 1,5 квадратных метра и установите рабочую частоту на 10 ГГц. Все расчеты выполняются в основной полосе частот.

c = physconst('LightSpeed');
pd = 0.9;
pfa = 1e-6;
max_range = 4000;
target_range = 3000.0;
range_res = 50;
tgt_rcs = 1.5;
fc = 10e9;
lambda = c/fc;

Любая симуляция, которая вычисляет PFA и pd, требует обработки многих сигналов. Чтобы снизить требования к памяти, обрабатывайте сигналы в фрагментах импульсов. Установите количество импульсов для обработки равным 45000 и установите размер каждого фрагмента равным 10000.

Npulse = 45000;
Npulsebuffsize = 10000;

Выберите параметры формы волны и сигнала

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

pulse_bw = c/(2*range_res);
prf = c/(2*max_range);
fs = 2*pulse_bw;
pulse_duration = 10/pulse_bw;
waveform = phased.LinearFMWaveform('PulseWidth',pulse_duration,...
    'SampleRate',fs,'SweepBandwidth',...
    pulse_bw,'PRF',prf);

Достижение конкретных Pd и PFA требует, чтобы достаточная степень сигнала поступала в приемник после того, как цель отражает сигнал. Вычислите минимальный ОСШ, необходимый для достижения заданной вероятности ложного предупреждения и вероятности обнаружения с помощью уравнения Альберсхайма.

snr_min = albersheim(pd,pfa);

Чтобы достичь этого ОСШ, достаточная степень должна быть передана к цели. Используйте основное уравнение радиолокации, чтобы оценить пиковую степень передачи, peak_power, требуемый для достижения заданного ОСШ в дБ для цели в области значений 3000 метров. Принятый сигнал также зависит от поперечного сечения радара цели (RCS). который принимается в соответствии с неколеблющейся моделью (Swerling 0). Установите радар на одинаковые коэффициенты усиления передачи и приема 20 дБ. Дано основное уравнение радиолокации

txrx_gain = 20;
peak_power = ((4*pi)^3*noisepow(1/pulse_duration)*target_range^4*...
    db2pow(snr_min))/(db2pow(2*txrx_gain)*tgt_rcs*lambda^2)

peak_power = 293.1830

Настройте системные объекты передатчика

Создайте Системные Объекты, которые составляют передаточную часть симуляции: радарную платформу, антенну, передатчик и излучателя.

antennaplatform = phased.Platform(...
    'InitialPosition',[0; 0; 0],...
    'Velocity',[0; 0; 0]);
antenna = phased.IsotropicAntennaElement(...
    'FrequencyRange',[5e9 15e9]);
transmitter = phased.Transmitter(...
    'Gain',txrx_gain,...
    'PeakPower',peak_power,...
    'InUseOutputPort',true);
radiator = phased.Radiator(...
    'Sensor',antenna,...
    'OperatingFrequency',fc);

Настройка целевого системного объекта

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

target{1} = phased.RadarTarget(...
    'MeanRCS',tgt_rcs,...
    'OperatingFrequency',fc);
targetplatform{1} = phased.Platform(...
    'InitialPosition',[target_range; 0; 0]);
target{2} = phased.RadarTarget(...
    'MeanRCS',0,...
    'OperatingFrequency',fc);
targetplatform{2} = phased.Platform(...
    'InitialPosition',[target_range; 0; 0]);

Настройка системных объектов распространения свободного пространства

Моделируйте окружение распространения от радара до целей и назад.

channel{1} = phased.FreeSpace(...
    'SampleRate',fs,...
    'TwoWayPropagation',true,...
    'OperatingFrequency',fc);
channel{2} = phased.FreeSpace(...
    'SampleRate',fs,...
    'TwoWayPropagation',true,...
    'OperatingFrequency',fc);

Настройка системных объектов приемника

Задал шум путем установки NoiseMethod свойство к 'Noise temperature' и ReferenceTemperature свойство до 290 K.

collector = phased.Collector(...
    'Sensor',antenna,...
    'OperatingFrequency',fc);
receiver = phased.ReceiverPreamp(...
    'Gain',txrx_gain,...
    'NoiseMethod','Noise temperature',...
    'ReferenceTemperature',290.0,...
    'NoiseFigure',0,...
    'SampleRate',fs,...
    'EnableInputPort',true);
receiver.SeedSource = 'Property';
receiver.Seed = 2010;

Задайте сетку быстрого времени

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

fast_time_grid = unigrid(0,1/fs,1/prf,'[)');
rangebins = c*fast_time_grid/2;

Создайте переданный импульс из формы волны

Создайте сигнал, который вы хотите передать.

wavfrm = waveform();

Создайте переданный сигнал, который включает в себя усиление переданной антенны.

[sigtrans,tx_status] = transmitter(wavfrm);

Создайте согласованный фильтр коэффициенты из Системного объекта формы волны. Затем создайте согласованный фильтр System object™.

MFCoeff = getMatchedFilter(waveform);
matchingdelay = size(MFCoeff,1) - 1;
filter = phased.MatchedFilter(...
    'Coefficients',MFCoeff,...
    'GainOutputPort',false);

Вычисление целевого интервала области значений

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

ant_pos = antennaplatform.InitialPosition;
ant_vel = antennaplatform.Velocity;
tgt_pos = targetplatform{1}.InitialPosition;
tgt_vel = targetplatform{1}.Velocity;
[tgt_rng,tgt_ang] = rangeangle(tgt_pos,ant_pos);
rangeidx = val2ind(tgt_rng,rangebins(2)-rangebins(1),rangebins(1));

Цикл по импульсам

Создайте цикл обработки сигналов. Каждый шаг выполняется путем выполнения системных объектов. Цикл обрабатывает импульсы дважды, один раз для условия присутствия цели и один раз для условия отсутствия цели.

rcv_pulses = zeros(length(sigtrans),Npulsebuffsize);
h1 = zeros(Npulse,1);
h0 = zeros(Npulse,1);
Nbuff = floor(Npulse/Npulsebuffsize);
Nrem = Npulse - Nbuff*Npulsebuffsize;
for n = 1:2 % H1 and H0 Hypothesis
    trsig = radiator(sigtrans,tgt_ang);
    trsig = channel{n}(trsig,...
        ant_pos,tgt_pos,...
        ant_vel,tgt_vel);
    rcvsig = target{n}(trsig);
    rcvsig = collector(rcvsig,tgt_ang);
    
    for k = 1:Nbuff
        for m = 1:Npulsebuffsize
            rcv_pulses(:,m) = receiver(rcvsig,~(tx_status>0));
        end
        rcv_pulses = filter(rcv_pulses);
        rcv_pulses = buffer(rcv_pulses(matchingdelay+1:end),size(rcv_pulses,1));
        if n == 1
            h1((1:Npulsebuffsize) + (k-1)*Npulsebuffsize) = rcv_pulses(rangeidx,:).';
        else
            h0((1:Npulsebuffsize) + (k-1)*Npulsebuffsize) = rcv_pulses(rangeidx,:).';
        end
    end
    if (Nrem > 0)
        for m = 1:Nrem
            rcv_pulses(:,m) = receiver(rcvsig,~(tx_status>0));
        end
        rcv_pulses = filter(rcv_pulses);
        rcv_pulses = buffer(rcv_pulses(matchingdelay+1:end),size(rcv_pulses,1));
         if n == 1
            h1((1:Nrem) + Nbuff*Npulsebuffsize) = rcv_pulses(rangeidx,1:Nrem).';
        else
            h0((1:Nrem) + Nbuff*Npulsebuffsize) = rcv_pulses(rangeidx,1:Nrem).';
         end
    end
end

Создайте гистограмму Согласованного фильтра Выходов

Вычислите гистограммы возвратов target-present и target-dissent. Используйте 100 интервалов, чтобы дать приблизительную оценку расширения значений сигналов. Установите область значений значений гистограммы от наименьшего сигнала до наибольшего сигнала.

h1a = abs(h1);
h0a = abs(h0);
thresh_low = min([h1a;h0a]);
thresh_hi  = max([h1a;h0a]);
nbins = 100;
binedges = linspace(thresh_low,thresh_hi,nbins);
figure
histogram(h0a,binedges)
hold on
histogram(h1a,binedges)
hold off
title('Target-Absent Vs Target-Present Histograms')
legend('Target Absent','Target Present')

Сравнение моделируемых и теоретических Pd и PFA

Чтобы вычислить Pd и PFA, вычислите количество образцов, когда возврат без целевого значения и возврат после целевого значения превышают заданный порог. Этот набор порогов имеет более мелкую гранулярность, чем интервалы, используемые для создания гистограммы в предыдущей симуляции. Затем нормализуйте эти счетчики по количеству импульсов, чтобы получить оценку вероятностей. Векторная sim_pfa - моделируемая вероятность ложного предупреждения как функция от порога, thresh. Векторная sim_pd является моделируемой вероятностью обнаружения, также функцией порога. Приемник устанавливает порог так, чтобы он мог определить, присутствует или отсутствует цель. Гистограмма выше предполагает, что лучший порог около 1,8.

nbins = 1000;
thresh_steps = linspace(thresh_low,thresh_hi,nbins);
sim_pd = zeros(1,nbins);
sim_pfa = zeros(1,nbins);
for k = 1:nbins
    thresh = thresh_steps(k);
    sim_pd(k) = sum(h1a >= thresh);
    sim_pfa(k) = sum(h0a >= thresh);
end
sim_pd = sim_pd/Npulse;
sim_pfa = sim_pfa/Npulse;

Чтобы построить график экспериментальной кривой ROC, необходимо инвертировать кривую PFA, чтобы можно было построить график Pd против PFA. Инвертировать кривую PFA можно только, когда можно выразить PFA как строго монотонную функцию уменьшения thresh. Чтобы выразить PFA таким образом, найдите все индексы массива, где PFA является константой по соседним индексам. Затем удалите эти значения из массивов Pd и PFA.

pfa_diff = diff(sim_pfa);
idx = (pfa_diff == 0);
sim_pfa(idx) = [];
sim_pd(idx) = [];

Ограничьте наименьший PFA, $$10^{-6}$$.

minpfa = 1e-6;
N = sum(sim_pfa >= minpfa);
sim_pfa = fliplr(sim_pfa(1:N)).';
sim_pd = fliplr(sim_pd(1:N)).';

Вычислите теоретические значения PFA и Pd от наименьших PFA до 1. Затем постройте график теоретической кривой PFA.

[theor_pd,theor_pfa] = rocsnr(snr_min,'SignalType',...
    'NonfluctuatingNoncoherent',...
    'MinPfa',minpfa,'NumPoints',N,'NumPulses',1);
semilogx(theor_pfa,theor_pd)
hold on
semilogx(sim_pfa,sim_pd,'r.')
title('Simulated and Theoretical ROC Curves')
xlabel('Pfa')
ylabel('Pd')
grid on
legend('Theoretical','Simulated','Location','SE')

Улучшите симуляцию с использованием миллиона импульсов

В предыдущей симуляции значения Pd при низком PFA не падают вдоль гладкой кривой и даже не растягиваются до заданного режима работы. Причиной этого является то, что при очень низких уровнях PFA очень немногие, если таковые имеются, выборки превышают порог. Чтобы сгенерировать кривые при низком PFA, необходимо использовать несколько выборок по порядку обратных PFA. Этот тип симуляции занимает много времени. Следующая кривая использует один миллион импульсов вместо 45 000. Чтобы запустить эту симуляцию, повторите пример, но установите Npulse до 1000000.