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

Установите желаемую вероятность обнаружения быть 0.9 и вероятность ложного предупреждения, чтобы быть $$10^{-6}$$. Установите максимальную область значений радара к 4 000 метров и разрешения области значений 50 метров. Установите фактический целевой диапазон на 3 000 метров. Установите эффективную площадь рассеивания цели на 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 и фунт, требует обработки многих сигналов. Чтобы поддержать требования к памяти на низком уровне, обработайте сигналы во фрагментах импульсов. Определите номер импульсов к процессу к 45 000 и установите размер каждого фрагмента к 10 000.

Npulse = 45000;
Npulsebuffsize = 10000;

Выберите Waveform и Signal Parameters

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

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 требует, чтобы достаточная степень сигнала прибыла в приемник после того, как цель отражает сигнал. Вычислите минимальный ОСШ, должен был достигнуть заданной вероятности ложного предупреждения и вероятности обнаружения при помощи уравнения Albersheim.

snr_min = albersheim(pd,pfa);

Чтобы достигнуть этого ОСШ, достаточная степень должна быть передана к цели. Используйте основное уравнение радиолокации, чтобы оценить пиковую мощность передачи, peak_power, требуемый достигнуть заданного ОСШ в дБ для цели в области значений 3 000 метров. Полученный сигнал также зависит от эффективной площади рассеивания цели (ЭПР). который принят, чтобы следовать не колеблющейся модели (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);

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

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

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);

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

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

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

Вычислите гистограммы целевых существующих и целевых отсутствующих возвратов. Используйте 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 является константой по соседним индексам. Затем удалите эти значения из массивов PFA и Pd.

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.