Обнаружение сигнала по нескольким выборкам

Открыть пример в Engee

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

Введение

В примере "Обнаружение сигнала в белом гауссовском шуме" представлена базовая задача обнаружения сигнала. В этом примере для обнаружения используется только одна выборка принятого сигнала. В данном примере в процессе обнаружения используется больше образцов для повышения эффективности обнаружения.

Как и в предыдущем примере, предположим, что мощность сигнала равна 1, а отношение сигнал/шум (SNR) одного образца составляет 3 дБ. Число испытаний Монте-Карло равно 100000. Желаемый уровень вероятности ложной тревоги (Pfa) - 0,001.

Ntrial = Int64(1e5); # количество испытаний Монте-Карло 
Pfa = 1e-3; 
snrdb = 3; # SNR в дБ 
snr = db2pow(snrdb); # SNR в линейном масштабе 
npower = 1/snr; # мощность шума 
namp = sqrt(npower/2); # амплитуда шума в каждом канале

Обнаружение сигнала с помощью длинной формы волны

Как уже говорилось в предыдущем примере, порог определяется на основе Pfa. Поэтому, пока выбирается порог, Pfa остается неизменным, и наоборот. Между тем, конечно, желательно иметь более высокую вероятность обнаружения (Pd). Один из способов добиться этого - использовать несколько образцов для обнаружения. Например, в предыдущем случае SNR на одном образце составляет 3 дБ. Если можно использовать несколько образцов, то согласованный фильтр может дать дополнительный выигрыш в SNR и тем самым улучшить производительность. На практике для достижения такого выигрыша можно использовать более длинную форму сигнала. В случае обработки сигналов дискретного времени множественные выборки также могут быть получены путем увеличения частоты дискретизации.

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

Nsamp = 2;
wf = ones(Nsamp,1);
mf = conj.(reverse(wf)); # matched filter

Для когерентного приемника сигнал, шум и порог определяются следующим образом

using Random
# fix the random number generator
seed = 2014
rstream = Xoshiro(seed)

s = wf*ones(1,Ntrial);
n = namp*(randn(rstream,Nsamp,Ntrial)+1im*randn(rstream,Nsamp,Ntrial));
snrthreshold = db2pow(npwgnthresh(Pfa, 1,"coherent"));
mfgain = mf'*mf;
threshold = sqrt.(npower*mfgain*snrthreshold);   # Final threshold T;

Если цель присутствует:

x = s + n;
y = mf'*x;
z = real(y);
Pd = sum(z.>threshold)/Ntrial
print("Pd = $(Pd)")

Pd = 0.39814

Если цель отсутствует

x = n;
y = mf'*x;
z = real(y);
Pfa = sum(z.>threshold)/Ntrial
print("Pfa = $(Pfa)")

Pfa = 0.00104

Обратите внимание, что SNR улучшается благодаря согласованному фильтру.

snr_new = snr*mf'*mf;
snrdb_new = pow2db.(snr_new)[1];
print("snrdb_new = $(round(snrdb_new;sigdigits=5))")

snrdb_new = 6.0103

Plot the ROC curve with this new SNR value.

rocsnr(snrdb_new,SignalType="NonfluctuatingCoherent",MinPfa=1e-4)

Из рисунка видно, что точка, заданная Pfa и Pd, попадает прямо на кривую. Таким образом, SNR, соответствующий ROC-кривой, - это SNR одного образца на выходе согласованного фильтра. Это показывает, что, хотя для обнаружения можно использовать несколько образцов, пороговое значение SNR для одного образца (snrthreshold в программе) не меняется по сравнению со случаем простого образца. Изменений нет, потому что пороговое значение, по сути, определяется Pfa. Однако конечный порог, T, изменяется из-за дополнительного усиления согласованного фильтра. Результирующее значение Pfa осталось прежним по сравнению со случаем, когда для обнаружения используется только один образец. Однако дополнительное согласованное усиление улучшило Pd с 0,1390 до 0,3947.

Для проверки зависимости между Pd, Pfa и SNR можно провести аналогичные расчеты для некогерентного приемника.

Обнаружение сигнала с помощью интегрирования импульсов

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

Предположим, что интегрирование составляет 2 импульса. Затем постройте принятый сигнал и примените к нему согласованный фильтр.

PulseIntNum = 2;
Ntotal = PulseIntNum*Ntrial;
s = wf*exp.(1im*2*pi*rand(rstream,1,Ntotal));  # noncoherent
n = sqrt(npower/2).*(randn(rstream,Nsamp,Ntotal) + 1im*randn(rstream,Nsamp,Ntotal));

Если цель присутствует получаем

