sqrt(2) 0 16 декабря, 2016 Опубликовано 16 декабря, 2016 (изменено) · Жалоба Здравствуйте. Есть задача обнаружения сигнала известного по форме, но его фаза и амплитуда - неизвестны. Как я понимаю, в таком случае мне нужен некогерентный приемник. Отсчеты поступают на два согласованных фильтра - у одного импульсная характеристика это зеркально отраженный во времени сигнал, а у второго - мнимая часть от преобразования Гильберта от импульсной характеристики первого фильтра. Далее выходы обоих фильтров возводятся в квадрат, вычисляется сумма квадратов выходов СФ и от всего этого берется квадратный корень - т.е., детектор огибающей. И вот тут начинается непонятное. По идее, надо как раз величину этого корня сравнивать с неким порогом. Поскольку в качестве критерия оптимальности используется критерий Неймана-Пирсона, то порог этот определяется исходя из заданной вероятности ложной тревоги (см. картинку). вероятность ложной тревоги = exp(-(нормированный порог)^2/2). Как найти нормированный порог - проблем нет (он на картинке обозначен как h0). Вопрос в том, как найти именно порог обнаружения h. Не совсем понимаю, как зная SNR, структуры фильтров, разрядности и тд (особенности архитектуры приемника у себя) определить то, что на картинке обозначено как No и E1. Видимо нужна еще какая-то априорная информация, но скорее всего я не понимаю какую-то простую вещь. N0 - очевидно, мощность шума в полосе сигнала. Но вот как посчитать энергию, да еще с учетом некой усредненной амплитуды... Изменено 16 декабря, 2016 пользователем sqrt(2) Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
andyp 9 16 декабря, 2016 Опубликовано 16 декабря, 2016 · Жалоба Обычно при проектировании обнаружителя закладываются на худшее возможное отношение сигнал/шум 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 Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
sqrt(2) 0 16 декабря, 2016 Опубликовано 16 декабря, 2016 · Жалоба Хм, возможно я что-не понимаю, но.... 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 (картинка к данному посту; маркером помечена регистрация сигнала на детекторе). Что не так? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
andyp 9 16 декабря, 2016 Опубликовано 16 декабря, 2016 · Жалоба Что не так? Для Вашего случая sigma_n^2 = 1/16 E = 1 A2 = 1+1/16; Относительный порог из формулы: >> h = sqrt(-2*log(1e-6)*1*1/16) h = 1.3141 (следите за руками, тут мог двойку потерять) Проблема скорее всего в том, что мощность полезного сигнала на выходе СФ не отнормирована к 1. Вобщем, нужно выход СФ на длину накопления поделить (ну или отсчеты ИХ). Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
sqrt(2) 0 16 декабря, 2016 Опубликовано 16 декабря, 2016 (изменено) · Жалоба А так получается очень большой порог. На выходе детектора получается 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 - не амплитуда, а энергетическая величина. Изменено 16 декабря, 2016 пользователем sqrt(2) Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
andyp 9 16 декабря, 2016 Опубликовано 16 декабря, 2016 (изменено) · Жалоба А так получается очень большой порог. На выходе детектора получается 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 всегда путаюсь. Е и у меня - энергетическая величина. Изменено 16 декабря, 2016 пользователем andyp Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
sqrt(2) 0 16 декабря, 2016 Опубликовано 16 декабря, 2016 (изменено) · Жалоба Формула та же, просто сомножители под корнем по другому сгруппированы. Для нормировки проще всего подать на вход повторения Вашего 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 не падает резко, то работает. Потом перестает. И для более-менее большой амплитуды тоже не работает. Изменено 16 декабря, 2016 пользователем sqrt(2) Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Милливольт 0 16 декабря, 2016 Опубликовано 16 декабря, 2016 (изменено) · Жалоба UPD: рано радовался. Пока величина 1/sigma_n^2 не падает резко, то работает. Потом перестает. И для более-менее большой амплитуды тоже не работает. С Вашего позволения: у нас с порогами (правда, для другого сигнала) были примерно такие же сложности. Описание здесь: http://www.tredex-company.com/ru/content/о...рного-излучения Может быть пригодится. Изменено 16 декабря, 2016 пользователем Милливольт Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
sqrt(2) 0 17 декабря, 2016 Опубликовано 17 декабря, 2016 · Жалоба Пишет "страница не найдена" Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Милливольт 0 17 декабря, 2016 Опубликовано 17 декабря, 2016 (изменено) · Жалоба Пишет "страница не найдена" Да, извините, слишком длинное название для этого движка. На главной странице сайта http://www.tredex-company.com - справа статьи, эта первая "Опыт проектирования и эксплуатации устройства обнаружения лазерного излучения." Изменено 17 декабря, 2016 пользователем Милливольт Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
andyp 9 19 декабря, 2016 Опубликовано 19 декабря, 2016 (изменено) · Жалоба 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; Вот что получается при его выполнении Изменено 19 декабря, 2016 пользователем andyp Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться