Определение требований к радару

Установите требуемую вероятность обнаружения равной 0,9, а вероятность ложного аварийного сигнала равной $$\mathrm{}{\mathrm{}}^{\mathrm{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, необходимый для достижения заданной вероятности ложной тревоги и вероятности обнаружения, используя уравнение Олберсхайма.

snr_min = albersheim(pd,pfa);

Для достижения этого SNR на цель должна передаваться достаточная мощность. Используйте уравнение радара для оценки пиковой мощности передачи, peak_power, требуется для достижения заданного SNR в дБ для цели на дальности 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 К.

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 формы сигнала. Затем создайте соответствующий фильтр System object™.

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

Вычислить ячейку целевого диапазона

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

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-absent. Используйте 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 $$\mathrm{}{\mathrm{}}^{\mathrm{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. Затем постройте график теоретической кривой Пфа.

[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.