Некогерентный обнаружитель: практика применения Опции

[sqrt(2)]

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

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

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

И вот тут начинается непонятное. По идее, надо как раз величину этого корня сравнивать с неким порогом. Поскольку в качестве критерия оптимальности используется критерий Неймана-Пирсона, то порог этот определяется исходя из заданной вероятности ложной тревоги (см. картинку).

вероятность ложной тревоги = exp(-(нормированный порог)^2/2).

Как найти нормированный порог - проблем нет (он на картинке обозначен как h0). Вопрос в том, как найти именно порог обнаружения h. Не совсем понимаю, как зная SNR, структуры фильтров, разрядности и тд (особенности архитектуры приемника у себя) определить то, что на картинке обозначено как No и E1. Видимо нужна еще какая-то априорная информация, но скорее всего я не понимаю какую-то простую вещь.

N0 - очевидно, мощность шума в полосе сигнала. Но вот как посчитать энергию, да еще с учетом некой усредненной амплитуды...

Сообщение отредактировал sqrt(2) - Dec 16 2016, 09:12

Эскизы прикрепленных изображений

 

[sqrt(2)]

А так получается очень большой порог. На выходе детектора получается 0.33 , а порог - 1.31.

Код

h = flipud(header);  %ИХ фильтра
h_hilbert = hilbert(h);
Y = length(h);
%прохождение сигнала через согласованный фильтр
mf_origin = filter(h,1,input_signal);
mf_origin = mf_origin./Y;
mf_hilbert = filter(imag(h_hilbert),1,input_signal);
mf_hilbert = mf_hilbert./Y;
%детектор
A = sqrt(mf_origin.^2 + mf_hilbert.^2);

Нашел вот что: http://www.mathworks.com/help/phased/examp...l#zmw57dd0e7431

Вроде о том же самом, но формула в итоге другая (или это я считаю её другой): threshold = sqrt(npower*mfgain*snrthreshold);

snrthreshold - это наш нормированный порог h0.

mfgain - усиление, которое дает согласованный фильтр (как я понял, этот множитель здесь для того, чтобы не нормировать выход СФ)

npower - мощность шума.

В целом, вроде похоже, но не очень.

Для Вашего случая
sigma_n^2 = 1/16

Кстати, а почему так? Я всю жизнь считал, что дисперсия шума (во всяком случае дисперсия AWGN, на счет остальных шумов не уверен) = спектральная плотность мощности/2. В отечественной литературе принято обозначение для этого дела N0/2, а SNR = 2E/N0, но здесь E - не амплитуда, а энергетическая величина.

Сообщение отредактировал sqrt(2) - Dec 16 2016, 12:43

[andyp]

А так получается очень большой порог. На выходе детектора получается 0.33 , а порог - 1.31.

Нашел вот что: http://www.mathworks.com/help/phased/examp...l#zmw57dd0e7431

Вроде о том же самом, но формула в итоге другая (или это я считаю её другой): threshold = sqrt(npower*mfgain*snrthreshold);

snrthreshold - это наш нормированный порог h0.

mfgain - усиление, которое дает согласованный фильтр (как я понял, этот множитель здесь для того, чтобы не нормировать выход СФ)

npower - мощность шума.

В целом, вроде похоже, но не очень.

Формула та же, просто сомножители под корнем по другому сгруппированы. Для нормировки проще всего подать на вход повторения Вашего header и подобрать нормировочный множитель так, чтобы максимальные пики на выходе детектора были равны 1.
Ну и да, стоит убедиться еще в одной вещи - mean(h)~=mean(h_hilbert)~=0.

Другими словами, проверить, что на выходе Вы имеете распределение Релея в случае отсутствия полезного сигнала и Райса в пиках при его наличии.

Кстати, а почему так? Я всю жизнь считал, что дисперсия шума (во всяком случае дисперсия AWGN, на счет остальных шумов не уверен) = спектральная плотность мощности/2. В отечественной литературе принято обозначение для этого дела N0/2, а SNR = 2E/N0, но здесь E - не амплитуда, а энергетическая величина.

Это все игры вокруг односторонней и двусторонней СПМ. Односторонняя вроде по определению равна удвоенной двусторонней для положительных частот. N_0 - обычно односторонняя СПМ. Я с этими 2 всегда путаюсь. Е и у меня - энергетическая величина.

Сообщение отредактировал andyp - Dec 16 2016, 13:35

[sqrt(2)]

Формула та же, просто сомножители под корнем по другому сгруппированы. Для нормировки проще всего подать на вход повторения Вашего header и подобрать нормировочный множитель так, чтобы максимальные пики на выходе детектора были равны 1.
Ну и да, стоит убедиться еще в одной вещи - mean(h)~=mean(h_hilbert)~=0.