x = s .+ n;
y = mf'*x;
y = reshape(y,Ntrial,PulseIntNum);  # reshape to align pulses in columns;

Интегрировать импульсы можно с помощью одного из двух возможных подходов. Оба подхода связаны с аппроксимацией модифицированной функции Бесселя первого рода, которая встречается при моделировании теста отношения правдоподобия (LRT) процесса некогерентного обнаружения с использованием нескольких импульсов. Первый подход заключается в суммировании abs(y)^2 по импульсам, что часто называют детектором по квадратичному закону. Второй подход заключается в суммировании abs(y) по всем импульсам, что часто называют линейным детектором. При малом SNR предпочтительнее использовать детектор с квадратным законом, в то время как при большом SNR выгоднее использовать линейный детектор. В данном моделировании мы используем детектор с квадратным законом. Однако разница между двумя видами детекторов обычно находится в пределах 0,2 дБ.

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

\begin{align} y=\sqrt{|x_1|^2+\cdots+|x_n|^2}\ . \end{align}

z = pulsint(y,"noncoherent");

Связь между порогом T и Pfa, с учетом новой достаточной статистики z, определяется следующим образом

\begin{align} P_{fa}=1-I\left(\frac{T^2/(NM)}{\sqrt{L}},L-1\right) =1-I\left(\frac{\rm SNR}{\sqrt{L}},L-1\right)\ . \end{align}

где \begin{align} I(u,K)=\int_0^{u\sqrt{K+1}}\frac{e^{-\tau}\tau^K}{K!}d\tau \end{align}

форма неполной гамма-функции Пирсона, а L - количество импульсов, используемых для интегрирования импульсов. Используя детектор с квадратичным законом, можно вычислить пороговое значение SNR, связанное с интегрированием импульсов, используя функцию npwgnthresh, как и раньше.

snrthreshold = db2pow(npwgnthresh(Pfa,PulseIntNum,"noncoherent"));

Результирующий порог для достаточной статистики, $z$, определяется следующим образом

mfgain = mf'*mf;
threshold = (sqrt.(npower*mfgain*snrthreshold))[1];

Вероятность обнаружения определяется следующим образом

Pd = sum(z.>threshold)/Ntrial
print("Pd = $(Pd)")

Pd = 0.52807

Затем вычислите Pfa, когда принимаемый сигнал является только шумом, используя некогерентный детектор с интегрированными 2 импульсами.

x = n;
y = mf'*x;
y = reshape(y,Ntrial,PulseIntNum);
z = pulsint(y,"noncoherent");
Pfa = sum(z.>threshold)/Ntrial
print("Pfa = $(Pfa)")

Pfa = 0.00101

Чтобы построить ROC-кривую с импульсным интегрированием, необходимо указать количество импульсов, используемых при интегрировании, в функции rocsnr:

rocsnr(snrdb_new,SignalType="NonfluctuatingNoncoherent",
    MinPfa=1e-4,NumPulses=PulseIntNum)

И снова точка, заданная Pfa и Pd, попадает на кривую. Таким образом, SNR на кривой ROC определяет SNR одного образца, используемого для детектирования по одному импульсу.

Такое значение SNR также может быть получено из Pd и Pfa с помощью уравнения Альбершема. Результат, полученный из уравнения Альбершейма, является лишь приближением, но он достаточно хорош в часто используемых диапазонах Pfa, Pd и интегрирования импульсов.

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

Чтобы рассчитать необходимое SNR для одной выборки для достижения определенных Pd и Pfa, используйте функцию Альбершема

snr_required = albersheim(Pd,Pfa,N=PulseIntNum)
print("snr_required = $(round(snr_required;digits=5))")

snr_required = 6.01593

Это рассчитанное требуемое значение SNR соответствует новому значению SNR 4 дБ.

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

rocsnr(snrdb_new,SignalType="NonfluctuatingNoncoherent",
    MinPfa=1e-4,NumPulses=1)

Из рисунка видно, что без интеграции импульсов Pd может составлять только около 0,24 при Pfa 1e-3. При интеграции двух импульсов, как показано в приведенном выше моделировании Монте-Карло, для того же Pfa Pd составляет около 0,53.

Заключение

Этот пример показал, как использование множественной выборки сигнала при обнаружении может повысить вероятность обнаружения при сохранении желаемого уровня вероятности ложной тревоги. В частности, было показано, как использование более длинной формы сигнала или метода интеграции импульсов позволяет улучшить Pd. На примере показана связь между Pd, Pfa, ROC-кривой и уравнением Альберсхайма. Эффективность рассчитывается с помощью моделирования методом Монте-Карло.