Здравствуйте. Есть задача обнаружения сигнала известного по форме, но его фаза и амплитуда - неизвестны. Как я понимаю, в таком случае мне нужен некогерентный приемник.
Отсчеты поступают на два согласованных фильтра - у одного импульсная характеристика это зеркально отраженный во времени сигнал, а у второго - мнимая часть от преобразования Гильберта от импульсной характеристики первого фильтра.
Далее выходы обоих фильтров возводятся в квадрат, вычисляется сумма квадратов выходов СФ и от всего этого берется квадратный корень - т.е., детектор огибающей.
И вот тут начинается непонятное. По идее, надо как раз величину этого корня сравнивать с неким порогом. Поскольку в качестве критерия оптимальности используется критерий Неймана-Пирсона, то порог этот определяется исходя из заданной вероятности ложной тревоги (см. картинку).
вероятность ложной тревоги = exp(-(нормированный порог)^2/2).
Как найти нормированный порог - проблем нет (он на картинке обозначен как h0). Вопрос в том, как найти именно порог обнаружения h. Не совсем понимаю, как зная SNR, структуры фильтров, разрядности и тд (особенности архитектуры приемника у себя) определить то, что на картинке обозначено как No и E1. Видимо нужна еще какая-то априорная информация, но скорее всего я не понимаю какую-то простую вещь.
N0 - очевидно, мощность шума в полосе сигнала. Но вот как посчитать энергию, да еще с учетом некой усредненной амплитуды...
Полная версия этой страницы: Некогерентный обнаружитель: практика применения
Форум разработчиков электроники ELECTRONIX.ru > Цифровая обработка сигналов - ЦОС (DSP) > Алгоритмы ЦОС (DSP)
Обычно при проектировании обнаружителя закладываются на худшее возможное отношение сигнал/шум 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
Т.е. при вычислении порога используют два уравнения:
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 (картинка к данному посту; маркером помечена регистрация сигнала на детекторе).
Что не так?
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 (картинка к данному посту; маркером помечена регистрация сигнала на детекторе).
Что не так?
Цитата(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. Вобщем, нужно выход СФ на длину накопления поделить (ну или отсчеты ИХ).
А так получается очень большой порог. На выходе детектора получается 0.33 , а порог - 1.31.
Нашел вот что: http://www.mathworks.com/help/phased/examp...l#zmw57dd0e7431
Вроде о том же самом, но формула в итоге другая (или это я считаю её другой): threshold = sqrt(npower*mfgain*snrthreshold);
snrthreshold - это наш нормированный порог h0.
mfgain - усиление, которое дает согласованный фильтр (как я понял, этот множитель здесь для того, чтобы не нормировать выход СФ)
npower - мощность шума.
В целом, вроде похоже, но не очень.
Кстати, а почему так? Я всю жизнь считал, что дисперсия шума (во всяком случае дисперсия AWGN, на счет остальных шумов не уверен) = спектральная плотность мощности/2. В отечественной литературе принято обозначение для этого дела N0/2, а SNR = 2E/N0, но здесь E - не амплитуда, а энергетическая величина.
Код
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);
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
sigma_n^2 = 1/16
Кстати, а почему так? Я всю жизнь считал, что дисперсия шума (во всяком случае дисперсия AWGN, на счет остальных шумов не уверен) = спектральная плотность мощности/2. В отечественной литературе принято обозначение для этого дела N0/2, а SNR = 2E/N0, но здесь E - не амплитуда, а энергетическая величина.
Цитата(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 - мощность шума.
В целом, вроде похоже, но не очень.
Вроде о том же самом, но формула в итоге другая (или это я считаю её другой): 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, 16:34)
Формула та же, просто сомножители под корнем по другому сгруппированы. Для нормировки проще всего подать на вход повторения Вашего header и подобрать нормировочный множитель так, чтобы максимальные пики на выходе детектора были равны 1.
Ну и да, стоит убедиться еще в одной вещи - mean(h)~=mean(h_hilbert)~=0.
Ну и да, стоит убедиться еще в одной вещи - 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);
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, 15:10)
UPD: рано радовался. Пока величина 1/sigma_n^2 не падает резко, то работает. Потом перестает. И для более-менее большой амплитуды тоже не работает.
С Вашего позволения: у нас с порогами (правда, для другого сигнала) были примерно такие же сложности.
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Описание здесь: http://www.tredex-company.com/ru/content/о...рного-излучения
Может быть пригодится.
Пишет "страница не найдена"
Цитата(sqrt(2) @ Dec 17 2016, 17:15)
Пишет "страница не найдена"
Да, извините, слишком длинное название для этого движка. На главной странице сайта http://www.tredex-company.com - справа статьи, эта первая "Опыт проектирования и эксплуатации устройства обнаружения лазерного излучения."
Цитата(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'])
%в двух случаях - известной и неизвестной амплитуды сигнала на входе
%входные параметры (предполагаем использование кода Баркера с длиной 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;
Вот что получается при его выполнении
Нажмите для просмотра прикрепленного файла
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.