реклама на сайте
подробности

 
 
 
Reply to this topicStart new topic
> Некогерентный обнаружитель: практика применения
sqrt(2)
сообщение Dec 16 2016, 09:05
Сообщение #1


Участник
*

Группа: Участник
Сообщений: 58
Регистрация: 20-07-16
Пользователь №: 92 633



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

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

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

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

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

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

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

Сообщение отредактировал sqrt(2) - Dec 16 2016, 09:12
Эскизы прикрепленных изображений
Прикрепленное изображение
 
Go to the top of the page
 
+Quote Post
andyp
сообщение Dec 16 2016, 09:48
Сообщение #2


Местный
***

Группа: Участник
Сообщений: 399
Регистрация: 23-07-08
Пользователь №: 39 163



Обычно при проектировании обнаружителя закладываются на худшее возможное отношение сигнал/шум 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
Go to the top of the page
 
+Quote Post
sqrt(2)
сообщение Dec 16 2016, 11:11
Сообщение #3


Участник
*

Группа: Участник
Сообщений: 58
Регистрация: 20-07-16
Пользователь №: 92 633



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

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 (картинка к данному посту; маркером помечена регистрация сигнала на детекторе).

Что не так?



Эскизы прикрепленных изображений
Прикрепленное изображение
 
Go to the top of the page
 
+Quote Post
andyp
сообщение Dec 16 2016, 11:36
Сообщение #4


Местный
***

Группа: Участник
Сообщений: 399
Регистрация: 23-07-08
Пользователь №: 39 163



Цитата(sqrt(2) @ Dec 16 2016, 14:11) *
Что не так?


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

Относительный порог из формулы:
>> h = sqrt(-2*log(1e-6)*1*1/16)
h = 1.3141

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

Проблема скорее всего в том, что мощность полезного сигнала на выходе СФ не отнормирована к 1. Вобщем, нужно выход СФ на длину накопления поделить (ну или отсчеты ИХ).
Go to the top of the page
 
+Quote Post
sqrt(2)
сообщение Dec 16 2016, 12:42
Сообщение #5


Участник
*

Группа: Участник
Сообщений: 58
Регистрация: 20-07-16
Пользователь №: 92 633



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

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

Цитата(andyp @ Dec 16 2016, 14:36) *
Для Вашего случая
sigma_n^2 = 1/16

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

Сообщение отредактировал sqrt(2) - Dec 16 2016, 12:43
Go to the top of the page
 
+Quote Post
andyp
сообщение Dec 16 2016, 13:34
Сообщение #6


Местный
***

Группа: Участник
Сообщений: 399
Регистрация: 23-07-08
Пользователь №: 39 163



Цитата(sqrt(2) @ Dec 16 2016, 15:42) *
А так получается очень большой порог. На выходе детектора получается 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
Go to the top of the page
 
+Quote Post
sqrt(2)
сообщение Dec 16 2016, 14:10
Сообщение #7


Участник
*

Группа: Участник
Сообщений: 58
Регистрация: 20-07-16
Пользователь №: 92 633



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

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

Цитата(andyp @ Dec 16 2016, 16:34) *
Ну и да, стоит убедиться еще в одной вещи - 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
Go to the top of the page
 
+Quote Post
Милливольт
сообщение Dec 16 2016, 18:31
Сообщение #8


Участник
*

Группа: Участник
Сообщений: 65
Регистрация: 17-05-15
Пользователь №: 86 729



Цитата(sqrt(2) @ Dec 16 2016, 15:10) *
UPD: рано радовался. Пока величина 1/sigma_n^2 не падает резко, то работает. Потом перестает. И для более-менее большой амплитуды тоже не работает.


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

Прикрепленное изображение


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

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

Сообщение отредактировал Милливольт - Dec 16 2016, 18:31
Go to the top of the page
 
+Quote Post
sqrt(2)
сообщение Dec 17 2016, 16:15
Сообщение #9


Участник
*

Группа: Участник
Сообщений: 58
Регистрация: 20-07-16
Пользователь №: 92 633



Пишет "страница не найдена"
Go to the top of the page
 
+Quote Post
Милливольт
сообщение Dec 17 2016, 16:36
Сообщение #10


Участник
*

Группа: Участник
Сообщений: 65
Регистрация: 17-05-15
Пользователь №: 86 729



Цитата(sqrt(2) @ Dec 17 2016, 17:15) *
Пишет "страница не найдена"

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

Сообщение отредактировал Милливольт - Dec 17 2016, 16:39
Go to the top of the page
 
+Quote Post
andyp
сообщение Dec 19 2016, 12:14
Сообщение #11


Местный
***

Группа: Участник
Сообщений: 399
Регистрация: 23-07-08
Пользователь №: 39 163



Цитата(sqrt(2) @ Dec 16 2016, 17:10) *
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
Go to the top of the page
 
+Quote Post

Reply to this topicStart new topic
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0

 


RSS Текстовая версия Сейчас: 30th April 2017 - 14:53
Рейтинг@Mail.ru


Страница сгенерированна за 0.0153 секунд с 7
ELECTRONIX ©2004-2016