Ну допустим пики будут иметь максимальное значение 1, но ведь порог получился 1.31.

Ну и да, стоит убедиться еще в одной вещи - mean(h)~=mean(h_hilbert)~=0.

mean(h) ~= 0
mean(h_hilbert) == 0

Короче, похоже дело вот в чем.

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

В связи с этим переписал код так:

Код

h = flipud(header);  %ИХ фильтра
h_hilbert = hilbert(h);
Y = length(h);
%прохождение сигнала через согласованный фильтр
mf_origin = filter(h,1,input_signal);
mf_origin = 4.*(mf_origin./Y);
mf_hilbert = filter(imag(h_hilbert),1,input_signal);
mf_hilbert = 4.*(mf_hilbert./Y);
%детектор
A = sqrt(mf_origin.^2 + mf_hilbert.^2);

С указанным вами ранее способом вычисления порога - заработало.

UPD: рано радовался. Пока величина 1/sigma_n^2 не падает резко, то работает. Потом перестает. И для более-менее большой амплитуды тоже не работает.

Сообщение отредактировал sqrt(2) - Dec 16 2016, 14:37

[andyp]

UPD: рано радовался. Пока величина 1/sigma_n^2 не падает резко, то работает. Потом перестает. И для более-менее большой амплитуды тоже не работает.

Непонятно, что у Вас не получается

Давайте сначала разберемся с тем, что у Вас очень большой порог получается для заданного ОСШ и вероятности ложной тревоги. Это связано с тем, что чем меньше вероятность ложной тревоги, тем больше получается вероятность пропуска сигнала, поэтому минимальное рабочее отношение сигнал-шум нужно брать не от балды, а с учетом вероятности пропуска преамбулы.

Обычно строят семейство т. н. ROC curves (можно найти в wiki) для разных значений ОСШ на входе согласованного фильтра и по ним выбирают тройку D, F, ОСШ

Вот скрипт как это сделать для разных входных ОСШ

Код

%строим кривые ROC для заданного Eb/N0 на выходе согласованного фильтра
%в двух случаях - известной и неизвестной амплитуды сигнала на входе

%входные параметры (предполагаем использование кода Баркера с длиной 11)
Preamble_len = 11;
SNRin_dB  = 3;
Pin = 1; %средняя мощность входного сигнала


E = Pin*Preamble_len;
%на входе шум с односторонней спектральной плотностью мощности N0
%(диспресия равна N0/2)
N0 = 2*10^(-SNRin_dB/10) * Pin;

%возможные вероятности ложной тревоги
p = -6:0;
FA_vals = 10.^p - 1e-10;

%в случае отсутствия сигнала на входе, на выходе СФ и детектора имеем
%распределение Релея с параметром sigma^2 = N0*E;
%находим все возможные значения порога для рассчета вероятности правильного
%обнаружения сигнала
h = icdf('rayl',1 - FA_vals,sqrt(N0*E/2));

%рассчитываем для них вероятности правильного обнаружения
D = 1 - cdf('rician',h, E,sqrt(N0*E/2));
D_unk_amp = 1 - cdf('rayl',h, sqrt((N0*E + E^2)/2));
%вот и кривая ROC
figure(1)
plot(log10(FA_vals),D)
hold on
plot(log10(FA_vals),D_unk_amp,'r')
hold off
xlabel('log10(F)');
ylabel('D');
grid on;
title(['ROC curve for ' num2str(SNRin_dB) ' dB input SNR'])

Вот что выходит



Вот простенький скрипт для тестирования рассчитанных порогов

Код

%---------Входные параметры--------------
%мощность полезного сигнала
Pin = 1;
%входной рабочий snr  (отношение  мощностей сигнала и шума)
snr_dB = 10;
%расчетный входной snr обнаружителя (худшее отношение  мощностей сигнала и шума)
targetSNR_dB = 3;
%порог, рассчитанный предыдущим скриптом для targetSNR_dB = 3, Pfa = 1e-4
th = 7.1258

%спм шума на входе
N0 = 2*10^(-snr_dB/10) * Pin;

%Полезный сигнал с фазовым сдвигом
signal = sqrt(Pin) * repmat(barker,1,5)  *exp(j*pi/4);
signal = [zeros(1, length(barker))  signal  zeros(1, length(barker))];
%белый шум c двусторонней спм N0/2
noise = sqrt(N0/2) *(randn(size(signal)) + j*randn(size(signal)));
%сигнал плюс шум
signal = signal + noise;


out = filter(fliplr(barker),1, signal);
plot(abs(out))
hold on;
plot(th*ones(size(out)),'r')
hold off;

Вот что получается при его выполнении



Сообщение отредактировал andyp - Dec 19 2016, 12:16