Перейти к содержанию
    

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

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

 

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

 

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

 

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

 

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

 

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

 

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

post-92633-1481878959_thumb.png

Изменено пользователем sqrt(2)

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Обычно при проектировании обнаружителя закладываются на худшее возможное отношение сигнал/шум A1 = (E/sigma_n^2) на входе, которое используют для вычисления порога.

Т.е. при вычислении порога используют два уравнения:

 

E/sigma_n^2 = A1 (1)

E+sigma_n^2 = A2 (2)

 

Если на входе есть АРУ, то эта система поддерживает A2 постоянным и известным, если нет - то параллельно с обнаружителем ставят схему оценивания A2 (это квадратор и усреднение).

 

Из двух уравнений находят E, sigma_n^2 и требуемый порог.

 

Разумеется (2), можно упростить для случаев высоких и низких рабочих ОСШ на входе. Например, для низких ОСШ (2) фактически станет

sigma_n^2 = A2

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Хм, возможно я что-не понимаю, но....

 

sigma_n^2 = N0/2 - если в обозначениях с картинки в первом посте.

 

Задано:

1) SNR = 12 dB (в разах - примерно 16);

2) вероятность ложной тревоги F = 1e-6.

 

Тогда, в тех же обозначениях с картинки:

 

h0 = sqrt[ -2*ln(F)]

 

В обозначениях из второго поста:

 

A1 = SNR = 16 => E/sigma_n^2 = 16

 

Допустим, что у меня сигнал единичной амлитуды (в среднем) и я не хочу, чтобы было выше 1.5. Тогда:

E + sigma_n^2 = 1.5

 

Тогда:

sigma_n^2 = 0.0938

E = 1.5

 

Отсюда порог h = sqrt[-2*ln(F)*E*sigma_n^2] = 2.

 

При моделировании вижу, что на деле мне нужен порог ориентировочно h >= 60 (картинка к данному посту; маркером помечена регистрация сигнала на детекторе).

 

Что не так?

 

 

post-92633-1481886683_thumb.png

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Что не так?

 

Для Вашего случая

sigma_n^2 = 1/16

E = 1

A2 = 1+1/16;

 

Относительный порог из формулы:

>> h = sqrt(-2*log(1e-6)*1*1/16)

h = 1.3141

 

(следите за руками, тут мог двойку потерять)

 

Проблема скорее всего в том, что мощность полезного сигнала на выходе СФ не отнормирована к 1. Вобщем, нужно выход СФ на длину накопления поделить (ну или отсчеты ИХ).

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

А так получается очень большой порог. На выходе детектора получается 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)

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

А так получается очень большой порог. На выходе детектора получается 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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Формула та же, просто сомножители под корнем по другому сгруппированы. Для нормировки проще всего подать на вход повторения Вашего 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)

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

 

С Вашего позволения: у нас с порогами (правда, для другого сигнала) были примерно такие же сложности.

post-86729-1481912848_thumb.png

post-86729-1481912875_thumb.png

 

Описание здесь: http://www.tredex-company.com/ru/content/о...рного-излучения

 

Может быть пригодится.

Изменено пользователем Милливольт

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Пишет "страница не найдена"

Да, извините, слишком длинное название для этого движка. На главной странице сайта http://www.tredex-company.com - справа статьи, эта первая "Опыт проектирования и эксплуатации устройства обнаружения лазерного излучения."

Изменено пользователем Милливольт

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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'])

 

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

 

post-39163-1482149353_thumb.jpg

 

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

 

 
%---------Входные параметры--------------
%мощность полезного сигнала
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;

 

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

 

post-39163-1482149342_thumb.jpg

Изменено пользователем andyp

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

Вы можете написать сейчас и зарегистрироваться позже. Если у вас есть аккаунт, авторизуйтесь, чтобы опубликовать от имени своего аккаунта.

Гость
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...