Analiza ukŕadˇw detekcji ze stabilizacjĐ poziomu faŕszywego alarmu

Politechnika Warszawska

Wydział Elektroniki i Technik Informacyjnych

Instytut Systemów Elektronicznych

PRACA MAGISTERSKA

Analiza układów detekcji ze stabilizacją poziomu

fałszywego alarmu

Autor:

Igor Kulkowski

Kierownik pracy:

dr inż. Krzysztof Kulpa

Warszawa 2000

SPIS TREŚCI

WSTĘP......................................................................................................................... 4

1. PODSTAWY I NAJWAŻNIEJSZE POJĘCIA TEORII DETEKCJI............................ 6

1.1. PARAMETRY OKREŚLAJĄCE JAKOŚĆ DETEKCJI ........................................ 6

1.2. POJĘCIE RYZYKA. OPTYMALNE KRYTERIA DECYZYJNE ........................... 9

1.3. PRZYKŁAD OPTYMALNEJ DETEKCJI........................................................... 14

1.4. WERYFIKACJA HIPOTEZ ZA POMOCĄ POMIARÓW WIELOKROTNYCH ... 21

1.5. PRZYKŁAD OPTYMALNEJ DETEKCJI NA PODSTAWIE WIELU POMIARÓW

........................................................................................................................ 25

2. ANALIZA UKŁADÓW OPTYMALNEJ DETEKCJI RZECZYWISTYCH SYGNAŁÓW

RADAROWYCH .................................................................................................... 30

2.1. DETEKCJA SYGNAŁU O NIEZNANYCH PARAMETRACH............................ 30

2.2. STOSUNEK WIARYGODNOŚCI DLA SYGNAŁU O NIEZNANYCH

PARAMETRACH............................................................................................. 34

2.3. STOSUNEK WIARYGODNOŚCI DLA SYGNAŁU O NIEZNANEJ FAZIE........ 36

2.4. ZASTOSOWANIE FILTRU DOPASOWANEGO W DETEKCJI SYGNAŁÓW Z

NIEZNANYM CZASEM OPÓŹNIENIA ............................................................ 40

2.5. ZASTOSOWANIE ZESPOLONYCH AMPLITUD W ANALIZIE PROCESU

DETEKCJI....................................................................................................... 42

3. INTEGRACJA IMPULSÓW I STABILIZACJA POZIOMU FAŁSZYWEGO ALARMU

.............................................................................................................................. 54

3.1. ZASADA DZIAŁANIA RZECZYWISTYCH SYSTEMÓW RADIOLOKACYJNYCH

........................................................................................................................ 54

3.2. STOSUNEK WIARYGODNOŚCI DLA SYGNAŁU W POSTACI PACZKI

IMPULSÓW W.CZ. O NIEZNANYCH FAZACH............................................... 59

3.3. DETEKTOR DLA SYGNAŁU W POSTACI PACZKI IMPULSÓW W.CZ. O

NIEZNANYCH FAZACH.................................................................................. 61

3.4. INTEGRATOR BINARNY................................................................................ 65

3.5. FAŁSZYWY ALARM........................................................................................ 66

4. ANALIZA JAKOŚCI UKŁADÓW OPTYMALNEJ DETEKCJI SYGNAŁÓW

RADIOLOKACYJNYCH ........................................................................................ 72

4.1. ROKŁADY OBWIEDNI SYGNAŁU .................................................................. 73

4.2. ROZKŁADY SYGNAŁU ZA UKŁADEM CFAR................................................. 80

4.3. ROZKŁADY SYGNAŁU ZA INTEGRATOREM BINARNYM ............................ 85

4.4. ROZKŁADY SYGNAŁU ZA INTEGRATOREM DZIAŁAJĄCYM

BEZPOŚREDNIO NA AMPLITUDZIE ............................................................. 90

DODATKI ................................................................................................................... 95

A: OPIS ZESTAWU FUNKCJI SŁUŻĄCYCH DO ANALIZY DETEKTORÓW

CYFROWYCH................................................................................................. 95

binint........................................................................................................... 96
dbloss......................................................................................................... 98
detect ....................................................................................................... 100
normaliz.................................................................................................... 102
pdftocdf .................................................................................................... 104
ricepdf ...................................................................................................... 106
sqrtpdf ...................................................................................................... 107
sumpdf ..................................................................................................... 109

B: KODY

ŹRÓDŁOWE FUNKCJI SŁUŻĄCYCH DO ANALIZY DETEKTORÓW

CYFROWYCH............................................................................................... 111

binint......................................................................................................... 111
dbloss....................................................................................................... 112
detect ....................................................................................................... 113
normaliz.................................................................................................... 114
pdftocdf .................................................................................................... 115
ricepdf ...................................................................................................... 116
sqrtpdf ...................................................................................................... 116
sumpdf ..................................................................................................... 117

PODSUMOWANIE ................................................................................................... 119

LITERATURA........................................................................................................... 120

WSTĘP

Minęło ponad pięćdziesiąt lat od skonstruowania pierwszego radaru. Od

tamtej chwili urządzenia te zostały bardzo unowocześnione. Obecnie operator

nie musi się już bez przerwy wpatrywać w ekran obrazujący przestrzeń dookoła

stacji czy też w przebieg sygnału echa i na ich podstawie samemu podejmować

decyzję o obecności celu. Dzięki zastosowaniu metod cyfrowej obróbki sygnału

możliwe stało się zwiększenie prawdopodobieństwa wykrycia obiektów, a co za

tym idzie, można w większości przypadków odciążyć operatora i zrzucić

podejmowanie decyzji na maszynę. Układy detekcyjne składają się z wielu

bloków odpowiednio przekształcających sygnał echa tak, aby wykrycie było

możliwe nawet w niesprzyjających warunkach. Jednak taka duża komplikacja

systemów powoduje, że na etapie projektowania ich analiza staje się bardzo

żmudna i czasochłonna

−

zwłaszcza przy stosowaniu metod analitycznych. Z

pomocą mogą tu przyjść komputery. Stosując metody numeryczne można

znacznie skrócić czas potrzebny na przebadanie układu detektora.

Głównym celem niniejszej pracy jest stworzenie narzędzi w postaci

zestawu programów, które ułatwiałyby analizę detektorów cyfrowych.

Pozwalałyby one na określenie rozkładów prawdopodobieństw sygnału echa

wraz z zakłóceniami zarówno na wyjściu anteny jak i po przejściu przez

poszczególne bloki obróbki sygnału. Czynnością końcową byłoby wykreślenie

dla konkretnego zestawu bloków funkcjonalnych charakterystyk detekcji. Celem

drugoplanowym, ale wynikającym bezpośrednio z celu głównego, jest

przeprowadzenie analizy statystycznej cyfrowych układów detekcji oraz podanie

metod pozwalających na wyznaczenie odpowiednich rozkładów

prawdopodobieństwa.

Niniejsza praca składa się z czterech rozdziałów, z czego trzy pierwsze

zawierają analizę teoretyczną problemu, natomiast w rozdziale czwartym

przeprowadzono analizę jakości układów detekcji, wykorzystując jako narzędzia

stworzone w ramach realizowanej pracy dyplomowej zestawy programów,

pracujące w środowisku MATLAB.

W rozdziale pierwszym podano najważniejsze pojęcia teorii wykrywania,

zaprezentowano wybrane kryteria decyzyjne oraz omówiono problem

optymalnej detekcji pojedynczych wartości oraz sygnałów całkowicie znanych.

Przytoczono także prosty przykład optymalnej detekcji.

Rozdział drugi porusza problem wykrywania sygnałów w sytuacji, gdy

część z jego parametrów jest nieznana i może zmieniać się losowo.

Przeprowadzono teoretyczne rozważania, których wynikiem jest schemat

optymalnego układu detekcji z podziałem na bloki funkcjonalne, wykonujące

różne operacje na sygnale.

W trzecim rozdziale omówiono problem stabilizacji poziomu

prawdopodobieństwa fałszywego alarmu oraz integracji wielu impulsów, czyli

podejmowania decyzji w oparciu o dużą liczbę odbieranych sygnałów. W tym

miejscu przedstawiono podstawowe typowe rozwiązania detektorów oraz

sposób ich działania.

Omówione trzy rozdziały stanowią swego rodzaju próbę

uporządkowanego i uszeregowanego opisu rozważań dotyczących problemów

detekcji sygnałów, przedstawianych w polskojęzycznej literaturze w sposób

wyrywkowy lub pobieżny.

Rozdział czwarty zawiera analizę jakości układów detekcji. W rozdziale

tym wyliczono rozkłady prawdopodobieństw sygnałów na wyjściu

poszczególnych bloków funkcjonalnych całego układu optymalnej detekcji oraz

wykreślono dla nich charakterystyki wykrywania. Na końcu pracy dyplomowej,

w dodatkach A i B, znajdują się dokładny opis do stworzonych narzędzi oraz

kod źródłowy.

Napisane funkcje, pracujące w środowisku MATLAB, pozwoliły na

numeryczne wyznaczenie rozkładów prawdopodobieństw, oraz na wykreślenie

zamieszczonych w rozdziale czwartym charakterystyk detekcji dla różnych

typów detektorów.

Przedstawione w pracy autorskie narzędzia do analizy całych układów

detekcji oraz ich poszczególnych bloków stanowią uwieńczenie rozważań

teoretycznych i wydaje się, że mogą być przydatne w procesie dydaktycznym, a

także w pracach naukowo-badawczych.

1. PODSTAWY I NAJWAŻNIEJSZE POJĘCIA

TEORII DETEKCJI

1.1. PARAMETRY OKREŚLAJĄCE JAKOŚĆ DETEKCJI

Celem radiolokacji jest uzyskanie informacji na temat obiektów latających

znajdujących się w przestrzeni powietrznej oraz ich współrzędnych. Informacje

te uzyskuje się poprzez sondowanie przestrzeni sygnałami radiowymi. Jeśli

gdzieś w powietrzu znajduje się samolot, to sygnał odbije się od niego i wróci

do stacji nadawczo-odbiorczej. Na wejściu odbiornika poza sygnałem odbitym

od rzeczywistego obiektu znajdują się także zakłócenia. Są to np. sygnały

pochodzące od innych stacji radiolokacyjnych, odbite od powierzchni wody,

gęstych chmur czy nieruchomych obiektów na ziemi. Powodują one błędy w

wykrywaniu obiektów i określaniu ich współrzędnych. Przypadkowość sygnałów

radiolokacyjnych i zakłóceń powoduje, że przy analizie działania radaru

konieczne jest zastosowanie metod statystycznych i używanie do oceny jego

jakości parametrów statystycznych.

Analiza sygnałów radiolokacyjnych docierających z przestrzeni

powietrznej powinna być zakończona podjęciem decyzji na temat obecności lub

nieobecności obiektu w określonym miejscu przeszukiwanej przez radar

przestrzeni. Decyzja może być podjęta przy spełnieniu dwóch wzajemnie

wykluczających się warunków:

warunek H

−

"obiekt jest"

warunek H

−

"obiektu nie ma".

W trakcie podejmowania decyzji nie wiadomo, który z warunków jest w

rzeczywistości spełniony. Stawiane są dwie hipotezy co do obecności celu.

Celem jest ustalenie, która jest prawdziwa, a którą należy odrzucić.

Każdej z powyższych hipotez można przypisać dwie różne decyzje:

decyzja d

−

"obiekt jest"

decyzja d

−

"obiektu nie ma".

Zakłada się, że po zakończeniu procesu wykrywania nie może mieć miejsca

decyzja "nie wiem".

Jeśli spełniona jest hipoteza H

, to znaczy gdy obiekt faktycznie znajduje

się w przestrzeni i zostanie podjęta decyzja d

−

"obiekt jest", to mówi się o

poprawnym wykryciu (detekcji). Jeżeli spełniona jest powyższa hipoteza, a

zostanie podjęta decyzja d

−

"obiektu nie ma", to mówi się o tzw. fałszywym

spokoju (ang. false dismissal)

−

obiekt w powietrzu pozostał nie zauważony.

Taki błąd zwany jest często błędem II rodzaju i jest ze wszech miar

niepożądany. Jakość detekcji przy hipotezie H

można określić poprzez

odpowiednie prawdopodobieństwa warunkowe:

prawdopodobieństwo poprawnego wykrycia

)

(

(1.1)

oraz prawdopodobieństwo przepuszczenia celu

)

(

(1.2)

Ponieważ decyzje d

i d

stanowią wzajemnie wykluczające się zdarzenia

losowe i związane są z tą samą hipotezą H

(obiekt jest), suma

prawdopodobieństw detekcji i przepuszczenia celu równa jest jedności

(1.3)

Znajomość warunkowego prawdopodobieństwa detekcji pozwala zawsze

określić warunkowe prawdopodobieństwo fałszywego spokoju. Jeśli np. dla

jednego cyklu wykrywania prawdopodobieństwo detekcji wynosi 0,9, to

prawdopodobieństwo przepuszczenia celu wynosi 0,1. Oznacza to, że radar

gwarantuje wykrycie obiektu średnio w 90% przypadków, natomiast w 10

przypadkach na 100 obiekt pozostanie nie zauważony.

Jeśli spełniona jest hipoteza H

, to znaczy gdy obiektu nie ma w

obserwowanej przestrzeni, poprawną jest decyzja d

−

"obiektu nie ma". Jeżeli z

powodu zakłóceń sygnał odebrany zostanie źle zinterpretowany i podejmie się

decyzję d

−

"obiekt jest", mówi się o tzw. fałszywym alarmie (ang. false alarm).

Fałszywy alarm określany jest często jako błąd I rodzaju i także jest bardzo

niepożądany. Nie jest tak groźny jak np. niewykrycie wrogiego samolotu, ale

błędne informacje mogą obciążać zbędnie system obróbki danych. Może to

doprowadzić do zakłóceń w przekazywaniu i analizie danych prawidłowych.

Jakość wykrywania obiektu przy jego nieobecności w przestrzeni (spełnieniu

hipotezy H

) także określa się przez prawdopodobieństwa warunkowe:

prawdopodobieństwo fałszywego alarmu

)

(

(1.4)

oraz prawdopodobieństwo poprawnego niewykrycia

)

(

(1.5)

Ponieważ decyzje d

i d

wzajemnie się wykluczają i odpowiadają tej samej

hipotezie H

(obiektu nie ma), zachodzi podobnie jak dla prawdopodobieństw

detekcji i przepuszczenia celu:

(1.6)

Znajomość prawdopodobieństwa fałszywego alarmu F pozwala zawsze

określić prawdopodobieństwo poprawnego niewykrycia.

Jeżeli radar obserwuje pewien fragment przestrzeni i np. F=10

-6

, to

oznacza iż na 10

sondowań tego fragmentu pojawi się średnio jeden fałszywy

alarm, a 999999 przypadkach poprawnie da informację, że ten obszar

przestrzeni jest pusty. Stacja radarowa przeszukuje z reguły w jednostce czasu

cały obszar powietrzny wokół siebie

−

dużą liczbę nie pokrywających się małych

fragmentów przestrzeni. Jeśli mamy n fragmentów i dla pojedynczego

to dla całej przestrzeni prawdopodobieństwo fałszywego alarmu F

rośnie

proporcjonalnie do n. Prawdopodobieństwo poprawnego niewykrycia

)

(

stąd

−

≈

−

)

(

)

(

(1.7)

Ten powyższy fakt powoduje, że w teorii detekcji sygnałów radarowych

przyjmuje się typowo bardzo małe dopuszczalne prawdopodobieństwo

fałszywego alarmu dla pojedynczego fragmentu przestrzeni.

ndop

dop

Natomiast wartość prawdopodobieństwa detekcji D powinna być jak

największa. Są to warunki trudne do zrealizowania, zwłaszcza gdy obiekt

znajduje się w dużej odległości od radaru i energia sygnału odbitego od celu

jest niska. Granice strefy wykrywania radaru określa się odległością, na której w

jednym cyklu obserwacji prawdopodobieństwo detekcji D nie jest mniejsze od

założonej wartości dopuszczalnej D

dop

Tak więc głównymi parametrami określającymi jakość detektora w

systemach radarowych są prawdopodobieństwa warunkowe detekcji D oraz

fałszywego alarmu F. W granicach strefy wykrywania musi być spełnione:

F<F

dop

oraz D

≥

dop

1.2. POJĘCIE RYZYKA. OPTYMALNE KRYTERIA DECYZYJNE

Wymagania stawiane radarom tzn. jak najmniejsze F oraz jak największe

D są często sprzeczne. Aby uniknąć przepuszczenia celu, a w konsekwencji

tego np. ataku, wskazane jest podejmowanie decyzji "obiekt jest" w przypadku,

gdy odebrany sygnał jest mocno zakłócony i zniekształcony i trudno z

pewnością stwierdzić, że gdzieś w przestrzeni jest samolot. Takie postępowanie

prowadzi do wzrostu częstości występowania fałszywego alarmu. Zachodzi

więc konieczność znalezienia kompromisu między tymi sprzecznymi

wymaganiami. Trzeba znaleźć taką metodę obróbki danych radiolokacyjnych,

aby była ona optymalna z punktu widzenia całego procesu analizy. Uzyskany

tak sposób podejmowania decyzji o obecności lub nieobecności celu może

okazać się w pewnych konkretnych sytuacjach nie najlepszy. Powinien być

jednak optymalny w sensie statystycznym tzn. uwzględniać rozkłady

prawdopodobieństwa możliwych sytuacji w procesie wykrywania.

Przy poszukiwaniu optymalnych rozwiązań dla tego rodzaju problemów

zasadne jest posługiwanie się pojęciem ryzyka [10]. Można się nim z

powodzeniem posługiwać w innych dziedzinach nie związanych z systemami

radarowymi np. w telekomunikacji czy w gospodarce

−

wszędzie tam, gdzie

trzeba podjąć jedną z wielu możliwych do wyboru decyzji. W radiolokacji

pozwala ono na jednolite i dostatecznie ogólne potraktowanie problemu detekcji

sygnałów radiolokacyjnych oraz estymacji ich parametrów.

Do zbioru możliwych sytuacji podczas wykrywania obiektów zaliczane są

wszystkie możliwe stany wynikające w procesie radiolokacji oraz

przyporządkowane im decyzje. Mogą wystąpić a priori cztery poniższe sytuacje:

−

poprawne niewykrycie obiektu

−

fałszywy alarm

−

fałszywy spokój

−

poprawne wykrycie obiektu.

Analogicznie wygląda to w przypadku estymacji jakiejś wielkości zmieniającej

się w sposób ciągły. Przykładowo każdy pomiar odległości obiektu od radaru

polega na tym, że estymowanemu parametrowi

przypisuje się różniącą się

od niego estymatę (inaczej podejmuje się decyzję odnośnie wartości

parametru)

. W wyniku tego powstaje błąd estymacji

. Zbiór

możliwych sytuacji składa się ze wszystkich możliwych par

Każdej z możliwych sytuacji odpowiada określone prawdopodobieństwo.

W przypadku detekcji można mówić o prawdopodobieństwach P

i-tej sytuacji

(i=1,2,3,4), których suma równa jest jedności:

)

(

)

(

)

(

)

(

(1.8)

W przypadku estymacji parametrów sytuacje (

) opisywane są

dwuwymiarową gęstością prawdopodobieństwa p(

). Analogicznie, jak dla

detekcji suma prawdopodobieństw wszystkich możliwych sytuacji równa jest

jedności:

∫ ∫

∫

∞

−

∞

−

)

(

)

(

)

(

Każdej możliwej sytuacji można także przyporządkować pewną liczbę

mającą sens kosztu (straty) związanego z popełnionym błędem. Wielkość straty

zależy od ważności popełnionego błędu tzn. większemu błędowi

przyporządkowywana jest z reguły większa strata. Zwykle przyjmuje się, że

koszt bezbłędnych decyzji jest równy zeru. Przyporządkowywanie

poszczególnym błędom odpowiednich strat stanowi istotny etap w ogólnej

analizie danego problemu. Dlatego powinno się uwzględnić w nim wszelkie

dane o tym, na ile niepożądany jest jakiś rodzaj błędu. Strata może być

określona np. przez przesłanki typu ekonomicznego, czy strategicznego.

Trzeba by się np. zastanowić, co mogłoby być gorsze i więcej kosztować:

fałszywy alarm i poderwanie samolotów w celu przechwycenia widma, czy

przepuszczenie wrogiego samolotu i narażenie na atak z powietrza. Metody

znajdowania straty nie należą do teorii związanej z wykrywaniem. W wielu

przypadkach postać optymalnej reguły obróbki sygnałów (w przypadku detekcji

tzw. reguły decyzyjnej) nie jest zależna od konkretnej postaci funkcji strat.

Strata jest określona dla każdej możliwej sytuacji w sposób niezależny od

reguły decyzyjnej. Te dwa pojęcia można jednak ze sobą powiązać i otrzymać

w ten sposób ocenę jakości reguły decyzyjnej. W przypadku, gdy stany

(hipotezy dotyczące obiektu w przestrzeni) i związane z nimi decyzje są

wielkościami przypadkowymi o ustalonych prawdopodobieństwach, jako miernik

jakości reguły decyzyjnej można przyjąć stratę uśrednioną na wszystkie stany i

wszystkie decyzje. Wprowadza się pojęcie średniej (w sensie statystycznym)

straty zwanej ryzykiem. Dla zmiennej losowej skokowej

−

którą jest zbiór

wszystkich możliwych sytuacji w procesie detekcji ryzyko liczy się następująco:

∑

(1.9)

gdzie:

−

strata dla i-tej sytuacji

−

prawdopodobieństwo występowania i-tej sytuacji.

Dla zmiennej losowej ciągłej (estymacja parametrów) ryzyko liczy się z

poniższego wzoru:

∫

⋅

)

(

)

(

)

(

Przyglądając się procesowi detekcji widać, że błędnymi decyzjami są

sytuacje: gdy podejmie się decyzję o obecności celu, jeśli w rzeczywistości go

nie ma oraz gdy pozostawimy cel niezauważonym. Straty dla poprawnych

decyzji są równe zeru, tak więc wystarczy tylko wyznaczyć straty dla dwóch

wspomnianych wyżej błędów:

)

(

−

strata dla błędu fałszywego alarmu

)

(

−

oraz stratę dla błędu przepuszczenia celu.

Korzystając ze wzoru na prawdopodobieństwo warunkowe można zapisać, że:

⋅

)

(

)

(

)

(

)

(

oraz

⋅

)

(

)

(

)

(

)

(

gdzie P(H

) oraz P(H

) są prawdopodobieństwami bezwarunkowymi a priori

stanów H

i H

tzn. nieobecności lub obecności obiektu w przestrzeni. Jeśli (jak

wspomniano wyżej) koszty poprawnych decyzji są zerowe, to ryzyko w procesie

wykrywania ma postać:

)

(

)

(

⋅

(1.10)

Porównując różne systemy obrabiające dane radarowe za lepsze uznaje się te,

dla których ryzyko jest mniejsze. Pozwala to zdefiniować kryterium optymalizacji

dla detektorów nazywanym kryterium minimum ryzyka, znane także pod nazwą

kryterium Bayesa [4]. Projektując optymalny system wykrywania należy dążyć

do minimalizacji wyrażenia (1.10).

W przypadku estymacji parametrów strata związana z błędem jest funkcją

dwóch zmiennych r(

). Jeżeli np. przy estymacji odległości założyć, że

koszt zależy tylko od różnicy faktycznej odległości i jej estymaty

, to

otrzymuje się funkcję jednej zmiennej r(

). Jeżeli strata równa jest kwadratowi

błędu tzn.

)

(

)

(

−

śrkw

ε =

czyli ryzyko równe jest błędowi średniokwadratowemu, a minimum ryzyka

minimum błędu średniokwadratowego.

Kryterium minimum ryzyka jest kryterium bardzo ogólnym i łatwo jest (np.

zmieniając straty) przechodzić od niego do bardziej szczegółowych i prostszych

kryteriów. Jeśli we wzorze (1.10) wstawi się

to otrzyma się wyrażenie na ryzyko, które równe jest sumarycznemu

prawdopodobieństwu błędów detekcji:

)

(

)

(

)

(

)

(

⋅

(1.11)

Kryterium minimalizujące tak zdefiniowane ryzyko nosi nazwę kryterium

idealnego obserwatora. Jest bardziej szczegółowe od wspomnianego wcześniej

(ogólniejszego) kryterium minimum ryzyka, gdyż nie pozwala uwzględnić różnic

w kosztach popełnienia błędów fałszywego alarmu i przepuszczenia obiektu.

Zastąpiwszy D przez 1-D w ogólnym wyrażeniu na ryzyko można je

zapisać jak poniżej:

[

]

)

(

)

(

⋅

−

⋅

(1.12)

gdzie

)

(

)

(

⋅

Ponieważ

)

P(H

, kryterium minimalizacji ryzyka sprowadza się do

maksymalizacji wyrażenia

−

(1.13)

Kryterium to nosi nazwę kryterium wagowego. Sprowadza się do

maksymalizacji prawdopodobieństwa detekcji D i jednoczesnej minimalizacji

prawdopodobieństwa fałszywego alarmu F. Współczynnik

pełni rolę

współczynnika wagowego. Na jego wartość wpływają zależności między

kosztami poszczególnych błędów oraz prawdopodobieństw obecności lub

nieobecności celu w pewnym fragmencie przestrzeni.

Porównując dwa detektory, z których jeden jest optymalny w sensie

kryterium wagowego można zapisać, że:

opt

−

≥

−

Stąd:

)

(

opt

−

⋅

≥

to znaczy, że dla F=F

opt

opt

≥

Gdy F

≤

opt

, można zapisać inaczej, że

opt

≤

Oznacza to, że optymalny detektor daje najmniejsze prawdopodobieństwo

przepuszczenia celu (czyli największe prawdopodobieństwo detekcji) ze

wszystkich detektorów, dla których prawdopodobieństwo fałszywego alarmu

jest nie większe niż dla urządzenia optymalnego.

W wielu sytuacjach nie tylko prawdopodobieństwo a priori, ale i koszty

błędnych decyzji są trudne do podania, czy nawet do zdefiniowania. Występuje

to np. w przypadku detekcji sygnałów radarowych, gdzie trudno określić koszty

nie wykrycia celu i gdzie prawdopodobieństwa a priori P(H

) i P(H

) mogą nie

mieć sensu. W przypadku, gdy hipoteza H

występuje skrajnie rzadko, głównym

czynnikiem w ogólnych średnich kosztach jest częstość F prób, w których jest

mylnie podjęta decyzja d

, co oznacza popełnienie błędu I rodzaju. W związku z

tym pewna kosztowna czynność wykonywana jest niepotrzebnie. Na przykład w

systemie obrony przeciwlotniczej połączonej z radarem fałszywy alarm może

doprowadzić do wystrzelenia drogiego pocisku atakującego nie istniejący cel. W

takiej sytuacji analizę (projektowanie detektora) rozpoczyna się od określenia

wartości prawdopodobieństwa F, na którą obserwator może sobie pozwolić, a

następnie poszukuje się takiego sposobu decydowania, aby osiągnąć tą

wartość i jednocześnie uzyskać minimum możliwego prawdopodobieństwa D

określającego błąd II rodzaju. Taki sposób podejścia do problemu określa

kryterium Neymana-Pearsona [1, 4]. W systemie radarowym odpowiada to

maksymalizacji prawdopodobieństwa detekcji dla założonego poziomu

prawdopodobieństwa fałszywego alarmu.

1.3. PRZYKŁAD OPTYMALNEJ DETEKCJI

Zasadność stosowania kryterium wagowego (ogólnie kryterium minimum

ryzyka) w procesie detekcji można przedstawić na prostym przykładzie miernika

wskazówkowego, mierzącego wartość napięcia. W swej istocie jest to zbliżone

do procesu wykrywania sygnałów radiolokacyjnych. Wskazania miernika można

określić liczbą y, która jest albo sumą napięcia x i zakłócenia n

albo tylko napięciem zakłócającym

Wielkości n, x i y nie zmieniają się podczas pomiaru. Wartość sygnału x jest

dokładnie znana. Znany jest też rozkład zmiennej losowej n. Na podstawie

zmierzonego napięcia y, należy podjąć decyzję d

lub d

odnośnie obecności

lub nieobecności sygnału x, która powinna być optymalna w sensie kryterium

minimum ryzyka (lub równoważnego mu kryterium wagowego).

Dla potrzeb dalszych rozważań zostanie przyjęty gaussowski rozkład

zakłóceń n. Na wykresie 1.1 na następnej stronie zostały pokazane warunkowe

gęstości prawdopodobieństw zmiennej losowej y przy nieobecności sygnału x

(warunek H

) i przy jego obecności (warunek H

)

(

)

(

(1.14)

)

(

)

(

(1.15)

Ponieważ rozkład zmiennej losowej n jest normalny, a x jest stałą wartością, to

krzywa p

x+n

(y) jest przesunięta względem krzywej p

(y) o wartość x.

)

(

)

(

−

(1.16)

Funkcja (reguła) decyzyjna rozwiązująca problem detekcji sygnału x będzie

zależeć od y i w zależności od niego przyjmować wartości 0 lub 1. Przykładowa

reguła decyzyjna d(y) (niekoniecznie optymalna) została pokazana na wykresie

1.1. W danym przypadku jeśli y

<y<y

, podejmowana jest decyzja "sygnał jest".

Dla pozostałych wartości y podejmowana jest decyzja o nieobecności sygnału

x. Prawdopodobieństwami D i F określa się prawdopodobieństwa, że sygnał y

trafi do przedziału (y

)

−

odpowiednio, gdy y=x+n (sygnał i zakłócenie) oraz

gdy y=n (samo zakłócenie). Dla reguły decyzyjnej z wykresu liczone są

następująco:

∫

)

(

(1.17)

∫

)

(

(1.18)

Gdyby znane były prawdopodobieństwa a priori hipotez H

i H

tzn. P(H

) i P(H

) można by

określić ogólną funkcję gęstości prawdopodobieństwa wyników y:

)

(

)

(

)

(

)

(

)

(

⋅

Rys. 1.1. Krzywe warunkowych rozkładów prawdopodobieństw zmiennej losowej y dla
prawdziwości hipotezy H

−

(y) i H

−

x+n

(y) oraz przykładowa reguła decyzyjna d(y)

Prawdopodobieństwom tym odpowiadają zakreskowane pola na pierwszym z

wykresów. Wprowadzając dowolną regułę decyzyjną do wyrażeń na D i F

można je zapisać następująco:

∫

∞

−

⋅

)

(

)

(

(1.19)

∫

∞

−

⋅

)

(

)

(

(1.20)

Na odcinkach, w których d(y)=0, obie całki będą równe zero. Natomiast tam,

gdzie d(y)=1, całki będą równe zakreskowanym na wykresie polom pod

krzywymi odpowiednio p

x+n

(y) i p

(y). Dwie powyższe zależności są prawdziwe

dla dowolnej postaci reguły decyzyjnej d(y).

Biorąc pod uwagę powyższe zależności, wyrażenie D-

F odpowiadające

kryterium wagowemu można przedstawić następująco:

[

]

∫

∞

−

⋅

−

)

(

)

(

)

(

(1.21)

gdzie

)

(

)

(

)

(

Zgodnie z kryterium wagowym optymalny detektor powinien

zagwarantować maksimum powyższej całki. Aby spełnić ten warunek, należy

dla każdego y maksymalizować wyrażenie podcałkowe. Można to uzyskać

przez odpowiednią modyfikację reguły decyzyjnej d(y). Funkcja d(y) przyjmuje

tylko wartości 0 lub 1, w związku z tym wyrażenie podcałkowe albo się zeruje,

albo jest mnożone przez 1. Ponieważ każda liczba dodatnia jest większa od

zera, a zero jest większe od dowolnej liczby ujemnej, łatwo jest wskazać

sposób maksymalizacji wyrażenia podcałkowego, a tym samym całego

wyrażenia: jeśli wyrażenie podcałkowe jest dodatnie przyjmuje się d(y)=1,

natomiast gdy jest mniejsze od zera zakłada się d(y)=0. Ponieważ gęstości

prawdopodobieństwa p

(y) i p

x+n

(y) są nieujemne, to optymalną regułę

decyzyjną dla procesu detekcji można zapisać jak niżej:







)

(

)

(

)

(

dla

opt

(1.22)

Wielkość

)

(

)

(

)

(

nosi nazwę stosunku wiarygodności [1, 3, 12]. Stosunek wiarygodności jest

liczbą nieujemną i wyraża się go ilorazem gęstości prawdopodobieństw tej

samej realizacji y przy dwóch różnych warunkach: kiedy występuje sygnał i

zakłócenie oraz gdy występuje jedynie zakłócenie. Jeśli dla konkretnej wartości

y wielkość

(y)>

, to znaczy, że sytuacja, gdy y=x+n jest bardziej

prawdopodobna od sytuacji, gdy y=n. Ogólnie można powiedzieć, że jeśli dane

są dwa rozkłady tej samej zmiennej losowej x, z tym że pierwszy zależy od

zbioru parametrów

, a drugi od zbioru

i iloraz wiarygodności

)

;

(

)

;

(

jest większy od założonego wcześniej progu

(w szczególnym przypadku

=1), to znaczy że zbiór parametrów

jest bardziej prawdopodobny niż zbiór

Z powyższych rozważań widać, że zamiast bardzo ogólnego kryterium

minimum ryzyka zasadniejsze jest stosowanie kryterium ilorazu wiarygodności.

Zgodnie z nim decyzję o obecności obiektu w określonym fragmencie

przestrzeni podejmuje się, gdy iloraz wiarygodności przekracza założoną

wartość progową

. Tak sformułowane kryterium optymalnego wykrywania

dobrze się nadaje do zastosowań praktycznych w systemach radarowych.

Ważne jest to, że można je stosować niezależnie od postaci rozkładu

prawdopodobieństwa zakłóceń.

W przykładzie miernika wskazówkowego przyjęto normalny rozkład

zakłóceń

)

;

(

, o zerowej wartości średniej i wariancji

. W sytuacji, gdy nie

ma sygnału użytecznego tzn. y=n, to:

)

(

−

W sytuacji gdy y=x+n, ponieważ p

x+n

(y)=p

(y-x), można zapisać:

)

(

)

(

−

Dla takich rozkładów iloraz wiarygodności ma postać:

)

(

)

(

)

(

−

⋅

Przebieg krzywej

(y) dla x>0 pokazany jest na rysunku 1.2 na następnej

stronie (na osi rzędnych zaznaczono wartość progową

). Dzięki ścisłej

monotoniczności funkcji

(y), warunek

(y)>

można sprowadzić do warunku

y>y

, a

(y)<

odpowiednio do y<y

. Wcześniejsze wyrażenie na optymalną

regułę decyzyjną można zapisać jak niżej:







)

(

dla

opt

(1.23)

Rys. 1.2. Krzywa ilorazu wiarygodności

)

w zależności od wartości zmiennej losowej y

Porównując je z zaproponowaną regułą decyzyjną widać, że tamta nie była

optymalna. Dalej na wykresie 1.3 pokazano przebiegi rozkładów

prawdopodobieństw p

(y) i p

x+n

(y) wraz z optymalną regułą decyzyjną.

Prawdopodobieństwo fałszywego alarmu F dla tej reguły odpowiada polu pod

krzywą p

(y) dla y

≥

. Wielkość y

nazywana jest progiem detekcji. Przy

założonym poziomie zakłóceń prawdopodobieństwo fałszywego alarmu F

zależy wyłącznie od y







 Φ

−













−

∫

−

∞

−

∞

)

(

)

(

gdzie

jest tzw. funkcją błędu. Oznacza to, że wartość progu y

można ustalić

bezpośrednio z zakładanego prawdopodobieństwa fałszywego alarmu.

Podejście takie odpowiada wspomnianemu wcześniej kryterium Neymana-

Pearsona [1, 4]. Jest wygodne w zastosowaniach praktycznych, ponieważ

pozwala uniknąć konieczności uwzględnienia prawdopodobieństw a priori

dotyczących obecności lub nieobecności sygnału. W systemach radarowych

odpowiada to prawdopodobieństwom stanów H

i H

−

odpowiednio P(H

) i

P(H

) czyli obecności i nieobecności obiektu w określonym fragmencie

przestrzeni powietrznej.

Prawdopodobieństwo detekcji D sygnału x odpowiada polu pod krzywą

x+n

(y) dla y

≥

. Przy założonym poziomie zakłóceń

, zależy nie tylko od progu

detekcji y

, ale i od oczekiwanej wartości sygnału x, co opisywane jest

zależnością:









−













−

∫

−

∞

−

∞

)

(

)

(

Rys. 1.3. Krzywe warunkowych rozkładów prawdopodobieństw zmiennej losowej y dla
prawdziwości hipotezy H

−

(y) i H

−

x+n

(y) oraz optymalna reguła decyzyjna d

opt

(y)

Na wykresie 1.4 na kolejnej stronie pokazano krzywe zwane charakterystykami

detekcji. Jak widać, im większa jest wartość progu y

, czyli mniejsze

prawdopodobieństwo fałszywego alarmu F, tym bardziej w prawo przesuwa się

wykres D(x). Oznacza to, że aby zagwarantować takie samo

prawdopodobieństwo detekcji D potrzebny jest wyższy poziom sygnału

użytecznego x.

Rys. 1.4. Charakterystyki detekcji dla przypadku wskaźnika wychyłowego

1.4. WERYFIKACJA HIPOTEZ ZA POMOCĄ

POMIARÓW WIELOKROTNYCH

Teoria i przykład przedstawiony we wcześniejszych punktach mogą być

łatwo rozszerzone na przypadki, w których wybór jednej z hipotez H

lub H

dokonywany jest na podstawie większej liczby pomiarów. Mierzone wielkości y

, ... y

N-1

mogą reprezentować kolejne pomiary tej samej wielkości fizycznej,

pomiar N różnych parametrów, albo dowolną kombinację tych dwóch

możliwości. Wielkości te mogą być opisane łącznymi, wielowymiarowymi

funkcjami gęstości prawdopodobieństwa p(y/H

)=p

(y)=p

, y

, ..., y

N-1

) i

p(y/H

)=p

(y)=p

, y

, ... y

N-1

) odpowiednio dla hipotez H

i H

. Zadanie polega

na wyborze jednej z tych dwóch funkcji gęstości prawdopodobieństwa, która

lepiej opisuje wynik pomiaru. Słowo „lepiej” należy rozumieć w sensie

odpowiedniego kryterium decyzyjnego wprowadzonego w poprzednich

punktach. Jeśli funkcje p

(y) i p

(y) nie zawierają nieznanych parametrów, to

odpowiadające im hipotezy można nazwać „prostymi”. Obie funkcje gęstości

pomnożone przez dy

...dy

N-1

określają prawdopodobieństwo zdarzenia

polegającego na tym, że zmienna losowa y

przyjmie wartość z przedziału

〈

+dy

〉

, y

z przedziału

〈

, y

+dy

〉

itd. Oczywiście dla obu funkcji gęstości

musi zachodzić związek ,że

...

)

,...,

(

...

⋅

∫∫

∞

−

W procesie detekcji każdej realizacji (y

, y

, ..., y

N-1

) należy

przyporządkować konkretną decyzję. Oznacza to, że na zbiorze realizacji może

być określona reguła decyzyjna d(y)=d(y

,...y

N-1

), przyjmująca w zależności

od swych argumentów dwie wartości: 0 lub 1, którym odpowiadają decyzje d

oraz d

o nieobecności lub obecności sygnału użytecznego.

Inaczej mówiąc każdemu punktowi wielowymiarowej przestrzeni należy

przypisać odpowiednią wartość: 0 lub 1. Regułę decyzyjną można wtedy

wyrazić jako podział wspomnianej przestrzeni na dwa obszary: R

oraz R

. Jeśli

wynik pomiaru (punkt) znajduje się w obszarze R

, to podejmowana jest

decyzja d

o słuszności hipotezy H

, zaś gdy punkt leży w obszarze R

wybierana jest hipoteza H

. Obszary R

i R

są rozdzielone powierzchnią d,

zwaną „powierzchnią decyzyjną”. Regułę decyzyjną można (jak w poprzednim

punkcie dla pojedynczego pomiaru

−

zależność 1.23) opisać zbiorem

nierówności dla N zmiennych y

≤

N-1).

Aby łatwiej zilustrować sytuację wygodnie jest przedstawić zbiór N

wyników pomiaru geometrycznie, jako jeden punkt w N-wymiarowej przestrzeni

kartezjańskiej o współrzędnych y=(y

, y

, ..., y

N-1

). Na rysunku 1.5 pokazano

przykładowe obszary decyzji i powierzchnię decyzyjną dla N=2. Położenie i

kształt powierzchni decyzyjnej d zależy od wybranego kryterium

−

inne będzie

dla ogólnego kryterium minimum ryzyka, a inne dla kryterium Neymana-

Pearsona.

Rys. 1.5. Obszary decyzyjne dla dwóch pomiarów (N=2) oddzielone powierzchnią decyzyjną d

Jeśli znane są koszty błędnych decyzji (dla poprawnych można przyjąć

zerowy koszt) tzn. straty dla fałszywego alarmu oraz przeoczenia obiektu w

przestrzeni, oraz bezwarunkowe prawdopodobieństwa a priori P(H

) i P(H

)

odpowiednio dla hipotez H

i H

, można zastosować ogólne kryterium Bayesa

polegające na minimalizacji ryzyka.

W sytuacji, gdy nie jest możliwa ocena prawdopodobieństw a priori P(H

) i

P(H

), a także nie są dokładnie znane koszty towarzyszące różnym decyzjom,

dogodne jest zastosowanie kryterium Neymana-Pearsona. Wymaga ono

takiego położenia powierzchni decyzyjnej d, aby zminimalizować

prawdopodobieństwo przepuszczenia obiektu D (czyli zmaksymalizować

prawdopodobieństwo detekcji D=1-D ), wtedy gdy narzucone jest dopuszczalne

prawdopodobieństwo fałszywego alarmu F. Obserwator w jakiś sposób

decyduje o odpowiednim doborze prawdopodobieństwa F. Zwiększanie go

oznacza zwiększenie prawdopodobieństwa D, czemu odpowiada dokładniejsze

stwierdzenie występowania hipotezy H

. Posługując się analogią do

jednowymiarowego przykładu detekcji łatwo można znaleźć optymalną

powierzchnię d i odpowiadającą jej regułę decyzyjną d(y)=d(y

, y

, ..., y

N-1

Prawdopodobieństwa fałszywego alarmu oraz detekcji dla dowolnej funkcji

decyzyjnej można zapisać jak dalej (analogicznie do wzorów 1.19, 1.20):

∫∫

∞

−

⋅

...

)

,...

(

,...)

(

...

(1.24)

∫∫

∞

−

⋅

...

)

,...

(

,...)

(

...

(1.25)

Dla wyboru optymalnej reguły można się posłużyć jak poprzednio wyrażeniem

na kryterium wagowe: D-

F. Można je przedstawić w poniższej postaci:

[

]

∫∫

∞

−

⋅

−

⋅

−

...

,...)

(

)

,...

(

,...)

(

...

(1.26)

gdzie

)

,...

(

)

,...

(

,...)

(

jest stosunkiem wiarygodności. Ponieważ optymalnym z punktu widzenia

kryterium wagowego jest rozwiązanie dające maksimum powyższej całki,

otrzymuje się je tak samo jak dla wykrywania jednowymiarowego w poprzednim

punkcie (wzór. 1.22):







,...)

(

,...)

(

,...)

(

dla

opt

(1.27)

Znając postać stosunku wiarygodności można przekształcić zapis powyższej

reguły decyzyjnej uzależniając ją od zmiennych (y

,...) i wartości progowej Y

uzyskując dzięki temu wyrażenie na powierzchnię decyzyjną d.

Tak więc optymalna powierzchnia decyzyjna d, dzieląca przestrzeń

realizacji na obszary R

i R

jest jedną z rodziny powierzchni, określonych przez

parametr

. Wybrana jest taka wartość parametru, dla której

prawdopodobieństwo fałszywego alarmu F równe jest wartości narzuconej.

Kryterium Neymana-Pearsona prowadzi do procedury, w której stosunek

wiarygodności porównywany jest z ustaloną wartością

; wybiera się hipotezę

jeżeli

(x)<

, a hipotezę H

, jeśli

(x)>

1.5. PRZYKŁAD OPTYMALNEJ DETEKCJI NA

PODSTAWIE WIELU POMIARÓW

Teorię przedstawioną w punktach 1.1

−

1.3 w połączeniu z rozważaniami z

punktu ostatniego 1.4 można zastosować rozważając przykładowy problem

detekcji sygnału x(t), którego wszystkie parametry (jak amplituda, faza, kształt

oraz czas pojawienia się) są znane.

Sygnał x(t) może występować w obecności białego szumu gaussowskiego

n(t) o wartości oczekiwanej równej zero. Obserwator odbiera sygnał wejściowy

y(t) będący np. echem radarowym odbitym od jakiegoś obiektu, obserwowanym

w czasie 0

T, zwanym przedziałem obserwacji.

Na podstawie przebiegu y(t) musi zostać podjęta decyzja odnośnie

wyboru jednej z dwóch możliwych hipotez prostych: H

−

że sygnał użyteczny

nie występuje i odebrany przebieg zawiera tylko szum tzn.

)

(

)

(

albo hipotezę H

−

że sygnał odebrany składa się z sumy oczekiwanego

sygnału i szumu:

)

(

)

(

)

(

Powyższe wyrażenia różnią się od tych występujących w przypadku

jednowymiarowym (wskaźnik wychyłowy) tym, że tutaj występują funkcje

czasu. Przyjęte zostało pewne kryterium, zapewniające sukces przy dużej

liczbie podejmowanych decyzji tego typu, będące miarą poprawności

działania (wybór kryterium ma wpływ na sposób oceny zbieranych

danych).

Przykładowym sygnałem użytecznym może być impuls prostokątny o

czasie trwania T

T, występującym w określonym (znanym) czasie wewnątrz

przedziału obserwacji:







∪

∈

dla

;

)

(

(1.28)

Po zakończeniu obserwacji (po czasie T) podejmowana jest decyzja o

obecności lub nieobecności sygnału. Istnienie szumu powoduje popełnianie

przypadkowych błędów i zadaniem obserwatora jest minimalizacja

prawdopodobieństwa ich popełnienia.

Rozważania i wyprowadzenie reguł decyzyjnych w poprzednich punktach

prowadzono przy założeniu skończonej liczby obserwacji. Jednak zgodnie z

twierdzeniem o próbkowaniu każdy sygnał o ograniczonym paśmie do jakiejś

częstotliwości f

max

można poddać dyskretyzacji, nie tracąc informacji w nim

zawartych, jeśli tylko częstotliwość próbkowania nie będzie mniejsza od 2f

max

[13]. Każdej chwili próbkowania t

≤

N-1) można przypisać trójkę liczb

określającą wartości sygnałów: x

=x(t

), n

=n(t

) oraz y

=y(t

) odpowiednio dla

sygnału użytecznego, zakłócenia i ich sumy. Istnieje zatem możliwość opisania

przypadkowych funkcji czasu za pomocą równoważnych im wielowymiarowych

zmiennych losowych. Obserwator wybiera jedną z dwóch hipotez na podstawie

wielu wartości y

. Wartości te charakteryzowane są przez łączne N-wymiarowe

funkcje gęstości prawdopodobieństwa:

)

,...,

(

−

oraz

)

,...,

(

−

pod warunkiem, że prawdziwa jest odpowiednio hipoteza H

albo H

Podejmując decyzję można się opierać na wartości stosunku wiarygodności

)

,...,

(

)

,...,

(

)

,...,

(

−

który porównywany jest z wartością progową

. Jeśli

(y)<

, podejmowana

jest decyzja d

, że słuszna jest hipoteza H

(odebrano tylko zakłócenia), jeśli

(y)

podejmowana jest decyzja d

tzn. że odebrano sygnał użyteczny.

Wielkość

zależy od wybranego kryterium decyzyjnego.

W przypadku, gdy n(t) jest gaussowskim szumem białym, to zmienne

losowe y

są statystycznie niezależne i jeśli prawdziwa jest hipoteza H

, to

łączna funkcja gęstości prawdopodobieństwa wynosi [3, 12]:

)

exp(

)

(

)

,...,

(

∑

−

⋅

πσ

(1.29)

Gdy istnieje sygnał użyteczny (prawdziwa hipoteza H

), to łączna funkcja

gęstości prawdopodobieństwa jest dana przez powyższe wyrażenie, z tym że y

zastąpiono przez (y

-x

Inaczej mówiąc wielkości y

są niezależnymi zmiennymi losowymi o wartości

oczekiwanej x

)

(

exp(

)

(

)

,...,

(

∑

−

⋅

πσ

(1.30)

Stosunek wiarygodności ma postać:

)

exp(

)

,...,

(

∑

−

(1.31)

Jeśli

(y)<

podejmowana jest decyzja d

. Logarytmując i przekształcając tą

nierówność dostaje się:

∑

−

(1.32)

Tak więc decyzję można podejmować opierając się na wartości poniższej

statystyki

∑

−

(1.33)

i porównywaniu jej z ustaloną wartością G

W N-wymiarowej przestrzeni realizacji zmiennych losowych y

powierzchnia decyzyjna d jest hiperpłaszczyzną

const

∑

−

do której jest prostopadły wektor o współrzędnych x

Zapisując wyrażenie na statystykę G

dla sygnałów z ciągłym czasem

dochodzi się do poniższej całki:

∫

)

(

)

(

zwanej całką korelacyjną. O systemie decyzyjnym, w którym statystyka G

lub

G jest porównywana z pewnym progiem G

, można powiedzieć, że wzajemnie

koreluje sygnał wejściowy y(t) z sygnałem oczekiwanym x(t). System taki

nazywa się detektorem korelacyjnym. Jeżeli G

, przyjmuje się, że

prawdziwa jest hipoteza H

−

"sygnału nie ma", jeśli G

wybierana jest

hipoteza H

−

"sygnał jest". Wielkość G

jest zależna od stosowanego kryterium

decyzyjnego. Przy zastosowaniu kryterium Neymana-Pearsona wartość

prawdopodobieństwa fałszywego alarmu F jest z góry ustalona, zwykle na

podstawie względnej liczby błędów I rodzaju, jaką obserwator może dopuścić.

Statystyka G

jest zmienną losową opisaną przez jedną z dwóch funkcji

gęstości prawdopodobieństwa: p

) i p

) odpowiednio gdy sygnał

użyteczny jest zawarty w sygnale y(t) oraz gdy go nie ma. Ponieważ G

jest

kombinacją liniową zmiennych losowej y(t) mających rozkład Gaussa, także jest

zmienną losową gaussowską. Jeśli prawdziwa jest hipoteza H

, to wartość

oczekiwana (średnia) statystyki G

równa jest

bo E[n(t)]=0. W przypadku prawdziwości hipotezy H

, wartość średnia sygnału

wejściowego y(t) w każdej chwili wynosi x(t), gdyż sygnał i szum się dodają.

Wartość oczekiwana EG

(łatwo ją wyliczyć korzystając z własności wartości

oczekiwanej) [12] wynosi wtedy

∑

−

(1.34)

i jest równa energii sygnału x(t). Co do wariancji, to są one takie same w

przypadku obu hipotez i wynoszą:

⋅

∑

−

(1.35)

Funkcje gęstości prawdopodobieństwa statystyki G

mogą być zapisane

według poniższych równań:

)

exp(

)

(

−

)

(

exp(

)

(

−

Jeżeli G

i sygnału nie ma popełnia się błąd I rodzaju (podnosi się

fałszywy alarm). Prawdopodobieństwo fałszywego alarmu wynosi

∫

∞

)

(

(1.36)

Jeśli podejmuje się decyzję d

tzn. stwierdza się istnienie sygnału, gdy

rzeczywiście występuje (G

, gdy hipoteza H

jest prawdziwa) oznacza to,

że wykryto sygnał x(t).

Prawdopodobieństwo detekcji (jego wykrycia) wynosi

∫

∞

)

(

(1.37)

Z drugiej strony, gdy G

, a sygnał występuje, popełniany jest błąd II rodzaju,

czyli przeocza się sygnał użyteczny. Prawdopodobieństwo błędu przeoczenia

wynosi 1-D.

2. ANALIZA UKŁADÓW OPTYMALNEJ DETEKCJI

RZECZYWISTYCH SYGNAŁÓW RADAROWYCH

2.1. DETEKCJA SYGNAŁU O NIEZNANYCH PARAMETRACH

W poprzednim rozdziale rozpatrywano optymalne metody detekcji i

problemy z nimi związane w przypadku, kiedy oczekiwany sygnał, który należy

wykryć jest całkowicie znany. Podejmowanie decyzji odnośnie dwóch hipotez

prostych H

i H

na podstawie ustalonej liczby pomiarów polegało na wyborze

jednej z dwóch funkcji gęstości prawdopodobieństwa p

(y) oraz p

(y), bardziej

odpowiadającej zaobserwowanym wartościom y=(y

, y

, ..., y

N-1

) N zmiennych

losowych. Wybór ten miał być optymalny ze względu na przyjęte kryterium

decyzyjne (Bayesa lub Neymana-Pearsona), zapewniające sukces

statystycznie w większości sytuacji. Zakładano przy tym, że w oba rozkłady

prawdopodobieństwa były całkowicie znane.

W rzeczywistych systemach takie idealne sytuacje należą do rzadkości.

Często zdarza się, że jedna lub obie funkcje gęstości prawdopodobieństwa

(y) i p

(y) zależą od pewnego zbioru M parametrów

, ...,

M-1

, których

wartości nie są znane obserwatorowi. Wtedy znalezienie optymalnych kryteriów

decyzyjnych staje się trudniejsze.

W niniejszym rozdziale rozpatrzona zostanie sytuacja, gdy rozkład p

(y)

odpowiadający hipotezie H

tzn. gdy y(t)=n(t) nie zależy od żadnych nieznanych

parametrów, natomiast p

(y) odpowiadający hipotezie H

tzn. y(t)=x(t)+n(t)

zależy od M parametrów, których zbiór można traktować jako wektor

ββββ

,...,

M-1

) w pewnej M-wymiarowej przestrzeni parametrów. Jeśli

składowe wektora

ββββ

są zmiennymi losowymi, to funkcja rozkładu p

(y/

ββββ

) jest

warunkową funkcją gęstości prawdopodobieństwa. Hipoteza H

zwana jest

wtedy prostą, a H

złożoną. Badaną regułę decyzyjną opisuje się jak

poprzednio tzn. przez podział N-wymiarowej przestrzeni kartezjańskiej wartości

próbek y

na obszary R

i R

. Wybiera się hipotezę H

, gdy punkt o

współrzędnych danych przez zmierzone wartości y=(y

,...,y

N-1

) leży na

obszarze R

, zaś H

−

gdy należy do R

. Powierzchnię d wybiera się tak, aby

reguła decyzyjna była w pewnym sensie optymalna.

Przykładowo w rzeczywistych systemach radarowych funkcja gęstości

prawdopodobieństwa p

(y) jest funkcją n próbek przebiegu wejściowego y(t)

zawierającego jedynie szum (zwykle biały o rozkładzie Gaussa), zaś funkcja

(y/

ββββ

) jest funkcją gęstości prawdopodobieństwa tych samych obserwacji

przebiegu wejściowego, gdy odbiera się także sygnał użyteczny x(t,

ββββ

), gdzie

składowe wektora

ββββ

reprezentują pewne parametry sygnału. Impulsowy sygnał

wąskopasmowy wysyłany przez radar da się zapisać jako:

)]

(

cos[

)

(

−

gdzie A jest znaną amplitudą,

−

pulsacją, a t

chwilą wysłania impulsu.

Jednak sygnał echa radarowego w zależności od tego, gdzie jest i jak się

zachowuje cel, ma już wiele nieznanych parametrów. Można go przedstawić

poniższym wyrażeniem:

)

(

)

(

)

(

)

(

)

(

−

⋅

−

⋅

gdzie:

=X(t)

−

amplituda sygnału;

−

opóźnienie wynikające z odległości celu od radaru;

−

pulsacja dopplerowska zależna od prędkości celu;

−

faza sygnału.

Obserwator może znać tylko niektóre z parametrów lub szerokie granice, w

jakich się one zmieniają. Można zatem powiedzieć, że próbuje wykryć jeden z

sygnałów należących do zbioru wykluczających się sygnałów o postaci jak

wyżej.

Najbardziej ogólnym, chociaż rzadko w praktyce spotykanym jest

przypadek, gdy wszystkie prawdopodobieństwa a priori: P(H

), P(H

), łączna

funkcja gęstości prawdopodobieństwa parametrów p(

ββββ

) oraz koszty są dobrze

określone. Można wtedy stosować kryterium Bayesa polegające na

minimalizacji ryzyka.

Jeśli P(H

), P(H

) i koszty decyzji nie są znane (sytuacja typowa w

systemach radarowych), stosuje się kryterium Neymana-Pearsona.

Prawdopodobieństwo błędu I rodzaju (fałszywego alarmu) wynosi

∫

−

...

)

,...,

(

gdzie R

jest obszarem przestrzeni próbek y

, w którym prawdziwa jest hipoteza

. Prawdopodobieństwo błędu II rodzaju (przepuszczenia obiektu) w

przypadku detekcji sygnałów z nieznanymi parametrami

ββββ

jest teraz funkcją

tych parametrów:

∫

−

...

)

,...,

(

)

(

podobnie R

jest obszarem, na którym wybierana jest hipoteza H

. Dla

zadanego zbioru wartości parametrów

ββββ

powierzchnia decyzyjna d

rozdzielająca obszary R

i R

powinna być taka, aby dla ustalonego poziomu F

prawdopodobieństwo

)

(

było minimalne.

Przy stosowaniu kryterium Neymana-Pearsona najwygodniejszy jest

przypadek, gdy usytuowanie powierzchni d minimalizujące

)

(

dla zadanego

F, nie zależy od wartości parametrów. Wtedy ta sama powierzchnia decyzyjna d

jest optymalna dla wszystkich wartości

ββββ

i dla funkcji gęstości

prawdopodobieństwa zmiennej

ββββ

mówi się, że strategia stanowi jednostajnie

najmocniejszy test hipotezy H

względem hipotezy H

. Zilustrować to można na

poniższym przykładzie: niech próbki y będą niezależnymi zmiennymi losowymi

o rozkładach Gaussa z wariancją

i niech wartość średnia x wynosi zero przy

hipotezie H

i x>0 przy hipotezie H

. Wtedy łączne funkcje gęstości

prawdopodobieństwa mają postać:

)

(

exp(

)

(

)

(

)

exp(

)

(

)

(

∑

−

⋅

−

⋅

πσ

O wartości średniej x wiadomo tylko tyle, że jest dodatnia. Dla ustalonej

wartości x wyliczywszy stosunek wiarygodności

(y), który porównywany byłby

z progiem

l , logarytmując następnie nierówność

(y)<

l stronami i

odpowiednio przekształcając uzyskuje się optymalną regułę decyzyjną, która

jest równoważna porównywaniu średniej

∑

−

z pewną ustaloną wartością graniczną Y

określoną przez poniższe wyrażenie:

Hipotezę H

wybiera się gdy Y<Y

, zaś H

gdy Y>Y

. Wartość graniczną ustala

się tak, aby uzyskać żądany poziom prawdopodobieństwa fałszywego alarmu

zgodnie ze wzorem, że

∫

∞

)

(

gdzie p

(Y) jest rozkładem Gaussa o zerowej średniej i wariancji

/N.

Powierzchnia decyzyjna d jest w tym przypadku płaszczyzną

∑

−

i jest niezależna od istniejącej przy hipotezie H

średniej x. Dlatego test jest

jednostajnie najmocniejszy, gdy występująca średnia x jest dodatnia i może być

stosowany nawet wtedy, gdy aktualna wartość x nie jest znana. Niestety raczej

wyjątkowo powierzchnia decyzyjna jest niezależna od parametrów

ββββ

. W wielu

przypadkach jedna i ta sama powierzchnia d(

ββββ

) nie jest optymalna dla

wszystkich możliwych wartości

ββββ

, a więc test jednostajnie najmocniejszy nie

istnieje. Tak jest na przykład, kiedy średnia x z powyższego przykładu może

być dodatnia lub ujemna. Należy wtedy skorzystać z danej funkcji gęstości

prawdopodobieństwa a priori p(

ββββ

) jeśli jest znana, lub wybrać pewną funkcję

gęstości prawdopodobieństwa p(

ββββ

) dobrze opisującą zjawisko i minimalizować

prawdopodobieństwo błędu II rodzaju uśrednione względem wszystkich

możliwych kombinacji parametrów:

śr

)

(

)

(

)

(

)

(

)

(

∫

⋅

Obserwator musi się zadowolić regułą, która jest poprawna tylko w odniesieniu

do zbioru możliwych wartości parametrów

ββββ

o rozkładzie zgodnym z przyjętym.

Może jednak istnieć jeden rozkład p

max

(

ββββ

), dla którego znaleziona wartość

)

(

śr

jest największa. Jest to najmniej korzystny rozkład przy kryterium

Neymana-Pearsona. W sytuacji, gdy wartości parametrów są nieznane

rozsądne jest stosowanie w kryterium Neymana-Pearsona najmniej

korzystnego rozkładu parametrów, ponieważ obserwator może być pewny, że

jego strategia daje założone przez niego prawdopodobieństwo fałszywego

alarmu F oraz prawdopodobieństwo przeoczenia obiektu nigdy nie przekroczy

wartości

)

(

)

(

)

(

max

śr

niezależnie od tego, który zbiór wartości parametrów rzeczywiście wystąpi w

następnych próbach [4].

2.2. STOSUNEK WIARYGODNOŚCI DLA SYGNAŁU

O NIEZNANYCH PARAMETRACH

W punkcie 1.5 wyprowadzono wyrażenie na stosunek wiarygodności dla

sygnału całkowicie znanego, pojawiającego się w obecności szumu białego o

rozkładzie Gaussa. W niniejszym punkcie wyprowadzone zostanie wyrażenie

na stosunek wiarygodności dla sygnału o nieznanych parametrach.

W sytuacji, gdy odbierany przebieg składa się z sygnału użytecznego o

nieznanych parametrach oraz zakłóceń, można go zapisać poniższą sumą:

)

(

)

,...,

(

)

(

−

Dyskretne wartości y

=y(t

) i zbiór parametrów

ββββ

można opisać łączną funkcją

gęstości prawdopodobieństwa. Jeśli założyć, że rozkład parametrów p(

ββββ

) jest

znany, to zgodnie z twierdzeniem Bayesa można zapisać jak niżej:

)

,...,

(

)

,...,

(

)

,...,

(

)

,...,

(

)

,...,

(

−

⋅

(2.1)

Całkowanie w nieskończonym przedziale bezwarunkowych lub warunkowych

gęstości prawdopodobieństwa daje zawsze jeden, tak więc:

∫ ∫

∞

−

∞

−

...

)

,...,

(

...

Całkując wyrażenie (2.1) dostaje się:

∫ ∫

∞

−

∞

−

⋅

...

)

,...,

(

)

,...,

(

...

)

,...,

(

(2.2)

Ponieważ p(y

,...,y

N-1

)=p

,...,y

N-1

), lewa strona wyrażenia (2.2) jest gęstością

prawdopodobieństwa realizacji wielkości y

, przy założeniu odebrania sygnału

wraz z zakłóceniami, mającymi rozkład a priori p(

ββββ

). Jeśli podzielić wyrażenie

(2.2) stronami przez gęstość prawdopodobieństwa tych samych wielkości y

przy założeniu odebrania jedynie zakłóceń p

,...,y

N-1

), to uzyska się poniższe

wyrażenie na stosunek wiarygodności dla sygnału o nieznanych parametrach:

∫ ∫

∞

−

∞

−

...

)

,...,

(

)

,...,

(

)

,...,

(

...

)

,...,

(

(2.3)

Mianownik został wprowadzony pod znak całki, ponieważ nie zależy od

zmiennych całkowania. Wielkość

)

,...,

(

)

,...,

(

)

,...,

(

−

przedstawia szczególną wartość stosunku wiarygodności dla

ustalonych wartości parametrów. Innymi słowy jest to stosunek

wiarygodności dla sygnału całkowicie znanego. Z powodu

występowania w nim warunkowej gęstości prawdopodobieństwa

często określa się go jako warunkowy stosunek wiarygodności.

Natomiast stosunek wiarygodności sygnału z nieznanymi

parametrami uzyska się w wyniku uśrednienia warunkowego

stosunku wiarygodności względem wszystkich parametrów

występujących w warunku:

∫ ∫

∞

−

∞

−

...

)

,...,

(

)

,...,

(

...

)

,...,

(

Jeśli sygnał użyteczny odbierany jest razem z addytywnym szumem

białym o rozkładzie Gaussa, to biorąc pod uwagę rozważania z

punktu 1.5, a w szczególności wyrażenie (1.31) na stosunek

wiarygodności oraz (1.33) na statystykę G

i (1.34) na energię E

sygnału x(t) warunkowy stosunek wiarygodności można zapisać

jako:

)

(

exp(

)

(

exp(

)

(

)

(

)

(

exp(

)

,...,

(

−

⋅

−

∑

−

gdzie

∑

−

)

(

)

(

)

(

)

(

)

(

w których to wyrażeniach wartości x(k) zależą od wektora

parametrów

ββββ

, zaś

jest odchyleniem standardowym szumu.

Ostatecznie stosunek wiarygodności dla sygnału o nieznanych

parametrach będzie miał poniższą postać:

∫ ∫

∞

−

∞

−

⋅

...

)

(

exp(

)

(

exp(

)

,...,

(

...

)

,...,

(

(2.4)

2.3. STOSUNEK WIARYGODNOŚCI DLA SYGNAŁU

O NIEZNANEJ FAZIE

Pokazane w punkcie 2.2 zależności, zwłaszcza wyrażenie (2.4) posłużą do

określenia stosunku wiarygodności w przypadku, gdy odbierany użyteczny sygnał

echa ma przypadkową fazę. Po dokonaniu próbkowania jego dyskretny

odpowiednik można zapisać jak poniżej:

]

)

(

cos[

)

(

)

(

(2.5)

Zgodnie z rozważaniami z poprzedniego punktu zostanie przyjęty najmniej

korzystny rozkład prawdopodobieństwa fazy tzn. równomierny:

)

;

)

(

)

(

max

∈

dla

(2.6)

Przekształcając postać (2.5) z zastosowaniem wzoru na sumę

cosinusa sumy, otrzymuje się, że:

sin

)]

(

sin[

)

(

cos

)]

(

cos[

)

(

)

(

−

Oznaczając

)]

(

sin[

)

(

)

(

)]

(

cos[

)

(

)

(

−

można zapisać x(n,

) jak niżej:

sin

)

(

cos

)

(

)

(

Wobec tego statystyka G

może być zapisana jako

∑

−

sin

cos

)

(

)

(

)

(

gdzie

∑

−

)

(

)

(

)

(

)

(

Jeśli potraktować z

i z

jako współrzędne punktu w prostokątnym

układzie współrzędnych, to można wprowadzić współrzędne

biegunowem Z i

określone jako

(2.7)

sin

cos

i wtedy otrzyma się ostateczną postać statystyki G

)

cos(

)

sin

cos

(cos

−

(2.8)

Pozostaje znależć jeszcze energię sygnału E

(

). Przyglądając się

wartości poniższej całki z funkcji cos

(ax+b):

∫

)

(

sin

)

(

cos

która odpowiada sumie dla sygnałów z ciągłym czasem, widać że

jeśli amplituda X(n) oraz faza

(n) sygnału x(n) zmienia się

nieznacznie podczas okresu wielkiej częstotliwości, to energia

sygnału nie zależy praktycznie od fazy i

)

(

Znając wartość G

(

) i E

(

) można korzystając z wyrażenia (2.4)

określić wartość stosunku wiarygodności dla sygnału o nieznanej

fazie:

∫

−

⋅

−

⋅

−

⋅

−

⋅

−

]

)

cos(

exp[

)

exp(

]

)

cos(

exp[

)

exp(

)

,...,

(

Wiedząc, że

)

(

)

cos(

∫

−

gdzie

)

(

jest zmodyfikowaną funkcją Bessela pierwszego rodzaju

zerowego rzędu (rys. 2.1) dostaje się ostateczne wyrażenie na

stosunek wiarygodności dla sygnału z nieznaną fazą:

)

(

)

exp(

)

,...,

(

⋅

−

(2.9)

Rys. 2.1. Wykres zmodyfikowanej funkcji Bessela pierwszego rodzaju zerowego rzędu

Ponieważ zmodyfikowana funkcja Bessela pierwszego rodzaju

zerowego rzędu jest ściśle monotoniczna, to decyzję o obecności

sygnału x w przebiegu y można podejmować porównując wartość

statystyki Z z ustaloną wartością Z

. Jeśli Z<Z

, to przyjmuje się

prawdziwość hipotezy H

, zaś gdy Z>Z

, to wybierana jest hipoteza

. Wielkość Z

ustalana jest w zależności od porządanego przez

obserwatora poziomu prawdopodobieństwa fałszywego alarmu F.

Schemat blokowy urządzenia decyzyjnego działającego w oparciu o

statystykę Z jest pokazany na rysunku 2.2. Na układy mnożące

podaje się jako napięcia odniesienia dwa przebiegi kwadraturowe:

(n) oraz x

(n) (przesunięte względem siebie w fazie o

/2).

Istnienie dwóch kanałów kwadraturowych wyklucza możliwość straty

sygnału użytecznego na skutek nieznajomości jego fazy. Jeśli

przykładowo nie powoduje on wzrostu wartości z

wskutek

przesunięcia fazy o

/2 względem przebiegu odniesienia w

pierwszym kanale, to na pewno wywoła wzrost wartości z

w drugim

kanale. Jak widać, rezultat odbioru sygnału w detektorze z rysunku

2.2 nie zależy od przypadkowej fazy początkowej sygnału.

Rys. 2.2. Schemat blokowy układu decyzyjnego z dwoma kanałami kwadraturowymi

Układ, o którym przed chwilą była mowa jest optymalny jedynie,

gdy znany jest czas pojawienia się (opóźnienia) sygnału

użytecznego. Decyzja o istnieniu sygnału z nieznanym czasem

pojawienia się może być podjęta, jeśli sprawdzone zostaną hipotezy

o istnieniu tego sygnału w różnych momentach czasu. W ten sposób

dochodzi się do układu wielokanałowego, w którym każdy kanał ma

postać, jak na rys. 2.2, lecz “dostrojony” jest do innego czasu

opóźnienia. Odstępy między tymi czasami nie przekraczają

rozdzielczości odległościowej radaru, jakiej rzyczy sobie obserwator.

Schemat blokowy układu pokazany jest na rysunku 2.3. Takie

wielokanałowe układy mogą być użyte nie tylko dla sygnałów echa

różniących się czasem opóźnienia (odległość obiektu od radaru), ale

także częstotliwością nośną na skutek efektu Dopplera (prędkość

radialna obiektu względem radaru).

Rys. 2.3. Schemat blokowy wielokanałowego układu detekcyjnego dla różnych czasów

opóźnienia sygnału echa

2.4. ZASTOSOWANIE FILTRU DOPASOWANEGO W DETEKCJI

SYGNAŁÓW Z NIEZNANYM CZASEM OPÓŹNIENIA

W poprzednich punktach zostało pokazane, że decyzja o tym, która z

dwóch hipotez jest prawdziwa tzn.:

: y(t)=n(t)

−

odebrano tylko zakłócenia w postaci białego szumu o rozkładzie

Gaussa, lub

: y(t)=x(t)+n(t)

−

na tle wyżej wspomnianych zakłóceń odebrano sygnał echa

odbity od obiektu w powietrzu, podejmowana jest (ogólnie) przez porównanie

wartości stosunku wiarygodności

(y) z pewnym, ustalonym progiem

l .

Ponieważ stosunek wiarygodności jest funkcją monotoniczną, to wystarczy

porównywać wartość statystyki G

o postaci

∑

−

)

(

)

(

(2.10)

z progiem G

, jeśli detekuje się sygnały całkowicie znane, lub wartość statystyki

w której

∑

−

)

(

)

(

)

(

)

(

z progiem Z

, jeśli trzeba wykryć sygnał o nieznanej fazie. Czynniki y(k)=y(t=t

)

są dyskretnymi wartościami sygnału y(t) obserwowanego w przedziale czasu

∈〈

0;T

〉

. Budowa układu wielokanałowego z rys. 2.3, w którym każdy kanał

wykrywa sygnał o różnym czasie opóźnienia n

jest nieekonomiczna.

Statystykę G

(czy z

i z

) dla dowolnego n

byłoby wygodnie otrzymać

przepuszczając próbki sygnału y(n) przez jeden filtr liniowy. Filtr powinien

obliczać statystykę G

dla dowolnego czasu opóźnienia sygnału echa x(n),

który można zapisać jak poniżej:







−

≤

)

(

)

(

dla

(2.11)

Ciąg próbek u(n) reprezentuje sygnał echa z zerowym czasem opóźnienia.

Zakłada się, że postać ciągu jest znana. Korzystając z u(n) można zapisać

wyrażenie na statystykę G

poniższą zależnością:

∑

−

)

(

)

(

(2.12)

W celu uproszczenia analizy (ale bez skutków dla końcowych wniosków) można

przyjąć zerowy czas opóźnienia n

=0. Ponieważ sygnał x(n) jest różny od zera

tylko w pewnym przedziale, to można zmodyfikować sumę G

)

(

)

(

)

(

)

(

−

∑

(2.13)

Sygnał na wyjściu filtru o skończonej odpowiedzi impulsowej (SOI) mającej

długość (N

+1) próbek jest opisany przez splot dyskretny:

∑

−

)

(

)

(

)

(

(2.14)

Wartość sygnału na wyjściu filtru w chwili N

wynosi:

∑

−

)

(

)

(

)

(

)

(

)

(

(2.15)

Poszukiwany jest taki filtr, aby próbka na jego wyjściu w chwili N

była równa

wartości statystyki G

. Porównując powyższe wyrażenie z sumą na statystykę

(wyr. 2.13) dla sygnału z zerowym czasem opóźnienia, widać ,że:

∑

−

)

(

)

(

)

(

)

(

)

(

)

(

dla

Wprowadzając nową zmienną n=N

-k otrzymuje się ostateczne wyrażenie na

odpowiedź impulsową filtru:

)

(

)

(

−

(2.16)

Filtr o takiej odpowiedzi impulsowej jest filtrem dopasowanym do sygnału

u(n). Jego odpowiedź impulsowa ma postać lustrzanego odbicia sygnału u(n)

względem chwili n=0 i przesunięta w prawo o N

, lub więcej próbek.

Przesunięcie nie może być jednak mniejsze od N

, ponieważ filtr będzie wtedy

nieprzyczynowy. Na wyjściu filtru będą pojawiać się opóźnione o N

próbek

wartości statystyki G

(lub z

czy z

, w zależności od tego, do jakiego sygnału

jest dopasowany filtr) dla kolejnych czasów opóźnienia n

. Im większe n

(czyli

większa odległość obiektu od radaru), tym później na wyjściu filtru pojawi się

wartość równa G

2.5. ZASTOSOWANIE ZESPOLONYCH AMPLITUD W

ANALIZIE PROCESU DETEKCJI

Sygnał echa odbity od obiektu w powietrzu i docierający do stacji

radarowej ma postać sumy

)

(

)

(

)

(

gdzie x(t) jest sygnałem sinusoidalnym wielkiej częstotliwości, zaś n(t) szumem

białym, którego widmo obejmuje szeroki zakres. Jednakże w zagadnieniach

związanych z detekcją sygnałów o wąskim widmie w białych szumach można

założyć, że sygnał i szum przechodzą przez filtr, którego pasmo przenoszenia

jest znacznie szersze niż widmo sygnału x(t). Zwykle rozważa się taki filtr jako

wąskopasmowy. Jego typ ma jednak niewielki wpływ na wykrywalność

sygnałów, ponieważ wycina on tylko składowe szumów z częstotliwościami nie

mającymi wpływu na sygnał użyteczny. Przebieg y(t) można traktować wtedy

jako wąskopasmowy tzn. jako przebieg sinusoidalny o częstotliwości takiej, jaką

ma x(t) i wolnozmiennej losowej aplitudzie i fazie [4]. Przykładowy szum

wąskopasmowy pokazano na poniższym wykresie 2.4.

Rys. 2.4. Przykładowy przebieg szumu wąskopasmowego

W przypadku sygnałów wąskopasmowych o wolnozmiennej amplitudzie i

fazie można znacznie uprościć analizę i obliczenia korzystając z pojęcia

sygnału analitycznego [2]. Niech dyskretny odpowiednik sygnału y(t) ma

poniższą postać:

)]

(

cos[

)

(

)

(

θ +

Tworzy się sygnał zespolony

(n), którego częścią rzeczywistą jest

y(n), zaś częścią urojoną transformata Hilberta sygnału y(n):

{ }

)

exp(

)

(

]

exp[

)]

(

exp[

)

(

)

(

)]

(

sin[

)

(

)]

(

cos[

)

(

)

(

)

(

)

(

)]

exp(

)

(

Re[

)

(

Y(n) jest amplitudą chwilową (obwiednią) sygnału y(n), a Y(n) zespoloną

amplitudą zawierającą fazę

(n):

)]

(

exp[

)

(

)

(

Analogicznie można stworzyć sygnał analityczny i zespolone amplitudy dla

innych przebiegów. Niech sygnał x(n) będzie oczekiwanym sygnałem echa o

pewnym czasie opóźnienia n

)

(

)

(

−

(2.17)

)]

(

)

cos[(

)

(

)]

(

cos[

)

(

)

(

−

Po wprowadzeniu zespolonych amplitud dostaje się:

)]

exp(

)

exp(

)

(

Re[

)]

exp(

)

(

Re[

−

Ponieważ równość ta jest spełniona dla wszytkich wartości zmiennej

zespolonej exp(jn

), to równe są także stojące przy tej zmiennej

współczynniki:

)

exp(

)

(

)

(

−

(2.18)

Analogicznie można zapisać zależności dla odpowiedzi impulsowej filtru

dopasowanego do sygnału u(n) o długości (N+1) próbek:

)

(

)

(

dop

−

(2.19)

Po wprowadzeniu zespolonych amplitud dostanie się że:

)]

exp(

)

exp(

)

(

Re[

)]

exp(

)

(

Re[

dop

−

Aby po prawej stronie powyższej równości znalazł się czynnik exp(jn

), należy

wprowadzić wartość sprzężoną, co nie wpłynie na porównywanie części

rzeczywistych:

)]

exp(

)

exp(

)

(

Re[

)]

exp(

)

(

Re[

dop

−

Ponieważ równość jest spełniona dla wszystkich wartości exp(jn

), to tak jak

poprzednio można zapisać:

)

exp(

)

(

)

(

dop

−

(2.20)

Powyższe wyrażenie określa związek między zespolonymi

amplitudami odpowiedzi impulsowej filtru dopasowanego i sygnału

oczekiwanego. Korzystając z zespolonych amplitud można zapisać

przebiegi na wejściu i wyjściu filtru dopasowanego:

∑

−

⋅

dop

)]

)

(

exp(

)

(

Re[

)]

exp(

)

(

Re[

)]

exp(

)

(

Re[

(2.21)

W celu przekształcenia tej zależności zostanie wykorzystany

poniższy związek:

]

Re[

⋅

Jeśli ponadto uwzględnić, że część rzeczywista jest równa sumie

części rzeczywistych, to zależność opisującą filtrację sygnału można

zapisać jak dalej:

)]

exp(

))

(

)

(

Re[(

)]

)

(

exp(

)

(

)

exp(

)

(

)

(

exp(

)

(

)

exp(

)

(

Re[

)]

exp(

)

(

Re[

dop

−

∑

(2.22)

gdzie:

)

(

exp(

)

(

)

(

)

(

)

(

)

(

)

(

dop

−

∑

Porównując ze sobą zespolone amplitudy wychodzi, że:

)

(

)

(

)

(

dop

(2.23)

Jeżeli częstotliwość sygnału użytecznego

jest dostatecznie duża w stosunku

do szerokości widma sygnału y(n), to amplitudy zespolone są funkcjami

wolnozmiennymi w stosunku do czynnika szybko zmieniającego się w czasie

exp(

−

jn2

). Podczas sumowania zauważa się, że W

(n)<< W

(n) i można

składową W

(n) pominąć:

∑

−

≅

dop

)

(

)

(

)

(

)

(

(2.24)

Korzystając z wcześniejszego przedstawienia h

dop

(n) za pomocą u(n)

(wyrażenia 2.19 i 2.20) dostaje się, że:

∑

−

≅

dop

)

(

)

exp(

)

(

)

(

(2.25)

Biorąc moduł zespolonej amplitudy W

dop

(n) otrzyma się wyrażenie na obwiednię

sygnału na wyjściu filtru dopasowanego:

∑

−

dop

)

(

)

(

)

(

)

(

(2.26)

Poczynione dotąd rozważania wraz z zależnościami z punktu 2.3 można

zastosować do zaprojektowania detektora z filtrem dopasowanym dla sygnałów

z nieznaną fazą. Decyzja odnośnie obecności sygnału echa radarowego będzie

podejmowana poprzez porównanie wartości statystyki Z i progu Z

określonego

przez poziom prawdopodobieństwa fałszywego alarmu F.

Wyrażenia na x

(n) oraz x

(n) z punktu 2.3 po wprowadzeniu zespolonych

amplitud można zapisać jak poniżej:

)]

exp(

)

(

Im[

)]

exp(

)

(

Im[

)]

(

sin[

)

(

)

(

)]

exp(

)

(

Re[

)]

exp(

)

(

Re[

)]

(

cos[

)

(

)

(

−

(2.27)

gdzie

)]

(

exp[

)

(

)

(

Wobec powyższego wyrażenia na z

i z

z punktu 2.3 będą miały postać (sygnał

y(n) obserwowany jest w zakresie n od 0 do N

-1):

∑

−

)]

exp(

)

(

)

(

Im[

)]

exp(

)

(

)

(

Re[

(2.28)

Wiedząc, że dla dowolnej liczby zespolonej a=Re[a]+jIm[a] słuszna jest

równość

(

) (

)

można zapisać wyrażenie na statystykę Z:

∑

−

)

exp(

)

(

)

(

(2.29)

Wykorzystując następujący związek, że

)

exp(

)

(

)

exp(

)

(

)]

(

cos[

)

(

)

(

−

uzyskuje się

∑

−













−

)

exp(

)

(

)

exp(

)

(

)

exp(

)

(

)

exp(

)

(

Drugi czynnik jest szybkozmienny (z częstotliwością 2

) w porównaniu do

pierwszego i podobnie jak poprzednio można go pominąć

∑

−

≅

)

(

)

(

(2.30)

Uwzględniając zależności (2.17) i (2.18), że x(n)=u(n-n

) dostaje się:

∑

−

≅

)

exp(

)

(

)

(

)

(

(2.31)

Po wyniesieniu przed sumę czynnika exp(jn

) i uwzględnieniu, że jego moduł

równy jest jeden, ostatecznie wychodzi, że:

∑

−

≅

)

(

)

(

)

(

(2.32)

Ponieważ sygnał x(n) jest różny od zera tylko w pewnym fragmencie całego

przedziału obserwacji, to w sumie na Z(n

) wystarczy wziąć pod uwagę tylko

(N+1) próbek sygnału użytecznego tzn. od n

do n

+N:

∑

−

≅

)

(

)

(

)

(

(2.33)

Wartość obwiedni sygnału na wyjściu filtru dopasowanego (patrz wzór 2.26) w

chwili n=(N+n

) przedstawia poniższa zależność:

∑

−

dop

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

(2.34)

z której porównawszy wyrażenia objęte znakiem sumy z (2.33) widać, że:

)

(

)

(

dop

(2.35)

tzn. amplituda chwilowa przebiegu na wyjściu filtru dopasowanego jest równa

wartości statystyki Z(n

) dla każdej wielkości czasu opóźnienia sygnału n

czyli inaczej: dla wszystkich odległości obiektu od stacji radarowej [5]. W ten

sposób dla utworzenia wielkości Z(n

) wystarczy jeden kanał z filtrem

dopasowanym, zaś w celu uzyskania amplitudy sygnału na wyjściu filtru należy

użyć detektora obwiedni. Schemat blokowy takiego urządzenia pokazany jest

na rysunku 2.5 na następnej stronie.

Rys. 2.5. Schemat blokowy jednokanałowego detektora z filtrem dopasowanym dla sygnału o

nieznanej fazie

Pozostaje jeszcze znalezienie funkcji gęstości prawdopodobieństwa

statystyki Z(n

), czyli rozkład wartości przebiegu będącego obwiednią sygnału

na wyjściu filtru dopasowanego.

W przypadku, gdy prawdziwa jest hipoteza H

, to sygnał na wejściu filtru

składa się tylko z szumu białego gaussowskiego o zerowej wartości

oczekiwanej i wariancji

. Przejście sygnału przez filtr liniowy nie zmienia typu

rozkładu, tak więc przebieg na wyjściu filtru dopasowanego również będzie

szumem gaussowskim (lecz już nie białym). Wartość oczekiwana pozostanie

zerowa, zaś zmieni się wariancja:

dop

)

(

)

(

]

[

−

∑

gdzie E

jest energią sygnału użytecznego. Obwiednia sygnału w(n) za filtrem i

jej rozkład prawdopodobieństwa zostanie wyznaczony z wykorzystaniem

pojęcia sygnału analitycznego. Sygnał analityczny skojarzony z przebiegiem

w(n) ma postać:

)}

(

{

)

(

)

(

Krzywa

)

(

)

(

będzie poszukiwaną obwiednią, a

))

(

)

(

poszukiwaną gęstością prawdopodobieństwa obwiedni. Z własności sygnału

analitycznego wynika, że zmienne losowe w(n) i H{w(n)} są nieskorelowane dla

danego n. Poza tym ponieważ H{w(n)} otrzymuje się w wyniku liniowej operacji

nad w(n

), to musi mieć taki sam rozkład, czyli rozkład Gaussa. W

szczególności widmo mocy i wariancja obu przebiegów będą jednakowe.

Skoro tak, to łączny rozkład gęstości prawdopodobieństwa zmiennych

losowych: w i H{w} można zapisać poniższym wyrażeniem

)

exp(

)

]

[

exp(

])

[

(

)

(

])

[

(

πσ

−

. (2.36)

Poszukiwany jest rozkład zmiennej losowej

]

[

dla Z

≥

0. Dystrybuanta zmiennej losowej Z równa jest masie

prawdopodobieństwa na obszarze będącym kołem o promieniu Z:

∫∫

≤

dwdH

]

[

]

[

]

[

])

[

(

)

(

Przechodząc na współrzędne biegunowe

sin

]

[

cos

≡

wychodzi, że

∫

−

]

[

)

exp(

)

exp(

]

[

)

(

)

(

πσ

(2.37)

Po zróżniczkowaniu uzyskanego wyrażenia na F

(Z) dostaje się wyrażenie na

gęstość prawdopodobieństwa statystyki Z w przypadku prawdziwości hipotezy

, czyli tzw. rozkład Rayleigha [9]:

exp(

)

(

≥

−

(2.38)

Krzywa rozkładu Rayleigha pokazana jest na wykresie 2.6 na kolejnej stronie.

Rys. 2.6. Gęstość prawdopodobieństwa obwiedni szumu wąskopasmowego czyli tzw. rozkład

Rayleya.

jest odchyleniem standardowym szumu

W przypadku, kiedy prawdziwa jest hipoteza H

, to sygnał na wejściu filtru

dopasowanego składa się z sygnału echa oraz addytywnego szumu białego

gaussowskiego o zerowej wartości oczekiwanej i wariancji

. Na wyjściu filtru

pojawi się sygnał będący sumą szumu takiego, jak w przypadku prawdziwości

hipotezy H

i przefiltrowanego sygnału użytecznego. Aby wyznaczyć rozkład

(Z(n)) należy (jak na początku niniejszego punktu 2.5) potraktować zakłócenia

na wyjściu filtru dopasowanego jako szum wąskopasmowy, o częstotliwości

nośnej

. Skoro tak, to można go zapisać poniższym wyrażeniem:

)

sin(

)

(

)

cos(

)

(

))

(

cos(

)

(

−

(2.39)

gdzie

)

(

)

(

]

[

)

(

)

(

= ψ

)

(

)

(

)

(

Sygnał wypadkowy na wyjściu filtru można przedstawić jako sumę prawej

strony wyrażenia (2.39) i Acos(n

) [6]:

)

sin(

)

(

)

cos(

)

(

)

cos(

))

(

cos(

)

(

−

(2.40)

Początek obserwacji przebiegu (2.40) przyjęto tak, aby faza była równa zero i

wtedy składnik użyteczny ma postać Acos(n

). Ponieważ (2.40) reprezentuje

sygnał na wyjściu filtru dopasowanego, to wartość A nie jest amplitudą echa

radarowego x(n). W szczególności jeśli x(n)=u(n-n

), to składnik użyteczny na

wyjściu filtru dopasowanego o długości (N+1) próbek w chwili n=(N+n

) będzie

równy energii sygnału u(n):

∑

−

dop

)

(

)

(

)

(

)

(

)

(

)

(

Porównując stronami wyrażenie (2.40) wychodzi, że:

)

(

))

(

sin(

)

(

)

(

)

(

))

(

cos(

)

(

)

(

−

(2.41)

Poszukiwany jest rozkład obwiedni sygnału za filtrem, czyli W(n), gdzie:

i W

≥

0. Podstawiając (2.41) do wyrażenia o postaci (2.36) otrzymuje się:

)

(

exp(

)

exp(

)

(

πσ

−

(2.42)

Gęstość prawdopodobieństwa zmiennej W: p

(W) można wyliczyć albo

różniczkując jej dystrybuantę F

(W), albo z poniższej zależności:

)

(

)

(

)

(

∫∫

∆

gdzie

∆

D jest obszarem płaszczyzny n

takim, że dla W

≥

Obszar ten jest pierścieniem kołowym o wewnętrznym promieniu W i

zewnętrznym W+dW. Wprowadzając współrzędne biegunowe

WdWd

≡

sin

cos

dostaje się:

∫

−

πσ

)

sin

(

)

cos

(

exp(

)

(

WdWd

a więc

πσ

)

cos

exp(

)

exp(

)

(

∫

⋅

−

(2.43)

Ostatnią całkę można wyrazić przez zmodyfikowaną funkcją Bessela rzędu

zerowego:

∫

cos

)

(

i ostatecznie rozkład prawdopodobieństwa obwiedni sygnału (czyli rozkład

statystyki Z) na wyjściu filtru dopasowanego przy prawdziwości hipotezy H

(tzw. rozkład Rice’a) ma poniższą postać:

(

)

exp(

)

(

≥

⋅

−

(2.44)

Krzywe rozkładu Rice’a dla różnych wartości energii sygnału użytecznego E

pokazane są na następnej stronie na wykresie 2.7.

Rys. 2.7. Gęstość prawdopodobieństwa obwiedni sumy sygnału i szumu wąskopasmowego

czyli tzw. rozkład Rice’a dla różnych wartości stosunku sygnału do szumu:

−

rozkład Rayleya

−

rozkład Rice’a

3. INTEGRACJA IMPULSÓW I STABILIZACJA

POZIOMU FAŁSZYWEGO ALARMU

3.1. ZASADA DZIAŁANIA RZECZYWISTYCH SYSTEMÓW

RADIOLOKACYJNYCH

W poprzednich rozdziałach opisano podstawy detekcji oraz przedstawiono

sytuację, w której do wykrywania obiektów w przestrzeni wykorzystuje się tylko

jeden sygnał echa. Przebieg wejściowy był obserwowany tylko w pojedynczym

przedziale czasu 0

T, w którym sygnał echa występował lub nie istniał.

Odpowiada to przypadkowi, gdy obserwator nakierowuje antenę radaru na jakiś

wybrany wycinek przestrzeni i w jego kierunku wysyła jeden impuls.

W wykorzystywanych rzeczywistych urządzeniach radiolokacyjnych

można nadawać wiele impulsów w kierunku obiektu i obserwować pewną liczbę

przedziałów, a w każdym z nich oczekiwać sygnału echa, jeśli tylko cel

rzeczywiście istnieje. Antena radaru obraca się ze stałą prędkością kątową i

analizuje otaczającą przestrzeń specjalnie ukształtowaną wiązką sondującą.

Przykład takiej wiązki jest pokazany na rysunku 3.1.

Rys. 3.1. Przykładowa wiązka radarowa (modulowana charakterystyką anteny) skanująca

przestrzeń dookoła stacji radiolokacyjnej

Stacja radarowa co pewien stały interwał czasowy generuje sygnał w postaci

impulsu wielkiej częstotliwości s

(t). Charakterystyka promieniowania anteny w

kierunku obiektu (przedstawiona na rys 3.1) jest określona pewną funkcją czasu

(t). Tak więc sygnał sondujący promieniowany w kierunku obiektu jest

opisany przez iloczyn H

(t)s

(t). Natężenie pola sygnału echa odbitego od

obiektu i przychodzącego do stacji radarowej z opóźnieniem t

wynosi z

dokładnością do stałej H

(t-t

), gdzie

przy czym r jest odległością celu od stacji, a c szybkością

rozchodzenia się fal elektromagnetycznych (czyli prędkością

światła).

W rezultacie ruchu anteny odbiorczej jej charakterystyka jest także funkcją

czasu. Oznaczona zostanie przez H

(t). Sygnał użyteczny x(t) przychodzący na

wejście odbiornika jest modulowany dodatkowo funkcją H

(t) i może być

zapisany w postaci:

)

(

)

(

)

(

)

(

−

⋅

−

⋅

−

Dla typowych jednoantenowych radiolokatorów H

(t)=H

(t)=H(t). Równocześnie

zmiana charakterystyki antenowej H(t) w czasie t

jest tak znikoma, że można

ją zaniedbać [5], a tym samym

)

(

)

(

)

(

−

⋅

Tak więc radiolokacyjny sygnał echa pochodzący od obiektu jest

modulowanym przebiegiem wielkiej częstotoliwości x(t), przy czym

modulacja ta jest wynikiem nie tylko celowego działania

obserwatora, ale i charakterystyki anteny nadawczo-odbiorczej.

Modulacja ta w ustalonych warunkach pracy radaru nie ulega

zmianie.

Modulowany, pojedynczy sygnał wielkiej częstotliwości x(t)

można zapisać w postaci

))

(

cos(

)

(

)

(

ω +

(3.1)

Sygnał taki, oraz powyższy zapis wykorzystywane były w analizie

detekcji w poprzednim rozdziale. W zapisie tym

jest pulsacją

nośną, zaś X(t) i

(t) funkcjami wolnozmiennymi w stosunku do

cos(

t) określającymi amplitudową i fazową modulację sygnału

echa. Rysunek 3.2 na następnej stronie pokazuje kształt takiego

impulsu, jego obwiednię X(t), oraz przebieg fazy początkowej

(t).

Ponieważ nie występuje modulacja fazy, to

(t)=const.

Na kolejnym rysunku 3.3 przedstawiono sygnał echa w postaci

paczki impulsów wielkiej częstotliwości, modulowanej

charakterystyką jednostajnie obracającej się anteny. Można go

zapisać jako sumę wyrażeń typu (3.1):

∑

))

(

cos(

)

(

)

(

(3.2)

Rysunek 3.3b przedstawia obwiednię paczki, zaś 3.3c przebieg fazy

początkowej jako funkcję czasu.

Rys. 3.2. Przykład impulsu wielkiej częstotliwości (a), jego obwiedni (b) oraz przebieg fazy

początkowej (c) w funkcji czasu

Rys. 3.3. Koherentna paczka impulsów echa modulowana charakterystyką anteny (a), jej

obwiednia (b) oraz przebieg fazy początkowej (c) w funkcji czasu dla nieruchomego obiektu

Jak widać z wykresu fazy, paczka jest echem odbitym od obiektu

nieruchomego, a impulsy są generowane przez nadajnik z taką

samą fazą początkową. Taka paczka, z jednoznacznie określonym

przebiegiem fazy

(t) dla każdego impulsu nosi nazwę paczki

koherentnej.

Inaczej wygląda sytuacja, gdy nadajnik generuje impulsy w.cz.

z przypadkową fazą początkową. Jeśli w chwili nadawania fazy nie

zostaną zapamiętane, tak aby umożliwić późniejszą korektę w

sygnale echa, to w paczkach sygnałów odbieranych trzeba

traktować je jako przypadkowe i z góry nieznane. Paczka taka

zwana jest paczką niekocherentną. Można ją zapisać wyrażeniem

(3.3):

∑

)

(

cos(

)

(

)

(

(3.3)

W wyrażeniu tym parametr

jest zmienną losową o rozkładzie

równomiernym w przedziale od 0 do 2

(punkty 2.1-2.3, rozdz. 2),

przy czym dla każdego impulsu zmienne te są niezależne od

zmiennych losowych z innych impulsów.

W związku z tym, że antena radaru obraca się nieprzerwanie i

co pewien okres czasu (czyli co obrót o jakiś ułamek kąta pełnego)

emituje w przestrzeń wokół siebie impuls sondujący, obszar

powietrzny otaczający stację radarową po jednym obrocie anteny o

pełny kąt może być przedstawiony jako wykres trójwymiarowy. Na

jednej z osi będzie kąt azymutalny

w zakresie od 0 do 360 stopni,

na drugiej zasięg radaru od 0 do r

max

, zaś na osi pionowej moc

odebranego sygnału wraz z zakłóceniami. Zgodnie z teorią z

rozdziałów pierwszego i drugiego jeśli moc sygnału przekroczy

ustalony próg, to znaczy, że gdzieś tam w powietrzu znajduje się

cel. Ponieważ stacja jest w stanie emitować impulsy tak często, że

nawet dla szybko lecącego obiektu duża ich liczba trafi weń i odbije

się od niego przy niewielkiej zmianie jego położenia, to można

zwiększyć prawdopodobieństwo detekcji poprzez podejmowanie

decyzji na podstawie większej niż jeden liczby impulsów.

Urządzenie takie, zwane integratorem będzie wykonywać różne

operacje na impulsach, a następnie na podstawie uzyskanych z

impulsów wartości (na przykład średnia arytmetyczna) poprzez

porównanie ich z pewnym progiem będzie podejmować decyzję

odnośnie obecności obiektu w przestrzeni. Najprostszy taki

detektor/integrator typu “ruchome okno” ma strukturę pokazaną na

rysunku 3.4 [11].

Rys. 3.4. Schemat blokowy integratora typu „ruchowe okno”. ‘C’ oznacza porównywanie z

progiem decyzyjnym T

Jego działanie polega na sumowaniu ustalonej przez

obserwatora liczby

impulsów dla każdej komórki odległościowej.

Suma ostatnich

impulsów jest następnie porównywana z pewnym

progiem. Jeśli suma przekroczy próg, to znaczy, że w wycinku

obszaru, w którym wyemitowano te

impulsów znajduje się cel.

Działanie integratora ze schematu 3.4 opisuje poniższe wyrażenie:

−

(3.4)

gdzie S

oznacza sumę wartości ostatnich

impulsów, a x

wartość

-tego impulsu. Realizuje on zatem “kroczącą” sumę amplitud

impulsów, dzięki czemu w ramach całej paczki następuje narastanie

wartości sygnału wyjściowego. Jednocześnie szum, jako sygnał

przypadkowy nie ulega tak efektywnemu sumowaniu, co prowadzi

do poprawy stosunku sygnału do szumu na wyjściu integratora.

Punktem wyjścia do analizy działania integratora i wykrywania

będzie tak, jak dla pojedynczego impulsu, stosunek wiarygodności.

3.2. STOSUNEK WIARYGODNOŚCI DLA SYGNAŁU W POSTACI

PACZKI IMPULSÓW W.CZ. O NIEZNANYCH FAZACH

Dyskretny odpowiednik sygnału w postaci paczki impulsów wielkiej

częstotliwości z przypadkowymi fazami początkowymi można zapisać

poniższym wyrażeniem:

∑

]

)

(

cos[

)

(

,...)

(

(3.5)

Amplitudy impulsów są wielkościami zdeterminowanymi,

zaś fazy są niezależnymi zmiennymi losowowymi o

rozkładach równomiernych w przedziale od 0 do 2

. Taki

najmniej korzystny rozkład fazy został przyjęty także w

punkcie 2.3, gdzie liczony był stosunek wiarygodności dla

pojedynczego impulsu. Łączna funkcja gęstości

prawdopodobieństwa fazy będzie równa

)...

(

)

(

,...)

(

(3.6)

gdzie

)

;

...

)

(

)

(

∈

dla

W celu obliczenia stosunku wiarygodności wykorzystane

zostaną ogólne wyrażenie na stosunek wiarygodności

sygnału o nieznanych parametrach (2.4), oraz zależności z

punktu 2.3.

Wykorzystując, podobnie jak w punkcie 2.3 wzór na

cosinus sumy wyrażenie (3.5) można sprowadzić do

postaci:

∑

]

sin

)

(

cos

)

(

[

,...)

(

w której to:

)]

(

sin[

)

(

)

(

)]

(

cos[

)

(

)

(

−

Statystyka G

może być zapisana jako

∑

−

]

sin

cos

[

)

(

,...)

(

,...)

(

przy czym

∑

−

)

(

)

(

)

(

)

(

Wprowadzając, analogicznie jak w punkcie 2.3 poniższe

oznaczenia

(3.7)

sin

cos

statystykę G

można przedstawić poniższym wyrażeniem:

∑

−

)]

cos(

[

,...)

(

(3.8)

Jeśli przyjąć, że rozpatrywana paczka składa się z

impulsów nie pokrywających się wzajemnie w czasie, to

energię całej paczki można obliczyć sumując energię

poszczególnych impulsów. Biorąc pod uwagę to, że

energia dla pojedynczego impulsu (punkt 2.3) nie jest

zależna od fazy można zapisać, że:

∑

)

(

,...)

(

(3.9)

gdzie E

jest energią

-tego impulsu. Teraz korzystając ze

wzoru (2.4) na stosunek wiarygodności dla sygnału o

nieznanych parametrach można określić wartość tego

ilorazu dla paczki impulsów o nieznanych fazach.

Można go zapisać w postaci poniższego iloczynu:

∏

∫

−

]

)

cos(

exp[

)

exp(

(3.10)

Wiedząc, że

)

(

)

cos(

∫

−

jest zmodyfikowaną funkcją Bessela pierwszego rodzaju

zerowego rzędu otrzymuje się ostateczne wyrażenie na

stosunek wiarygodności:

∏

−

)

(

)

exp(

(3.11)

Porównując wyrażenie (3.11) z wyrażeniem (2.9) z punktu

2.3 widać, że stosunek wiarygodności dla sygnału w

postaci paczki impulsów nie pokrywających się

wzajemnie w czasie z przypadkowymi i niezależnymi

fazami jest iloczynem stosunków wiarygodności

obliczonych dla każdego impulsu z paczki.

3.3. DETEKTOR DLA SYGNAŁU W POSTACI PACZKI

IMPULSÓW W.CZ. O NIEZNANYCH FAZACH

W punkcie tym zostanie przeanalizowany integrator odbierający

sygnał w postaci paczki impulsów wielkiej częstotliwości o

przypadkowych (i niezależnych) fazach początkowych, a następnie

sumujący

impulsów i porównujący sumę z pewnym progiem

Punktem wyjścia dla rozważań będzie wyprowadzone w poprzednim

punkcie 3.2 wyrażenie (3.11) na stosunek wiarygodności:

∏

−

)

(

)

exp(

Po zlogarytmowaniu powyższego wyrażenia otrzyma się, że:

∑

−

)

(

(3.12)

Operacje mnożenia zostały sprowadzone do operacji dodawania.

Ponieważ funkcja logarytmiczna jest ściśle monotoniczna, to

zamiast liczyć stosunek wiarygodności

i porównywać go z progiem

można określić wartość

i porównać ją z odpowiednim progiem

ln l

. Jeśli

, znaczy że prawdziwa jest hipoteza H

i wykryto

obiekt. W przeciwnym przypadku prawdziwa jest hipoteza H

Podejście takie pozwala na znaczne uproszczenie budowy

detektorów. Detektor taki powinien wykonać operacje matematyczne

opisane poniższym wyrażeniem:

∑

)

(

Pierwszym etapem jest znalezienie wartości statystyki Z

dla

każdego

-tego impulsu. Zgodnie z rozważaniami z punktu 2.5 jeśli

sygnał y(n) jest obserwowany w zakresie n od 0 do N

-1, to można

ją zapisać wyrażeniem (2.30):

∑

−

≅

)

(

)

(

gdzie X

(n) i Y(n) są zespolonymi amplitudami przebiegów x

(n)

oraz

y(n). W powyższym wyrażeniu

)]

(

exp[

)

(

)

(

jest znaną funkcją zależną od położenia na osi czasu. Aby

zniwelować modulację amplitudy paczki przez charakterystykę

anteny radaru zostaną wprowadzone współczynniki wagowe S

. Dla

największego impulsu z paczki zostanie przyjęte S

max

=1. Ponadto

oznaczając przez n

moment wysłania

-tego impulsu sondującego,

a przez n

opi

czas jego opóźnienia (związany z odległością obiektu od

stacji radarowej) można zapisać funkcję X

(n) w postaci:

)

(

exp(

)

(

)

(

opi

−

⋅

−

⋅

(3.13)

Wtedy

gdzie dla każdego impulsu

∑

−

≅

)

(

)

(

opi

(3.14)

Porównując wyrażenie (3.14) z wyrażeniem (2.34) z punktu 2.5

widać, że wszystkie wielkości Z

można otrzymać za pomocą

detektora przedstawionego na rysunku 2.5. Składa się on z filtru

dopasowanego do impulsu u(n) oraz detektora obwiedni. Aby

otrzymać wielkość Z

należy go poszerzyć o układ wprowadzający

współczynniki wagowe S

. Tak więc w wyniku pierwszego etapu

obróbki paczki niekoherentnej impulsów wielkiej częstotliwości

otrzyma się paczkę obwiedni impulsów branych z odpowiednią

wagą.

Kolejny etap obróbki polega na obliczeniu wartości

)

(

dla każdego

oraz ich zsumowaniu. Wynik sumowania nie zależy

od początkowych faz impulsów w.cz. gdyż dodaje się ich obwiednie.

Całkowity układ obróbki jest przedstawiony na poniższym rysunku

3.5.

Rys. 3.5. Schemat blokowy układu optymalnej detekcji paczki impulsów wielkiej częstotliwości o

nieznanych fazach

Składa się on z filtru dopasowanego do impulsu u(n), detektora

obwiedni, układu mnożącego przez współczynniki obwiedni paczki

, bloku nieliniowego o charakterystyce

)

(

oraz sumatora, który

skupia w czasie niejednoczesne obwiednie impulsów, a następnie

dodaje je do siebie.

W sytuacjach krańcowych, gdy poziomy sygnałów (impulsów)

są albo bardzo słabe, albo bardzo silne można schemat z rysunku

3.5 znacznie uprościć. Jak wiadomo z właściwości funkcji Bessela,

dobrą aproksymacją funkcji

)

(

dla małych wartości argumentu

u=S

jest pierwszy człon szeregu potęgowego, na który można

ją rozłożyć [5]:

)

(

≅

dla

(3.15)

Oznacza to, że początkowa część jej przebiegu ma charakter

paraboliczny. Natomiast dla dużych wartości argumentu u funkcję

można zastąpić przebiegiem asymptotycznym:

)

(

≅

dla

(3.16)

co wskazuje na liniowy kształt jej wykresu na tym odcinku. Wykres

funkcji

)

(

jest pokazany na poniższym wykresie 3.6.

Rys. 3.6. Wykres funkcji

)

(

W związku z tym dla paczki impulsów o amplitudach małych w

porównaniu z szumami wychodzi, że:

∑

≅

)

(

(3.17)

zaś dla impulsów o dużych amplitudach

∑

≅

)

(

(3.18)

W ten sposób dodawanie logarytmów zostało zastąpione przez

sumowanie liniowych lub kwadratowych funkcji wielkości Z

. W

detektorze z rysunku 3.5 w przypadku dużych sygnałów w ogóle

zniknie blok nieliniowy o charakterystyce

)

(

, zaś dla małych

sygnałów zostanie zastąpiony układem mnożącym, na którego oba

wejścia zostaną podane kolejne S

3.4. INTEGRATOR BINARNY

W poprzednim punkcie przedstawiono detektor, który kumuluje odebrane

impulsy echa radarowego w jednej wartości. W zależności od poziomu sygnału

sumował obwiednie impulsów lub ich kwadraty i porównywał z wartością

progową. Jednak przedstawione w punkcie 3.3 rozwiązanie ma pewną wadę.

Jest wrażliwe na sygnały o bardzo dużej amplitudzie, nie będące rzeczywistym

echem. Mogą to być np. aktywne zakłócenia impulsowe wygenerowane celowo

przez wrogi obiekt, aby uniemożliwić pracę radaru [2]. Przykładowo integrator

kumulujący

impulsów w momencie, gdy cztery odebrane będą tylko

szumem, a piąty odpowiednio wysokim impulsem zakłócającym wywoła

fałszywy alarm. Wrażliwość na zbyt duże impulsy nie jest tylko cechą

integratora z punktu 3.3, ale wszystkich, w których kumulant, na podstawie

którego podejmuje się decyzję tworzony jest bezpośrednio z amplitudy sygnału.

Aby pozbyć się wspomnianego problemu można zastosować tzw.

integrator dwuprogowy. Zwany jest często integratorem binarnym lub

detektorem M-z-N, gdzie N oznacza ilość integrowanych impulsów w paczce.

Integrator ten ma wiele zalet. Po za tym, że nie jest wrażliwy na duże skoki

amplitudy pojedynczych impulsów, łatwo go zaimplementować oraz działa

bardzo dobrze gdy rozkład zakłóceń jest różny od rozkładu Rayleya [2, 11].

Jego schemat blokowy pokazany jest na poniższym rysunku 3.7.

Rys. 3.7. Schemat blokowy integratora binarnego (dwuprogowego). ‘C’ oznacza porównywanie

z progami decyzyjnymi

Jak to jest pokazane na schemacie impulsy wejściowe są kwantowane do

wartości 1 lub 0, w zależności od tego, czy przekroczony lub nie zostanie próg

. Zakłada się przy tym prostokątną obwiednię paczki impulsów tzn. każdy

impuls został przemnożony przez odpowiedni współczynnik S

niwelujący wpływ

charakterystyki anteny na obwiednię paczki. Ostatnie N jedynek i zer jest

sumowane i porównywane z drugim progiem T

=M. Mówiąc inaczej sumowanie

liniowe lub kwadratowe występujące w detektorach operujących bezpośrednio

na amplitudzie sygnału zostało zastąpione zliczaniem impulsów, które

przekroczyły próg T

3.5. FAŁSZYWY ALARM

W analizie detektorów otoczenie stacji radarowej przyjmuje się zwykle

jako znane i homogeniczne (jednorodne). Przy takim założeniu można

stosować w detektorach stały poziom progu wykrywania. Jednakże w

rzeczywistym środowisku, znajdują się liczne pofałdowania terenu, morza,

padają deszcze i występują gęste chmury. Wszystkie te elementy środowiska

odbijają sygnały radiolokacyjne, pochodzące zarówno z własnej stacji

radarowej, jak i z obcych. Przestrzeń wokół radaru pełna jest biernych zakłóceń

(ang. clutter) i przy stosowaniu stałego progu detekcji może pojawić się

olbrzymia liczba fałszywych alarmów. Te z kolei mogą nadmiernie przeciążyć

cyfrowy system obróbki danych. Rysunek 3.8 na następnej stronie pokazuje

wzrost poziomu fałszywego alarmu spowodowany wzrostem poziomu zakłóceń.

Jedną z metod ograniczenia fałszywych alarmów w systemach ze stałym

progiem jest jego podniesienie. Niestety takie rozwiązanie zmniejszy

wykrywalność obiektów znajdujących się w obszarze o niskim poziomie

zakłóceń biernych. Zamiast tego stosuje się inne metody stabilizacji poziomu

prawdopodobieństwa fałszywego alarmu (ang. Constant False Alarm Ratio

−

CFAR) polegające na stabilizacji poziomu zakłóceń lub na automatycznej

zmianie poziomu progu w zależności od poziomu zakłóceń albo, na kombinacji

obu tych metod.

Rys. 3.8. Wpływ wielkości zakłóceń na liczbę przekroczeń progu decyzyjnego i poziom

fałszywego alarmu: (a)

−

niski poziom zakłóceń, (b)

−

wysoki poziom zakłóceń

Obie metody opierają się na założeniu, że środowisko jest jednorodne na

małym obszarze wokół komórki odległościowej, która aktualnie jest testowana.

Próbki sygnału z komórek odległościowych otaczających komórkę testowaną

(zwane komórkami odniesienia) są statystycznie niezależne i mają takie same

rozkłady.

W przypadku układów CFAR z progiem automatycznie adaptującym się

do poziomu szumów przyjmuje się, że funkcja gęstości prawdopodobieństwa

szumu jest znana poza kilkoma parametrami. Komórki odniesienia otaczające

komórkę testową są wykorzystywane do estymacji tych nieznanych

parametrów, a na podstawie ich estymat ustalany jest poziom progu detekcji.

Decyzję o obecności obiektu podejmuje się jeśli poziom sygnału z komórki

testowej jest dostatecznie duży w porównaniu z komórkami referencyjnymi

Najprostszy i najszerzej stosowany jest układ CFAR z rysunku 3.9 na

stronie dalej, bazujący na uśrednianiu odległościowym. Detektor ten był szeroko

badany przez Finna i Johnsona [11].

Rys. 3.9. Schemat blokowy układu CFAR z uśrednianiem odległościowym i progiem

automatycznie adaptującym się do poziomu szumów. ‘C’ oznacza porównywanie z progiem

decyzyjnym. Przed układem CFAR umieszczono integrator

Jeśli zakłócenia mają rozkład Rayleya tzn.

)

exp(

)

(

−

tylko jeden parametr

musi być estymowany (

jest mocą szumu). Próg

detekcji można określić z poniższego wyrażenia:

⋅

∑

(3.19)

gdzie

jest estymatą parametru

, zaś n ilością komórek odniesienia.

Ponieważ próg T jest ustalany przez estymatę

która jest określona z pewnym

błędem, trzeba go podnieść nieco wyżej „na wszelki wypadek”, niż gdyby był

dokładnie znany a priori. Dodatkowy zapas poziomu progu powoduje

pogorszenie czułości detektora i jest traktowany jako straty wprowadzone przez

CFAR. Jakość estymacji parametru

zależy od tego, na bazie ilu komórek

referencyjnych obliczona zostanie estymata

. Dla małej ich ilości strata CFAR

jest duża z powodu kiepskiej jakości estymacji. Z kolei zbyt duża ich liczba

może spowodować naruszenie założenia o jednorodności obszaru tylko w

małym otoczeniu komórki testowej. Dobrą praktyczną zasadą jest użycie takiej

liczby komórek odniesienia, aby utrzymać straty CFAR poniżej 1dB [11]. W

przypadku wątpliwości, czy zakłócenia faktycznie mają rozkład Rayleya należy

zmodyfikować układ z rys. 3.9 usuwając integrator sprzed bloku CFAR. Lepiej

jest uśredniać odległościowo i porównywać z progiem pojedyncze impulsy oraz

zastosować opisany w punkcie 3.4 integrator binarny (dwuprogowy). Schemat

układu pokazany jest na rysunku 3.10.

Rys. 3.10. Schemat blokowy układu CFAR z uśrednianiem odległościowym i progiem

automatycznie adaptującym się do poziomu szumów wraz z integratorem binarnym. ‘C’

oznacza porównywanie z progiem decyzyjnym

Zastosowanie takiego rozwiązania z integratorem binarnym jest także

przydatne, gdy moc zakłóceń zmienia się z każdym impulsem.

Układy CFAR, których działanie polega na stabilizacji poziomu zakłóceń

także bazują na uśrednianiu odległościowym impulsów. Jednak w

przeciwieństwie do układów opisanych wyżej tu poziom progu jest cały czas

stały. Przykład takiego układu jest pokazany na kolejnej stronie na rysunku

3.11.

Rys. 3.11. Schemat blokowy układu CFAR z uśrednianiem odległościowym i stałym progiem

Tak jak w układach z rysunków 3.9 i 3.10 próbki sygnału radiolokacyjnego

wprowadzane są do rejestru przesuwnego. Następnie liczona jest wartość

średnia z próbek pochodzących z komórek odległościowych otaczających

obustronnie komórkę testową. W wyniku podzielenia wartości sygnału S

komórki testowej przez wartość średnią S

otrzymuje się w sygnale wyjściowym

wyj

stabilizowany poziom szumu. Przez pojęcie „stabilizowany” należy

rozumieć, że zmienia się on w stopniu mniejszym niż poziom szumu

wejściowego. Analogicznie, jak w poprzednim układzie stopień stabilizacji

szumu zależy od tego, na bazie ilu komórek odległościowych otaczających

komórkę testową liczy się wartość średnią.

Przedstawione na rysunkach 3.9

−

3.11 schematy blokowe układów CFAR

są najbardziej podstawowymi ich realizacjami. W praktyce układy te są w różny

sposób modyfikowane. Jedną z możliwych zmian jest wstawienie na wejściu

bloku przetwornika liniowo-logarytmicznego. Dzięki temu zapewnia się

poprawną pracę układu w szerokim zakresie wartości sygnału wejściowego i

zastępuje się operację dzielenia próbek prostszą operacją odejmowania [2].

We współczesnych radarach klasyczne układy CFAR stabilizujące poziom

szumów są wspomagane przez układy bazujące na zasadzie uśredniania

powierzchniowego. Znane są one pod nazwą map zakłóceń biernych (ang.

clutter map). Obszar wokół stacji radarowej dzieli się na komórki azymutalno-

odległościowe i dla każdej zapamiętuje poziom zakłóceń. Zakładając, że w

długim okresie czasu rozkład przestrzenny zakłóceń nie ulega zmianie, można

obliczyć średni poziom zakłóceń za okres pewnej liczby sondowań

otaczającego obszaru. Uzyskane tak wartości można wykorzystać do

ustawienia poziomu progu decyzyjnego. W ten sposób próg nie ma stałego

poziomu dla całego obszaru obserwowanego przez radar, ale jest dostosowany

do poziomu zakłóceń w każdej komórce rozróżnialności.

4. ANALIZA JAKOŚCI UKŁADÓW OPTYMALNEJ

DETEKCJI SYGNAŁÓW RADIOLOKACYJNYCH

Najlepszą w sensie statystycznym jakość wykrywania i estymacji

parametrów sygnałów radiolokacyjnych można uzyskać tylko poprzez ich

optymalną obróbkę. Jak wspomniano w rozdziale pierwszym, jakość detekcji

określa się poprzez prawdopodobieństwo poprawnego wykrycia D przy

ustalonym prawdopodobieństwie fałszywego alarmu F. Prawdopodobieństwo

detekcji określa się dla różnych stosunków energii sygnału użytecznego do

mocy zakłóceń. Zgodnie z kryterium Neymana-Pearsona dla stałego poziomu F

dąży się do maksymalizacji D. Obliczanie jakości układów optymalnych ma

znaczenie, ponieważ jakość ta ukazuje granicę, do której dąży się przybliżając

układ rzeczywisty do optymalnego. W poprzednich rozdziałach opisana została

budowa i sposób działania układów detekcji, które realizują obróbkę sygnałów

zbliżoną do optymalnej. W tym rozdziale zostaną przeanalizowane

poszczególne bloki systemu, wykonujące różne operacje na sygnale echa

radarowego zakłóconego szumem. Pokazane zostaną metody pozwalające

wyliczyć rozkłady prawdopodobieństw przebiegów na ich wyjściu w przypadku

odbierania jedynie zakłóceń oraz dla sygnału użytecznego wraz z addytywnymi

zakłóceniami. Ich zastosowanie wymaga wielu obliczeń, które w postaci

analitycznej są żmudne i skomplikowane. Metody te można jednak z

powodzeniem zastosować do obliczeń numerycznych. W oparciu o nie został

stworzony zbiór funkcji pracujących w środowisku MATLAB. Za ich pomocą

wyliczono przedstawione dalej rozkłady prawdopodobieństw i wykreślone na ich

podstawie charakterystyki detekcji dla bloków systemu radarowego. Schemat

systemu przetwarzania sygnału w radarach z podziałem na bloki funkcjonalne

pokazano na następnej stronie na rysunku 4.1.

Rys. 4.1. Schemat blokowy sytemu przetwarzania danych w systemach radarowych

4.1. ROZKŁADY OBWIEDNI SYGNAŁU

W najprostszym przypadku w systemie nie występuje blok stabilizacji

poziomu fałszywego alarmu (CFAR) oraz integrator paczki impulsów. Decyzję o

obecności sygnału w szumie podejmuje się porównując amplitudę chwilową

(obwiednię) sygnału za detektorem fazy z odpowiednio dobranym progiem.

Zgodnie z powyższym schematem obwiednia równa jest pierwiastkowi

kwadratowemu z sumy kwadratów składowych I i Q :

)

(

)

(

)

(

(4.1)

Składowe te można zapisać jak poniżej:

))

(

sin(

)

(

)

(

))

(

cos(

)

(

)

(

(4.2)

gdzie

jest pulsacją dopplerowską (czyli różnicą pulsacji wzorcowej nadajnika

i pulsacji sygnału echa) określającą prędkość radialną obiektu, od którego

pochodzi echo. Przyglądając się postaci składowych I (n) oraz Q (n) widać, że

można je potraktować jako składowe (odpowiednio rzeczywistą i urojoną)

pewnego sygnału analitycznego. Skoro tak, to (wiadomo to z własności sygnału

analitycznego) zmienne losowe I (n) i Q (n) są dla danego n nieskorelowane,

oraz ich rozkłady prawdopodobieństwa są identyczne. Jednakowe są także

widmo mocy i wariancja obu przebiegów [6, 13].

W przypadku gdy prawdziwa jest hipoteza H

, to do odbiornika docierają

jedynie zakłócenia o zerowej wartości średniej. Pierszym etapem obliczeń jest

znalezienie łącznego rozkładu zmiennych losowych obu składowych I (n) i

Q (n) tzn.

( Q

Jeśli zmienne losowe I i Q są niezależne, to ich łączny

rozkład równy jest iloczynowi rozkładów obu zmiennych losowych:

)

(

)

(

)

(

(4.3)

Przykładowo gdy zakłócenia mają rozkład normalny, to łączny rozkład można

zapisać jak poniżej:

)

exp(

)

exp(

)

(

πσ

−

(4.4)

gdzie

jest wariancją szumu za filtrem dopasowanym. Niech zmienna losowa

Z i dodatkowa zmienna

będą wynikiem działania pewnego przekształcenia na

zmiennych losowych I i Q :

































(4.5)

gdzie funkcja g:R

→

określona jest poniższym wzorem (4.6):

































)

(

arctan

(4.6)

Przekształcenie odwrotne g

-1

i jej pochodna mają postać:













−

































−

)

cos(

)

sin(

)

sin(

)

cos(

)

(

)

sin(

)

cos(

(4.7)

natomiast modół jakobianu przekształcenia g

-1

równy jest

(4.8)

Znając go można określić łączny rozkład gęstości prawdopodobieństwa

zmiennych Z oraz

wyrażony za pomocą zmiennych wejściowych I i Q :

⋅

))

sin(

(

))

cos(

(

)

sin(

cos(

(

)

(

(4.9)

Rozkład brzegowy zmiennej losowej Z uzyska się całkując rozkład p

, y

)

względem zmiennej y

∫

)

(

)

(

(4.10)

W przypadku, gdy prawdziwa jest hipoteza H

, to odbierany sygnał jest

sumą szumu o zerowej średniej (takiego jak w przypadku prawdziwości

hipotezy H

) i sygnału użytecznego. Niestety nie można postąpić identycznie,

jak dla hipotezy H

, zmieniając jedynie wartość średnią obu zmiennych I i Q

na równą energii lub amplitudzie sygnału użytecznego w zależności od tego,

czy w systemie jest filtr dopasowany, czy nie. Aby wyznaczyć rozkład obwiedni

p(Z), należy potraktować zakłócenia jako szum wąskopasmowy o

częstotoliwości centralnej

. Można go wtedy zapisać poniższym wyrażeniem:

)

sin(

)

(

)

cos(

)

(

))

(

cos(

)

(

−

(4.11)

gdzie

)

(

)

(

)

(

)

(

)

(

)

(

(4.12)

)

(

)

(

)

(

Jeśli przyjąć początek obserwacji przebiegu wejściowego tak, aby faza sygnału

użytecznego była zerowa, to suma zakłócenia i sygnału zgodnie z wyrażeniami

(4.11) i (4.12) równa jest

)

sin(

)

(

)

cos(

)

(

)

cos(

))

(

cos(

)

(

−

(4.13)

gdzie A jest energią lub amplitudą sygnału użytecznego. Porównując stronami

to wyrażenie wychodzi, że

)

(

))

(

sin(

)

(

)

(

)

(

))

(

cos(

)

(

)

(

−

(4.14)

czyli obwiednia sumy zakłócenia i sygnału będzie równa

)

(

)

(

)

(

)

(

)

(

−

(4.15)

Jak widać, jej rozkład można policzyć korzystając z wyżej przedstawionej

metody wstawiając w miejsce zmiennych losowych I i Q odpowiednio n

i n

z tym, że jedna ze zmiennych (jeśli ich rozkłady są takie same, to nie ma

znaczenia która) ma wartość średnią równą energii/amplitudzie sygnału echa.

Jeśli zmienne losowe n

i n

mają rozkład Gaussa, to ich łączny rozkład można

zapisać poniższą zależnością:

)

(

exp(

)

exp(

)

(

πσ

−

(4.16)

Należy zauważyć, że uzyskane rozkłady dla przypadku, gdy prawdziwa

jest hipoteza H

są tak naprawdę rozkładami warunkowymi, gdyż ich postać

zależy od rozkładu fazy sygnału użytecznego

(ogólnie od wszystkich

losowych parametrów rozkładu). Aby uniezależnić się od fazy należy uśrednić

rozkłady po wszystkich wartościach

. Niestety w tym momencie niezbędna jest

znajomość rozkładu prawdopodobieństwa fazy p(

). Znając go, aby uzyskać

rozkład bezwarunkowy należy obliczyć poniższą całkę:

∫

⋅

)

(

)

(

)

(

(4.17)

Jeśli przyjąć najmniej korzystny, czyli równomierny rozkład fazy, to

)

(

)

(

)

(

)

(

)

(

∫

⋅

(4.18)

i wtedy postać rozkładu zmiennej Z nie ulega zmianie.

Opierając się na powyższej metodzie napisano funkcję działającą w

środowisku MATLAB, obliczającą rozkład prawdopodobieństwa obwiedni

sygnału za detektorem fazy, czyli rozkład zmiennej losowej równej

pierwiastkowi kwadratowemu z sumy kwadratów dwóch zmiennych losowych X

i Y. Prototyp funkcji zapisany jest poniżej:

[p]=sqrtpdf(x, 'fun1', 'fun2', arg1_fun1, arg1_fun2, arg2_fun1, ...)

Zmienna

jest wektorem wartości zmiennych losowych X i Y, zaś

fun1

fun2

nazwami funkcji generujących rozkłady tych dwóch zmiennych:.p

(x) i p

(x).

Dalej podawane są dodatkowe argumenty dla

fun1

fun2

jak np. średnia i

wariancja. Przykładowe poniższe wywołanie wyliczy gęstość rozkładu

prawdopodowbieństwa Rice’a o parametrach

=1 oraz A=10:

x=[0:0.125:30]’;

[p]=sqrtpdf(x,'normpdf','normpdf',0,1,10,1);

Funkcja

normpdf

jest częścią pakietu statystycznego i wylicza wartości

rozkładu normalnego zgodnie z ogólnie znanym wzorem. Na wykresie 4.2

pokazano przykładowe krzywe rozkładów Rice’a/Rayleya uzyskane funkcją

sqrtpdf

dla różnych wartości parametru A. Jeśli w układzie jest filtr

dopasowany i przyjąć moc szumu za nim

równą 1, to wartości A równe są

napięciowemu stosunkowi Sygnał/Szum sygnału odebranego przez antenę

radaru. Na kolejnym wykresie 4.3 pokazano błąd bezwzględny wartości

rozkładów uzyskanych funkcją

sqrtpdf

i bezpośrednio ze wzoru na rozkład

Rice’a.

0.1

0.2

0.3

0.4

0.5

0.6

p(x

)

S/N=0,0
S/N=0,5
S/N=1,0
S/N=1,5
S/N=2,0
S/N=4,0
S/N=5,0

Rys. 4.2. Krzywe rozkładów prawdopodobieństwa uzyskane funkcją

sqrtpdf

dla wybranych

wartości stosunku Sygnał/Szum. Dla S/N=0 rozkład Rayleya

-40

-35

-30

-25

-20

-15

-10

)-p

1(x

S/N=0,0
S/N=0,5
S/N=1,0
S/N=1,5
S/N=2,0
S/N=4,0
S/N=5,0

Rys. 4.3. Błąd bezwzględny wartości rozkładów uzyskanych funkcją

sqrtpdf

i bezpośrednio

ze wzoru na rozkład Rice’a

Błąd wynika głównie z tego, że aby uzyskać szukany rozkład należy wykonać

całkowanie rozkładu dwuwymiarowego, która to operacja wprowadza błędy w

obliczeniach. Wszystkie całki były liczone z zastosowaniem metody trapezów

standardową funkcją MATLABA

trapz

[7, 14].

Mając rozkłady prawdopodobieństwa dla różnych wartości S/N można

wykreślić charakterystyki detekcji. Wykonuje to funkcja

detect

, której sposób

wywołanie pokazany jest poniżej:

[d, sn]=detect(x, p, pfalse, snmax)

Jako argumenty przyjmuje:

, czyli wektor wartości zmiennej losowej, macierz

, w której każda kolumna zawiera wartości krzywych rozkładu dla różnych

wartości S/N, prawdopodobieństwo fałszywego alarmu (

pfalse

), oraz

maksymalną wartość S/N (

snmax

). Funkcja najpierw wylicza “odwrotne”

dystrybuanty rozkładów, czyli prawdopodobieństwa, że sygnał przekroczy

pewną wartość progową, a następnie wylicza i zwraca na wyjście wektor

wartości S/N (równomiernie rozłożone wartości od zera do

snmax

) oraz wektor

prawdopodobieństw detekcji dla założonej wartości

pfalse

. Na kolejnych

wykresach 4.4a,b pokazano uzyskane tą funkcją przykładowe krzywe detekcji

dla różnych wartości F w skali liniowej (a) i logarytmicznej (b).

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-2
F=1e-4
F=1e-6
F=1e-8
F=1e-10

Rys. 4.4a. Charakterystyki detekcji pojedynczego impulsu dla różnych wartości F (skala liniowa)

-10

-8

-6

-4

-2

S/N

F=1e-2
F=1e-4
F=1e-6
F=1e-8
F=1e-10

Rys. 4.4b. Charakterystyki detekcji pojedynczego impulsu dla różnych wartości F (skala
logarytmiczna)

4.2. ROZKŁADY SYGNAŁU ZA UKŁADEM CFAR

Obliczenie rozkładów prawdopodobieństw sygnału i charakterystyk

detekcji za układem stabilizacji fałszywego alarmu jest dość trudne. Niezależnie

od tego, którą z implementacji pokazanych w poprzednim rozdziale wybrać,

obliczenia polegają na znalezieniu najpierw funkcji gęstości

prawdopodobieństwa średniej z N próbek S

(czyli S

), a następnie rozkład

zmiennej losowej S

wyj

(gdzie S

ma rozkład Rayleya dla szumu lub

rozkład Rice’a dla sumy szumu i sygnału).

Pierwszym etapem jest znalezienie rozkładu zmiennej losowej S

, która z

dokładnością do stałych współczynników jest sumą N zmiennych losowych S

Rozkład sumy zmiennych losowych równy jest splotowi rozkładów sumowanych

zmiennych. Aby szybko go znaleźć dla dużych N, to zamiast liczyć splot można

znaleźć funkcję charakterystyczną rozkładu p(S

), która jest jego transformatą

Fouriera, podnieść ją do N-tej potęgi, a następnie obliczyć odwrotną

transformatę [1, 8, 12]. Po odpowiednim przeskalowaniu uzyska się rozkład

p(S

). Do obliczenia funkcji charakterystycznych metodami numerycznymi

można wykorzystać szybki algorytm FFT. Napisana dla MATLABA funkcja

sumpdf

w zależności od sposobu wywołania tzn.

[py y]=sumpdf(x, px, N, ’sum’)

[py y]=sumpdf(x, px, N, ’mean’)

oblicza albo rozkład sumy albo rozkład estymatora wartości średniej. Na

wykresie 4.5 pokazano przykładowe krzywe rozkładów wartości średniej z sumy

rozkładów Rayleya dla różnych N.

Znając rozkład średniej można określić rozkład zmiennej losowej S

wyj

Najpierw należy znaleźć łączny rozkład zmiennych S

i S

. Ponieważ są one

niezależne (bo w liczeniu wartości S

nie bierze udział komórka testowa), to

łączny rozkład równy jest

)

(

)

(

)

(

Ssr

⋅

(4.19)

Wprowadzona zostanie nowa zmienna losowa

i wtedy można zapisać, że

































wyj

(4.20)

0.5

1.5

2.5

0.5

1.5

2.5

3.5

4.5

N=2
N=5
N=8
N=16
N=64

Rys. 4.5. Rozkłady prawdopodobieństw estymatorów wartości średnich dla rozkładu Rayleya o
parametrze

gdzie funkcja g:R

→

określona jest wzorem

































(4.21)

Przekształcenie odwrotne g

-1

i jej pochodna mają postać jak poniżej [8]:























 ⋅





















−

)

(

(4.22)

W związku z tym moduł jakobianu przekształcenia g

-1

równy jest

(4.23)

i ostatecznie łączny rozkład gęstości prawdopodobieństwa zmiennych losowych

wyj

i S

można zapisać:

Ssr

Swyj

⋅

)

(

)

(

)

(

)

(

(4.24)

Rozkład brzegowy zmiennej losowej S

wyj

uzyska się całkując gęstość

Swyj,Ssr

) względem y

∫

)

(

)

(

Ssr

Swyj

(4.25)

Opierając się na powyższej metodzie napisano funkcję

normaliz

działającą w środowisku MATLAB, obliczającą rozkład prawdopodobieństwa

ilorazu zmiennych losowych, czyli też rozkład zmiennej losowej S

wyj

. Prototyp

funkcji zapisany jest poniżej:

[p]=normaliz(x, m, ’fun’, arg1_fun, arg2_fun, ...)

Zmienna

jest wektorem wartości obu dzielonych przez siebie zmiennych

losowych,

wektorem wartości rozkładu zmiennej z mianownika (np. uzyskany

za pomocą funkcji

sumpdf

), zaś

fun

nazwą funkcji generującej rozkład

zmiennej z licznika. Dalej podawane są dodatkowe argumenty dla

fun

Uzyskawszy rozkłady zmiennych losowych S

wyj

dla różnych wartości

stosunku Sygnał/Szum można za pomocą funkcji

detect

wyliczyć

charakterystyki wykrywania dla różnych wartości F. Przykładowe krzywe

pokazane są na czterech kolejnych wykresach 4.6a,b,c,d. W tabeli 4.1 poniżej

zestawiono straty, jakie wprowadza układ stabilizacji fałszywego alarmu dla

wybranych wartości F i różnych ilości N referencyjnych komórek

odległościowych dla D=0,9.

Tab. 4.1. Straty w dB wprowadzane przez układ CFAR dla D=0,9. Puste miejsce oznacza, że
dana krzywa w zakresie S/N=[0;20] nie osiągnęła wartości 0,9

-2

3,198

1,173

0,718

0,350

0,085

-4

−

2,506

1,497

0,719

0,177

-6

−

3,956

2,330

1,100

0,265

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-2

N=0
N=2
N=5
N=8
N=16
N=64

Rys. 4.6a. Charakterystyki detekcji pojedynczego impulsu za układem CFAR dla F=10

-2

różnych wartości N. N=0 oznacza brak układu CFAR

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-4

N=0
N=2
N=5
N=8
N=16
N=64

Rys. 4.6b. Charakterystyki detekcji pojedynczego impulsu za układem CFAR dla F=10

-4

różnych wartości N. N=0 oznacza brak układu CFAR

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-6

N=0
N=2
N=5
N=8
N=16
N=64

Rys. 4.6c. Charakterystyki detekcji pojedynczego impulsu za układem CFAR dla F=10

-6

różnych wartości N (skala liniowa). N=0 oznacza brak układu CFAR

-6

-5

-4

-3

-2

-1

S/N

F=1e-6

N=0
N=2
N=5
N=8
N=16
N=64

Rys. 4.6d. Charakterystyki detekcji pojedynczego impulsu za układem CFAR dla F=10

-6

różnych wartości N (skala logarytmiczna). N=0 oznacza brak układu CFAR

Jak widać z wartości w tabeli im większa ilość N komórek referencyjnych brana

jest do liczenia średniej, tym straty CFAR są mniejsze i charakterystyka detekcji

bliższa jest tej , gdy nie ma układu CFAR. Im większe N, tym rozkład zmiennej

losowej S

bardziej jest skupiony wokół swojej wartości centralnej (patrz wykres

4.5). Dla dużych wartości N dzielenie zmiennych losowych S

sprowadza się

praktycznie do przeskalowania zmiennej S

przez stałą wartość, a ta operacja

nie zmieni kształtu rozkładu zmiennej S

wyj

w stosunku do S

. Co za tym idzie

charakterystyki detekcji będą takie, jak przy braku układu CFAR. Niestety nie

można brać zbyt dużego N, ponieważ równałoby się to przyjęciu błędnego

założenia, że poziom clutteru jest stały w dużym obszarze przestrzeni, a nie

tylko lokalnie wokół testowanej komórki odległościowej. Dobrą praktyczną

zasadą jest użycie takiej liczby komórek odniesienia, aby utrzymać straty CFAR

poniżej 1dB [11]. Przyglądając się tabeli 4.1 widać, że dla prawdopodobieństwa

fałszywego alarmu F=10

-2

wystarczy już 8 komórek referencyjnych, natomiast

dla F=10

-4

i F=10

-6

wystarczy N równe około 16.

4.3. ROZKŁADY SYGNAŁU ZA INTEGRATOREM BINARNYM

Aby określić prawdopodobieństwo detekcji oraz fałszywego alarmu

integratora binarnego lub rozkłady za integratorem, trzeba znać odpowiednie

wartości i rozkłady dla pojedynczego impulsu. Integrator binarny integruje

paczki po N impulsów, zaś wykryje tylko te paczki, dla których liczba impulsów k

przekraczających próg będzie równa lub większa od założonej przez

obserwatora liczby M. Dla każdego k możliwych jest

(

−









(4.26)

takich paczek. Jeśli prawdopodobieństwo przekroczenia

wartości progowej przez dowolny impuls oznaczyć przez

, zaś prawdopodobieństwo fałszywego przekroczenia

(alarmu) przez F

, to prawdopodobieństwo, że dokładnie k

impulsów z N przekroczy próg (mówiąc ogólnie będzie k

sukcesów w N niezależnych próbach) dane jest wzorem na

prawdopodobieństwo w tzw. schemacie Bernouliego [12]:

−









)

(

)

(

(4.27)

Liczba D

jest prawdopodobieństwem detekcji pojedynczego impulsu, czyli

sukcesu w pojedynczej próbie, a (1-D

) prawdopodobieństwem przeoczenia

impulsu, czyli inaczej mówiąc porażki w pojedynczej próbie. W związku z

powyższym prawdopodobieństwo detekcji (przekroczenia progu) przez M i

więcej impulsów wynosi:

∑

−













−









)

(

(4.28)

Dla M=0 zawsze D=1, ponieważ ze wzoru (4.28) wynika, że

[

]

)

(

−

Analogicznie całe prawdopodobieństwo fałszywego

alarmu to znaczy prawdopodobieństwo, że M i więcej

zakłóceń przekroczy próg) wynosi:

∑

−













−









)

(

(4.29)

jednak musi być M>0. W przeciwnym przypadku F wyniesie 1. Korzystając z

powyższego wzoru i znając prawdopodobieństwa zdarzeń polegających na tym,

że sygnał przekroczy pewną wartość progową X

tzn. P(x>X

) dla pojedynczego

impulsu, można także określić prawdopodobieństwo, że M i więcej impulsów

przekroczy próg:

∑

−













−









))

(

)

(

)

(

(4.30)

Wartości D

i F

występujące w poprzednich wyrażeniach można znaleźć z

krzywych detekcji dla pojedynczego impulsu z przypadkową fazą początkową

(wykresy 4.4a, 4.4b). Przy czym jeżeli energia całej paczki wynosi E, to w celu

wyznaczenia D

z charakterystyk detekcji należy wziąć energię E/N.

Mając rozkłady prawdopodobieństw pojedynczego impulsu dla różnych

wartości S/N można wykreślić charakterystyki detekcji. Wykonuje to funkcja

binint

, której sposób wywołanie pokazany jest poniżej:

[d, sn]=binint(x, p, N, M, pfalse, snmax)

Funkcja jest oparta na wcześniej opisanej funkcji

detect

. Jako argumenty

przyjmuje:

, czyli wektor wartości zmiennej losowej, macierz

, w której każda

kolumna zawiera wartości krzywych rozkładu dla różnych wartości S/N,

parametry integratora binarnego (

), prawdopodobieństwo fałszywego

alarmu (

pfalse

), oraz maksymalną wartość S/N (

snmax

). Funkcja najpierw

wylicza odwrotne dystrybuanty rozkładów, czyli prawdopodobieństwa, że sygnał

przekroczy pewną wartość progową, potem każdą wartość

prawdopodobieństwa modyfikuje w oparciu o wzór na schemat Bernouliego, a

następnie wylicza i zwraca na wyjście wektor wartości S/N (równomiernie

rozłożone wartości od zera do

snmax

) oraz wektor prawdopodobieństw detekcji

dla założonej wartości

pfalse

. Dalej na kolejnych wykresach 4.7a,b,c,d

pokazano uzyskane tą funkcją przykładowe krzywe detekcji dla różnych

wartości F, N oraz M. W tabelach 4.2a,b dla przykładu zestawiono zysk z

integracji dla integratora binarnego w stosunku do detekcji pojedynczego

impulsu dla D=0,9. Pogrubieniem zaznaczono najlepszy integrator na tym

poziomie prawdopodobieństwa detekcji. Na wykresie 4.8 dla lepszego

uwidocznienia została pokazana zawartość tabel 4.2a,b.

Tab. 4.2a. Zysk w dB z zastosowania integratora binarnego dla N=10, D=0,9

Zysk w dB dla N=10, D=0,9

F=10

-6

F=10

-4

1.321

1.513

2.474

2.556

2.986

2.992

3.236

3.192

3.362

3.261

3.390

3.254

3.346

3.178

3.204

3.004

2.934

2.686

2.331

2.020

Tab. 4.2b. Zysk w dB z zastosowania integratora binarnego dla N=5, D=0,9

Zysk w dB dla N=5, D=0,9

F=10

-6

F=10

-4

1.000

1.141

1.965

1.990

2.273

2.213

2.256

2.124

1.835

1.616

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-6, N=10

N=1, M=1
M=1
M=2
M=3
M=4
M=5
M=6
M=7
M=8
M=9
M=10

Rys. 4.7a. Charakterystyki detekcji integratora binarnego dla F=10

-6

i N=10. Krzywa N=1,M=1

jest charakterystyką dla pojedynczego impulsu

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-4, N=10

N=1, M=1
M=1
M=2
M=3
M=4
M=5
M=6
M=7
M=8
M=9
M=10

Rys. 4.7b. Charakterystyki detekcji integratora binarnego dla F=10

-4

i N=10. Krzywa N=1,M=1

jest charakterystyką dla pojedynczego impulsu

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-6, N=5

N=1, M=1
M=1
M=2
M=3
M=4
M=5

Rys. 4.7c. Charakterystyki detekcji integratora binarnego dla F=10

-6

i N=5. Krzywa N=1,M=1

jest charakterystyką dla pojedynczego impulsu

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-4, N=5

N=1, M=1
M=1
M=2
M=3
M=4
M=5

Rys. 4.7d. Charakterystyki detekcji integratora binarnego dla F=10

-4

i N=5. Krzywa N=1,M=1

jest charakterystyką dla pojedynczego impulsu

1.5

2.5

3.5

Zysk w dB z integracji w stosunku do detekcji jednego impulsu dla D=0,9

F=1e-6, N=10
F=1e-4, N=10
F=1e-6, N=5
F=1e-4, N=5

Rys. 4.8. Zysk w dB z zastosowania integratora binarnego w stosunku do detekcji
pojedynczego impulsu dla D=0,9 i różnych wartości F, N

4.4. ROZKŁADY SYGNAŁU ZA INTEGRATOREM DZIAŁAJĄCYM

BEZPOŚREDNIO NA AMPLITUDZIE

Analiza integratora działającego bezpośrednio na amplitudzie sygnału i

sumującego impulsy będzie podzielona na dwie części, gdyż w zależności od

poziomu sygnału można mówić o dwóch różnych urządzeniach.

Zgodnie z tym, co powiedziano w punkcie 3.3 jeśli amplituda impulsu jest

duża w porównaniu szumami, integrator sumuje impulsy zgodnie z poniższym

wyrażeniem:

∑

Jeśli amplituda impulsu jest skrajnie mała lepiej integrator liniowy zastąpić

kwadratowym:

∑

W obu wyrażeniach

jest odchyleniem standardowym szumu na wejściu

integratora, zaś S

współczynnikami wagowymi mającymi zniwelować

modulację paczki przez charakterystykę anteny radaru.

W przypadku detektora liniowego sytuacja jest dość prosta. Aby znaleść

rozkład prawdopodobieństwa sumy zmiennych losowych o takich samych

rozkładach, który równy jest splotowi rozkładów sumowanych zmiennych,

należy postąpić podobnie, jak przy analizie układu CFAR, gdy liczyło się rozkład

średniej. Zamiast wykonywać wiele operacji splotów, należy skorzystać z

przekształcenia rozkładu na funkcję charakterystyczną. Takie sumowanie

wykonuje funkcja

sumpdf

Aby znaleźć rozkład sygnału za integratorem kwadratowym najpierw

należy obliczyć rozkład zmiennej losowej Y=X

, gdzie X o rozkładzie p

(x) jest

sygnałem przed integratorem. Skoro przekształcenie g:R

→

R ma postać

( )

(4.31)

to przekształcenie odwrotne i jej pochodna mogą być zapisane jak poniżej [5]:

( )

)

(

−

(4.32)

Rozkład prawdopodobieństwa zmiennej losowej Y równy jest:

)

(

)

(

⋅

(4.33)

Można go uzyskać np. za pomocą funkcji

sqrtpdf

podstawiając zamiast x

pierwiastek kwadratowy z x i wynik funkcji mnożąc jak w wyrażeniu powyżej.

Mając rozkład zmiennej losowej Y=X

należy najpierw wektor x podzielić, a

wektor krzywej rozkładu pomnożyć przez 4, a następnie metodami takimi, jak

dla integratora liniowego znaleźć rozkład sumy zmiennych losowych.

Znając rozkłady dla różnych wartości S/N można za pomocą funkcji

detect

stworzyć charakterystyki wykrywania. Na kolejnych wykresach 4.9a,b

znajdują się charakterystyki detekcji dla integratora liniowego i kwadratowego,

oraz dla porównania dla najlepszego integratora binarnego. Uzyskano je za

pomocą wyżej omówionych funkcji. Dalej w tabeli 4.3 zestawiono zysk dla

integratorów w stosunku do detekcji pojedynczego impulsu, dla D=0,9.

Tab. 4.3. Zysk w dB z zastosowania różnych typów integratorów dla D=0,9 i różnych M w
stosunku do pojedynczego impulsu

Zysk w dB, D=0,9

Typ

INT.

F=10

-6

F=10

-4

Lin

−

2.929

2.812

Lin

−

4.035

3.875

Sqr

−

2.840

2.725

Sqr

−

3.952

3.800

Bin

2.273

2.213

Bin

3.390

−

Bin

−

3.261

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-6

N=1, M=1
lin: N=5
lin: N=10
sqr: N=5
sqr: N=10
bin: N=5, M=3
bin: N=10, M=6

Rys. 4.9a. Charakterystyki detekcji dla różnych integratorów dla F=10

-6.

Krzywa N=1,M=1 jest

charakterystyką dla pojedynczego impulsu. Wartości S/N są dla pojedynczego impulsu

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

F=1e-4

S/N

N=1, M=1
lin: N=5
lin: N=10
sqr: N=5
sqr: N=10
bin: N=5, M=3
bin: N=10, M=6

Rys. 4.9b. Charakterystyki detekcji dla różnych integratorów dla F=10

-4.

Krzywa N=1,M=1 jest

charakterystyką dla pojedynczego impulsu. Wartości S/N są dla pojedynczego impulsu

W rozdziale trzecim (punkt 3.3) wspomniano, iż w krańcowej sytuacji, gdy

poziomy sygnałów (impulsów) są bardzo małe w porównaniu z zakłóceniami,

detektor kwadratowy będzie lepszy od liniowego. Tak przynajmniej wynika z

kształtu krzywej

)

(

. Jak widać na przedstawionych wykresach 4.9a i 4.9b

integrator liniowy nie jest gorszy od kwadratowego. Jest nawet w pokazanych

przypadkach w niewielkim stopniu lepszy. To, jaki typ integratora korzystniej

zastosować zależy od energii pojedynczego impulsu w stosunku do mocy

szumu, żądanego poziomu prawdopodobieństwa fałszywego alarmu czy liczby

integrowanych impulsów. Generalnie dla małych prawdopodobieństw

fałszywego alarmu i dużych prawdopodobieństw detekcji poziom sygnału

progowego dla obu typów integratorów jest taki sam [5]. Na wykresie 4.10

pokazano skrajny przypadek, gdzie integrator kwadratowy jest nieco lepszy od

liniowego. Zysk jest tak minimalny, że jego stosowanie zamiast liniowego i

niepotrzebne komplikowanie układu nie ma sensu.

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

S/N

F=1e-1

lin: N=10

sqr: N=10

Rys. 4.10. Charakterystyki detekcji dla integratorów liniowego i kwadratowego dla F=10

-1

i N=10

PODSUMOWANIE

W pracy dokonano analizy statystycznej różnego rodzaju detektorów

cyfrowych, stosowanych we współczesnych systemach radarowych.

Przedstawiono podstawy teoretyczne procesu optymalnej detekcji oraz

rozważania prowadzące od teoretycznego optymalnego detektora do jak

najbliższych mu praktycznych realizacji. Dla rzeczywistych rozwiązań opisano

metody pozwalające określić rozkłady prawdopodobieństw sygnałów w

przypadkach odbierania jedynie zakłóceń oraz gdy odbiera się sumę zakłóceń i

sygnału użytecznego.

W oparciu o te metody napisano zestaw funkcji pracujących w środowisku

MATLAB, pozwalających na numeryczne wyznaczenie rozkładów

prawdopodobieństw, a także na wykreślenie na ich podstawie charakterystyk

detekcji dla różnych typów detektorów. Zdecydowano się na to środowisko z

uwagi na jego uniwersalność, a przede wszystkim na łatwość implementacji

różnorodnych algorytmów obliczeń numerycznych. MATLAB, jako środowisko i

język programowania wysokiego poziomu, idealnie nadaje się do tego typu

zadań obliczeniowych. Co najważniejsze, środowisko to jest otwarte i pozwala

na integrację własnych procedur z już istniejącymi.

W dodatkach A i B znajdują się opis stworzonych funkcji, służących do

analizy detektorów cyfrowych oraz wydruk ich kodu źródłowego. Z ich

wykorzystaniem dokonano przedstawionej w rozdziale czwartym analizy

wybranych detektorów w przypadku, gdy zakłócenia mają postać szumu białego

o rozkładzie Gauss’a oraz dzięki możliwościom graficznym MATLABA

zobrazowano wyniki na wykresach.

Uzyskane wykresy i wyniki świadczą o przydatności opracowanego

zestawu oprogramowania do analizy detektorów cyfrowych. Napisane funkcje

pracujące w środowisku MATLAB mogą służyć jako baza przy tworzeniu

bardziej rozbudowanych i wyspecjalizowanych aplikacji.

100

DODATEK A:

OPIS ZESTAWU FUNKCJI SŁUŻĄCYCH DO

ANALIZY DETEKTORÓW CYFROWYCH

Niniejszy dodatek stanowi dokumnetację zestawu funkcji służących do analizy detektorów

cyfrowych. Funkcje te zostały napisane i przetestowane w środowisku MATLAB w wersji 5.2. Niektóre z
nich wewnątrz siebie wywołują zarówno inne funkcje z pakietu, jak i te należące do standardowej
dystrybucji MATLAB’a. Z tego powodu istnieje możliwość, że w starszych wersjach środowiska
zadziałają niepoprawnie. W skład pakietu wchodzą następujące funkcje:

binint

wylicza charakterystyki detekcji dla integratora binarnego

dbloss

liczy różnicę w decybelach między dwoma charakterystykami detekcji dla określonego
poziomu prawdopodobieństwa detekcji

detect

wylicza charakterystyki detekcji dla pojedynczego impulsu

normaliz

wylicza rozkład prawdopodobieństwa zmiennej losowej Z=X/Y

pdftocdf

wylicza dystrybuantę zmiennej losowej P(x) lub krzywą 1-P(x)

ricepdf

wylicza krzywą rozkładu prawdopodobieństwa Rice’a

sqrtpdf

wylicza rozkład prawdopodobieństwa zmiennej losowej

sumpdf

generuje rozkład prawdopodobieństwa zmiennej losowej będącej

sumą N niezależnych zmiennych losowych o takich samych

rozkładach.

Na kolejnych stronach znajdują się dokładne opisy wyżej wymienionych funkcji,

natomiast dodatek B zawiera ich kod źródłowy.

binint

SKŁADNIA:

[d, sn]=binint(x, p, N, M, pfalse, snmax)

OPIS:
Funkcja generuje charakterystyki detekcji dla układu integratora binarnego (zwanego też dwuprogowym
lub M z N) przy zadanym poziomie prawdopodobieństwa fałszywego alarmu. Na podstawie rozkładów
prawdopodobieństw sygnału dla różnych wartości stosunków SYGNAŁ/SZUM

)

( E

gdzie E jest energią sygnału użytecznego najpierw oblicza “odwrotną”

dystrybuantę

)

(

−

101

czyli prawdopodobieństwo, że sygnał przekroczy pewien poziom x. Robi to wywoływana wewnątrz

binint

funkcja

pdftocdf

(jej opis znajduje się kilka stron dalej). Następnie zgodnie ze wzorem na

schemat Bernouliego liczy prawdopodobieństwo, że M i więcej impulsów z N w paczce przekroczy
poziom x. Z tak uzyskanych danych dla zadanego poziomu prawdopodobieństwa fałszywego alarmu F
wyliczane są krzywe detekcji. Ponieważ rozkłady prawdopdobieństw mają postać dyskretną (zarówno
zmienna losowa X, jak i odpowiadające jej wartości prawdopodobieństw zmieniają się w sposób
skokowy), to aby określić prawdopodobieństwo dla x

o wartości z pomiędzy danych x-ów tzn.

zastosowano metodę aproksymacji liniowej.

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy zmiennej losowej X (inaczej wektor wartości,

jakie może przyjmować sygnał i szum w pojedynczym impulsie)

−

macierz, w której każda kolumna jest rozkładem

prawdopodobieństwa zmiennej losowej X dla pojedynczego

impulsu, przy czym kolumna pierwsza jest rozkładem tylko dla

szumu, zaś kolejne kolumny rozkładami sumy sygnał+szum dla

różnych coraz większych wartości stosunku SYGNAŁ/SZUM

−

parametry integratora

pfalse

−

ustalone dopuszczalne prawdopodobieństwo fałszywego

alarmu

snmax

−

wartość maksymalnego stosunku SYGNAŁ/SZUM (dla rozkładu

z ostatniej kolumny macierzy

)

WYJŚCIE:

−

wektor kolumnowy stosunku SYGNAŁ/SZUM (wektor

równomiernie rozłożonych wartości od zera do

snmax

−

wektor kolumnowy prawdopodobieństw detekcji

PRZYKŁAD:

Poniższa sekwencja wyliczy charakterystykę detekcji dla integratora binarnego

typu 5 z 10 przy założeniu, że pojedynczy impuls ma rozkład Rice’a,

maksymalny stosunek SYGNAŁ/SZUM wynosi 20, a prawdopodobieństwo

fałszywego alarmu F=10

-6

. Funkcja

ricepdf

nie jest częścią Matlaba. Jej opis

znajduje się kilka stron dalej.

E=[0:0.5:20]’; x=[0:0.1:40]’;
for i=1:length(E)

p(:,i)=ricepdf(x, 1, E(i));

end
[d sn]=binint(x, p, 10, 5, 1e-6, 20);

102

dbloss

SKŁADNIA:

d=dbloss(x1, f1x, x2, f2x, l)

OPIS:
Funkcja liczy różnicę w decybelach między dwoma charakterystykami detekcji dla określonego poziomu
prawdopodobieństwa detekcji tzn. różnicę między stosunkami SYGNAŁ/SZUM, dla których obie krzywe
mają taką samą wielkość prawdopodobieństwa detekcji. Funkcja także działa dobrze dla wszystkich
krzywych niemalejących. Jeśli po wywołaniu funkcji pojawi się komunikat o błędzie i niezdefiniowaniu
zmiennych

lub

, oznacza to, że któryś z przebiegów dla zadanego poziomu

nie ma określonej

wartości.

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy dziedziny (wartości stosunków

SYGNAŁ/SZUM) pierwszej krzywej

fx1

−

wektor kolumnowy zbioru wartości (prawdopodobieństw

detekcji) pierwszej krzywej

fx2

−

wektory kolumnowe dziedziny i zbioru wartości drugiej krzywej

−

wartość, dla której liczona jest różnica stosunków

SYGNAŁ/SZUM

WYJŚCIE:

−

różnica w dB między stosunkami SYGNAŁ/SZUM, dla których

obie krzywe mają taką samą wartość prawdopodobieństwa

detekcji

PRZYKŁAD:

Liczona jest różnica dla dwóch różnych wartości

x1=[1 2 3 4 5]';

fx1=[1 2 3 4 5]';

x2=[6 7 8 9 10]';

fx2=[1 3 3 7 10]';

dbloss(x1,fx1,x2,fx2,3)

% ans = 3.67976785294594

dbloss(x1,fx1,x2,fx2,6)

% błąd, bo dla l=6 fx1 nie ma
określonej wartości

103

detect

SKŁADNIA:

[d, sn]=detect(x, p, pfalse, snmax)

OPIS:
Funkcja generuje charakterystyki detekcji dla pojedynczego impulsu przy zadanym poziomie
prawdopodobieństwa fałszywego alarmu. Na podstawie rozkładów prawdopodobieństw sygnału dla
różnych wartości stosunków SYGNAŁ/SZUM

)

( E

gdzie E jest energią sygnału użytecznego, oblicza najpierw “odwrotną”

dystrybuantę

)

(

−

czyli prawdopodobieństwo, że sygnał przekroczy pewien poziom x. Robi to wywoływana wewnątrz

detect

funkcja

pdftocdf

(jej opis znajduje się kilka stron dalej). Z tak uzyskanych danych dla

zadanego poziomu prawdopodobieństwa fałszywego alarmu F wyliczane są krzywe detekcji. Ponieważ
rozkłady prawdopdobieństw mają postać dyskretną (zarówno zmienna losowa X, jak i odpowiadające jej
wartości prawdopodobieństw zmieniają się w sposób skokowy), to aby określić prawdopodobieństwo dla
x

o wartości z pomiędzy danych x-ów tzn.

zastosowano metodę aproksymacji liniowej.

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy zmiennej losowej X (inaczej wektor wartości,

jakie może przyjmować sygnał i szum w pojedynczym impulsie)

−

macierz, w której każda kolumna jest rozkładem

prawdopodobieństwa zmiennej losowej X dla pojedynczego

impulsu, przy czym kolumna pierwsza jest rozkładem tylko dla

szumu, zaś kolejne kolumny rozkładami sumy sygnał+szum dla

różnych coraz większych wartości stosunku SYGNAŁ/SZUM

pfalse

−

ustalone dopuszczalne prawdopodobieństwo fałszywego

alarmu

snmax

−

wartość maksymalnego stosunku SYGNAŁ/SZUM (dla rozkładu

z ostatniej kolumny macierzy

)

WYJŚCIE:

−

wektor kolumnowy stosunku SYGNAŁ/SZUM (wektor

równomiernie rozłożonych wartości od zera do

snmax

−

wektor kolumnowy prawdopodobieństw detekcji

PRZYKŁAD:

104

Poniższa sekwencja wyliczy charakterystykę detekcji dla opisanego w rozdziale

pierwszym przypadku miernika wychyłowego, gdzie mierzona wielkość zmienia

się od zera do 5, błąd pomiaru ma rozkład Gaussa o

=0,2, F=10

-4

X=[0:0.5:5]’;
y=[-10:0.1:10]’;
for i=1:length(X)

p(:,i)=normpdf(x, X(i), 0.2);

end
[d sn]=detect(x, p, 1e-4, 5);

105

normaliz

SKŁADNIA:

[p]=normaliz(x, m, 'fun', varargin)

OPIS:
Funkcja generuje rozkład prawdopodobieństwa zmiennej losowej Z:

gdzie X jest zmienną losową, której rozkład tworzony jest na bazie argumentu

przez dowolną funkcję oznaczoną w składni jako

fun

, zaś wartości

prawdopodobieństw rozkładu zmiennej Y zawiera wektor

Wektor

może

mieć inną długość niż

, ale przyjmuje się, że zmienna Y tak jak i X należy do

przedziełu

max

min

; x

Funkcja określona jako

fun

musi być tak napisana, aby jej pierwszym

argumentem był wektor

. Szczegóły na temat metody zastosowanej do

generacji rozkładu zmiennej Z można znaleźć w rozdziale czwartym.

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy zmiennej losowej X

−

wektor kolumnowy rozkładu zmiennej losowej Y

fun

−

nazwa funkcji generującej rozkład zmiennej X

varargin

−

zamiast tego napisu należy podać wszystkie pozostałe

argumenty (oddzielone przecinkami) dla funkcji

fun

WYJŚCIE:

−

wektor kolumnowy opisujący rozkład zmiennej losowej Z

PRZYKŁAD:

Sekwencja na następnej stronie wyliczy rozkłady prawdopodobieństwa sygnału

za układem CFAR dla N=16. Moc sygnału zmienia się od zera do 20.

x=linspace(0,40,1024)';

m=sumpdf(linspace(0,40,8192)',

raylpdf(linspace(0,40,8192)',1),16);

for A=0:0.25:20

n(:,A)=normaliz(x, m,'ricepdf', 1, A);

106

end

107

pdftocdf

SKŁADNIA:

[c]=pdftocdf(x, p, 'fwd')
[c]=pdftocdf(x, p, 'inv') lub [c]=pdftocdf(x, p)

OPIS:

Funkcja na podstawie rozkładu prawdopodobieństwa generuje dystrybuantę

zmiennej losowej X:

)

lub krzywą

)

(

−

określającą prawdopodobieństwo, że zmienna losowa będzie większa od

wartości x. Z uwagi na to, że zarówno zmienna losowa x, jak i odpowiadające

jej wartości prawdopodobieństw zmieniają się w sposób skokowy dystrybuanta

liczona jest nie poprzez całkowanie (np. metodą trapezów) ale jako suma

kumulatywna wartości wektora p(x) pomnożona przez krok, co jaki zmieniają się

wartości wektora

. Sumę kumulatywną liczy bardzo szybka wbudowana

funkcja MATLAB’a

cumsum

. W celu szybkiego wyliczenia krzywej 1-P(x)

zastosowano przed i po wywołaniu

cumsum

funkcję

flipud

odwracającą

kolejność danych w wektorze.

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy zmiennej losowej X

−

wektor kolumnowy funkcji gęstości prawdopodobieństwa p(x).

Jeśli argument

jest macierzą, to jej kolumny są traktowane

jako rozkłady warunkowe tej samej zmiennej losowej X

‘

fwd

’

−

liczona jest dystrybuanta zmiennej losowej X: P(x)

‘

inv

’

−

liczona jest "odwrotna" dystrybuanta tzn. 1-P(x). Niepodanie

parametru jest równoważne podaniu '

inv

WYJŚCIE:

108

−

"prosta" lub "odwrotna" dystrybuanta zmiennej X (wektor

kolumnowy lub macierz, jeśli argument wejściowy

jest

macierzą)

PRZYKŁAD:

Poniższa sekwencja wyliczy i wykreśli dystrybuantę standardowego rozkładu

normalnego:

x=[0:0.5:5]’;
p=normpdf(x, 0, 1);

c=pdftocdf(x, p,

fwd

);

109

ricepdf

SKŁADNIA:

[p]=ricepdf(x, sigma, A)

OPIS:

Funkcja generuje rozkład prawdopodobieństwa Rice’a o parametrach

oraz A zgodnie z poniższym wzorem:

(

)

exp(

)

(

⋅

−

gdzie

jest zmodyfikowaną funkcją Bessela pierwszego rodzaju zerowego rzędu. W funkcji zmieniono

odpowiednio kolejność obliczeń tak, aby dla dużych x-ów nie dochodziło do liczenia iloczynu

∞

⋅

przypadku, gdy A=0 wzór znacznie się upraszcza:

exp(

)

(

−

Opisuje on funkcję gęstości prawdopodobieństwa Rayleya.

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy zmiennej losowej X

sigma

−

parametry rozkładu:

, A

WYJŚCIE:

−

wektor kolumnowy opisujący rozkład zmiennej losowej X

PRZYKŁAD:
Następująca sekwencja wygeneruje i wykreśli rozkład Rice’a o parametrach

=1 oraz A=5:

x=[0:0.01:10]’;
p=ricepdf(x, 1, 5);
plot(x,p);

110

sqrtpdf

SKŁADNIA:

[p]=sqrtpdf(x, 'fun1', 'fun2', varargin)

[p]=sqrtpdf(x, '', '', mu1, mu2, sigma1, sigma2)

OPIS:

Funkcja generuje rozkład prawdopodobieństwa zmiennej losowej Z:

gdzie X jest zmienną losową, której rozkład tworzony jest na bazie argumentu

przez dowolną funkcję oznaczoną w składni jako

fun1

, zaś Y analogicznie

przez funkcję oznaczoną w składni jako

fun2

. Obie funkcje muszą być tak

napisane, aby ich pierwszym argumentem był wektor

. W przypadku

niepodania nazw funkcji (drugi sposób wywołania) zmienne X i Y będą miały

rozkład Gaussa o parametrach:

mu1

sigma1

mu2

sigma2

Szczegóły na temat metody zastosowanej do generacji rozkładu zmiennej Z

można znaleźć w rozdziale czwartym.

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy zmiennej losowej X i Y

fun1

−

dowolna nazwa istniejącej funkcji generującej rozkład zmiennej

fun2

−

jak wyżej dla rozkładu zmiennej Y

varargin

−

zamiast tego napisu należy podać wszystkie pozostałe

argumenty (oddzielone przecinkami) dla funkcji

fun1

fun2

kolejno: arg1_fun1, arg1_fun2, arg2_fun1, arg2_fun2, ...

WYJŚCIE:

−

wektor kolumnowy wartości prawdopodobieństw dla rozkładu

zmiennej losowej Z (wartości zmiennej Z opisuje wektor

)

PRZYKŁAD:
Poniższa sekwencja wyliczy i wykreśli rozkład Rice’a o parametrach

=1 oraz A=5. Jest równoważna

przykładowi podanemu dla funkcji

ricepdf

x=[0:0.01:10]’;

111

p=sqrtpdf(x, 'normpdf', 'normpdf', 0, 5, 1, 1);
plot(x,p);

112

sumpdf

SKŁADNIA:

[py y]=sumpdf(x, px, N, 'mean') lub [py y]=sumpdf(x, px, N)
[py y]=sumpdf(x, px, N, 'sum')

OPIS:

Funkcja generuje rozkład prawdopodobieństwa zmiennej losowej będącej sumą N
niezależnych zmiennych losowych X o takich samych rozkładach. W zależności od
tego, jaki podano czwarty argument funkcji, generuje ona rozkład średniej
opisany poniższym wzorem

∑

−

gdy podano

'mean'

lub rozkład sumy zmiennych X:

∑

−

dla czwartego argumentu

'sum'

W celu przyspieszenia obliczeń zamiast N-razy liczyć splot zastosowano

przejście przez funkcję charakterystyczną, która równa jest transformacie

Fouriera z rozkładu prawdopodobieństwa. Funkcja

sumpdf

najpierw za

pomocą metody FFT liczy dyskretną transformatę Fouriera o ilości

punktów

razy długość wektora

, następnie podnosi do N-tej potęgi i

liczy odwrotną DTF. W zależności od tego, czy liczony jest rozkład

średniej, czy sumy wykonywane jest odpowiednie skalowanie wyników

przez

PARAMETRY WEJŚCIOWE:

−

wektor kolumnowy zmiennej losowej X

−

wektor kolumnowy funkcji gęstości prawdopodobieństwa p(x)

−

liczba sumowanych zmiennych losowych

‘

mean

’

−

funkcja wygeneruje rozkład średniej

‘

sum

’

−

funkcja wygeneruje rozkład sumy. Parametr ten jest opcjonalny.

Niepodanie go jest równoważne podaniu '

mean

WYJŚCIE:

−

wektor kolumnowy opisujący rozkład zmiennej losowej Y

−

wektor kolumnowy zmiennej losowej Y

113

PRZYKŁAD:
Następująca sekwencja wygeneruje rozkład estymatora wartości średniej dla N=10 pomiarów, gdzie
mierzona wartość przyjmuje wartości od zera do 10:

x=[0:0.01:10]’;
px=ones(length(x),1).*0.1

[py y]=sumpdf(x, px, 10,

mean

);

plot(x,px, y, py);

114

DODATEK B:

KODY ŹRÓDŁOWE FUNKCJI SŁUŻĄCYCH DO

ANALIZY DETEKTORÓW CYFROWYCH

binint.m

function [d, sn]=binint(x, p, N, M, pfalse, snmax)
%
%[d, sn]=binint(x, p, N, M, pfalse, snmax)
%
% Generuje charakterystyke detekcji, czyli wykres prawdopodobienstwa
% detekcji od stosunku SYGNAL/SZUM dla ukladu integratora binarnego
% (integratora typu M z N), przy zadanym poziomie prawdopodobienstwa
% falszywego alarmu
%
% PARAMETRY WEJSCIOWE:
% x - zmienna losowa X (wektor kolumnowy), czyli wektor
% wartosci, jakie moze przyjmowac sygnal i szum w
% pojedynczym impulsie
% p - macierz, w ktorej kazda kolumna jest rozkladem
% prawdopodobienstwa zmiennej losowej X dla pojedynczego
% impulsu, przy czym kolumna pierwsza jest rozkladem
% tylko dla szumu, zas kolejne kolumny rozkladami
% sumy sygnal+szum dla roznych coraz wiekszych wartosci
% stosunku SYGNAL/SZUM
% N, M - parametry integratora binarnego (M z N)
% pfalse - ustalone dopuszczalne prawdopodobienstwo falszywego
% alarmu
% snmax - warotsc maksymalnego stosunku SYGNAL/SZUM (dla
% rozkladu z ostatniej kolumny macierzy 'p')
% WYJSCIE:
% sn - wektor kolumnowy stosunku SYGNAL/SZUM (wektor
% rownomiernie rozlozonych wartosci od zera do 'snmax'
% d - wektor kolumnowy prawdopodobienst detekcji
%
% Sprawdzone w srodowisku MATLAB 5.x
%

%tic
if (N<0 | M<0 | M>N)
display('Blad danych wejsciowych');
return;
end

p=pdftocdf(x, p);

siz=size(p);

115

for i=1:siz(1)
for j=1:siz(2)
p_out=0;
%p_out=sum(binopdf([M:1:N]', N, p(i,j)));
for k=M:N
p_out=p_out+(nchoosek(N,k).*(p(i,j).^k).*((1-p(i,j)).^(N-
k)));

end

p(i,j)=p_out;
%i
end
end

for i=1:siz(1)-1
if p(i,1)>pfalse & p(i+1,1)<pfalse
i1=i;
i2=i+1;
break;
end
end

a=(p(i2,1)-p(i1,1))./(i2-i1);
b=p(i1,1)-a.*i1;
x=(pfalse-b)./a;
d(1,1)=pfalse;

for i=2:siz(2)
a=(p(i2,i)-p(i1,i))./(i2-i1);
b=p(i1,i)-a.*i1;
d(i,1)=a.*x+b;
end

sn=linspace(0,snmax,siz(2))';
%toc

dbloss.m

function d=dbloss(x1, f1x, x2, f2x, l)
%
%d=dbloss(x1, f1x, x2, f2x, l)
%
% Funkcja liczy roznice w decybelach miedzy dwoma charakterystykami
% detekcji dla okreslonego poziomu prawdopodobienstwa detekcji tzn.
% roznice miedzy stosunkami SYGNAL/SZUM, dla ktorych obie krzywe
% maja taka sama wielkosc prawdopodobienstwa detekcji.
%
% PARAMETRY WEJSCIOWE:
% x1 - wektor kolumnowy dziedziny (wartosci stosunkow
% SYGNAL/SZUM) pierwszej krzywej
% f1x - wektor kolumnowy zbioru wartosci (prawdopodobienstw
% detekcji) pierwszej krzywej
% x2, f2x - wektory kolumnowe dziedziny i zbioru wartosci drugiej
% krzywej
% l - wartosc, dla ktorej liczona jest roznica stosunkow
% SYGNAL/SZUM
% WYJSCIE:
% d - roznica w dB miedzy stosunkami SYGNAL/SZUM, dla
% ktorych obie krzywe maja taka sama wartosc
% prawdopodobienstwa detekcji
%
% Jesli funkcja wyswietli

blad, że zmienna 'i1' lub 'i2' jest

116

% niezdefiniowana, to oznacza, ze ktorys z przebiegow
% dla zadanego poziomu 'l' nie ma okreslonej wartosci
%
% Sprawdzone w srodowisku MATLAB 5.x
%

len=length(f1x);
for i=1:len-1
if f1x(i,1)<l & f1x(i+1,1)>=l
i1=i;
i2=i+1;
break;
end
end
a=(f1x(i2,1)-f1x(i1,1))./(i2-i1);
b=f1x(i2,1)-a.*i2;
xx=(l-b)./a;
a=(x1(i2,1)-x1(i1,1))./(i2-i1);
b=x1(i1,1)-a.*i1;
d1=a.*xx+b;

len=length(f2x);
for i=1:len-1
if f2x(i,1)<l & f2x(i+1,1)>=l
i1=i;
i2=i+1;
break;
end
end
a=(f2x(i2,1)-f2x(i1,1))./(i2-i1);
b=f2x(i2,1)-a.*i2;
xx=(l-b)./a;
a=(x2(i2,1)-x2(i1,1))./(i2-i1);
b=x2(i1,1)-a.*i1;
d2=a.*xx+b;

d1, d2
d=abs(10*log10(d2)-10*log10(d1));

detect.m

function [d, sn]=detect(x, p, pfalse, snmax)
%
%[d, sn]=detect(x, p, pfalse, snmax)
%
% Generuje charakterystyke detekcji, czyli wykres prawdopodobienstwa
% detekcji od stosunku SYGNAL/SZUM dla sygnalu ukrytego w szumie,przy
% zadanym poziomie prawdopodobienstwa falszywego alarmu
%
% PARAMETRY WEJSCIOWE:
% x - zmienna losowa X (wektor kolumnowy), czyli wektor
% wartosci, jakie moze przyjmowac sygnal i szum
% p - macierz, w ktorej kazda kolumna jest rozkladem
% prawdopodobienstwa zmiennej losowej X, przy czym
% kolumna pierwsza jest rozkladem tylko dla szumu, zas
% kolejne kolumny rozkladami sumy sygnal+szum dla
% roznych coraz wiekszych wartosci stosunku SYGNAL/SZUM
% pfalse - ustalone dopuszczalne prawdopodobienstwo falszywego
% alarmu
% snmax - warotsc maksymalnego stosunku SYGNAL/SZUM (dla
% rozkladu z ostatniej kolumny macierzy 'p')

117

% WYJSCIE:
% sn - wektor kolumnowy stosunku SYGNAL/SZUM (wektor
% rownomiernie rozlozonych wartosci od zera do 'snmax'
% d - wektor kolumnowy prawdopodobienst detekcji
%
% Sprawdzone w srodowisku MATLAB 5.x
%

tic
p=pdftocdf(x, p);

siz=size(p);
for i=1:siz(1)-1
if p(i,1)>pfalse & p(i+1,1)<pfalse
i1=i;
i2=i+1;
break;
end
end

a=(p(i2,1)-p(i1,1))./(i2-i1);
b=p(i1,1)-a.*i1;
xx=(pfalse-b)./a;
d(1,1)=pfalse;

for i=2:siz(2)
a=(p(i2,i)-p(i1,i))./(i2-i1);
b=p(i1,i)-a.*i1;
d(i,1)=a.*xx+b;
end

sn=linspace(0,snmax,siz(2))';
toc

normaliz.m

function [p]=normaliz(x, m, fun, varargin)
%
%[p]=normaliz(x, m, 'fun', varargin)
%
% Generuje funkcje gestosci prawdopodobienstwa zmiennej losowej
% Z = X ./ Y, gdzie:
% X: zmienna losowa, ktorej rozklad tworzy na bazie parametru
% wejsciowego 'x' wywolywana wewnatrz funkcja, ktorej nazwe
% przekazuje sie przez argument 'fun'
% (X nalezy do przedzialu <xmin..xmax>)
% Y: zmienna losowa o rozkladzie opisanym wektorem 'm'
% (Y jak i X takze nalezy do przedzialu <xmin..xmax>, ale wektor
% 'm' moze miec inna dlugosc niz 'x')
% Zastosowano metode przejscia przez rozklad dwuwymiarowy tzn.
% /
% f(z)=|f(x,y)dy, zas
% /
% f(x,y)=f(x1*y1,y1)*|y1|=f(x1*y1)*f(y1)*|y1|,
% jesli zmienne sa niezalezne. Szczegoly w ksiazce A. Pacuta
% "Prawdopodobienstwo, teoria, modelowanie
% probabilistyczne w technice".
%
% PARAMETRY WEJSCIOWE:
% x - wektor kolumnowy zmiennej losowej X
% m - wektor kolumnowy rozkladu zmiennej losowej Y

118

% fun - nazwa funkcji generujacej rozklad zmiennej X
% varargin - dodatkowe argumenty funcji 'fun'
% (pierwszym parametrem, ktorego sie tu nie podaje jest
% zawsze 'x')
% WYJSCIE:
% p - wektor kolumnowy opisujacy rozklad zmiennej losowej Z
%
% Sprawdzone w srodowisku MATLAB 5.x
%

X=linspace(x(1),x(length(x)),length(m))';
MX=m.*X;
Z=zeros(length(m),1);
p=zeros(length(x),1);

%tic
for i=1:length(x)
Z=feval(fun, x(i,1).*X, varargin{:});
Z=Z.*MX;
p(i,1)=trapz(X,Z);
%i
end
%toc

pdftocdf.m

function [c]=pdftocdf(x, p, todo)
%
%[c]=pdftocdf(x, p, 'fwd')
%[c]=pdftocdf(x, p, 'inv') lub [c]=pdftocdf(x, p)
%
% Generuje dystrybuante zmiennej losowej X: P(x) lub
% "odwrotna" dystrybuante tzn. 1-P(x)
%
% PARAMETRY WEJSCIOWE:
% x - wektor kolumnowy zmiennej losowej X
% p - wektor kolumnowy funkcji gestosci
% prawdopodobienstwa p(x).
% Jesli argument 'p' jest macierza, to jej kolumny sa
% traktowane jako rozklady warunkowe tej samej
% zmiennej losowej X
% 'fwd' - liczona jest dystrybuanta zmiennej losowej X: P(x)
% 'inv' - liczona jest "odwrotna" dystrybuanta tzn. 1-P(x).
% Niepodanie parametru jest rownowazne podaniu 'inv'
% WYJSCIE:
% c - "prosta" lub "odwrotna" dystrybuanta zmiennej X
% (wektor kolumnowy lub macierz, jesli argument
% wejsciowy 'p' jest macierza)
%
% Sprawdzone w srodowisku MATLAB 5.x
%

%tic
if exist('todo')==0 | isempty(todo)
todo='inv';
end

step=abs(x(2)-x(1));

if strcmp(todo, 'inv')

119

p=flipud(p);
c=cumsum(p).*step;
c=flipud(c);
else
c=cumsum(p).*step;
end
%toc

ricepdf.m

function [p]=ricepdf(x, sigma, A)
%
%[p]=ricepdf(x, sigma, A)
%
% Generuje funkcje gestosci prawdopodobienstwa Rice'a
% (lub Rayley'a, jesli A=0) zgodnie ze wzorem:
%
% p(x)=x/sigma^2 * exp(- x^2+A^2/2sigma^2) * I (xA/sigma^2)
% 0
% PARAMETRY WEJSCIOWE:
% x - wektor kolumnowy zmiennej losowej X
% sigma, A - parametry rozkladu
%
% WYJSCIE:
% p - wektor kolumnowy opisujacy rozklad zmiennej losowej X
%
% Sprawdzone w srodowisku MATLAB 5.x
%

%tic
%p=(x./(sigma.^2)).*exp(-
((x.^2)+(A.^2))./(2.*sigma.^2)).*besseli(0,(x.*A)./(sigma.^2));
p=(x./(sigma.^2)).*exp((-(x.^2)+2.*A.*x)./(2.*sigma.^2)).*exp(-
(A.^2)./(2.*sigma.^2));
p=p.*besseli(0,(x.*A)./(sigma.^2),1);
%for i=1:length(p)
% if isnan(p(i)) | isinf(p(i))
% p(i)=0;
% end
%end
%toc

sqrtpdf.m

function [p]=sqrtpdf(x, fun1, fun2, varargin)
%
%[p]=sqrtpdf(x, 'fun1', 'fun2', varargin)
%[p]=sqrtpdf(x, '', '', mu1, mu2, sigma1, sigma2)
%
% Generuje funkcje gestosci prawdopodobienstwa zmiennej losowej
% Z = sqrt( X.^2 + Y.^2 ), gdzie:
% X: zmienna losowa o rozkladzie tworzonym na bazie parametru 'x'
% przez funkcje 'fun1'
% Y: jak wyzej, rozklad tworzony przez 'fun2' na bazie
% parametru 'x'
% Zastosowano metode przejscia przez rozklad dwuwymiarowy tzn.
% /
% f(z)=|f(x,y)dy, zas
% /
% f(x,y)=f(x1*cos y1)*f(x1*sin y1)*|x1| jesli zmienne sa niezalezne.

120

% Skorzystano przy tym z przeksztalcenia g:R^2-->R^2 jak ponizej:
% z = sqrt(x.^2+y.^2) ---\ x = z*cos fi
% fi = arctan(y/x) ---/ y = z*sin fi
%
% PARAMETRY WEJSCIOWE:
% x - wektor kolumnowy zmiennej losowej X i Y
% fun1 - nazwa funkcji generujacej rozklad zmiennej X
% fun2 - nazwa funkcji generujacej rozklad zmiennej Y
% varargin - dodatkowe parametry dla 'fun1' i 'fun2'
% (kolejno arg1_fun1, arg1_fun2, arg2_fun1, ...)
% pierwszym parametrem 'fun1' i 'fun2' jest zawsze 'x'
% Jesli skorzysta sie z drugiej metody wywolania funkcji tzn. nie
% poda sie 'fun1' i 'fun2', wowczas zmienne losowe X i Y beda mialy
% rozklad Gaussa o parametrach kolejno mu1, sigma1 oraz mu2, sigma2
%
% WYJSCIE:
% p - wektor kolumnowy opisujacy rozklad zmiennej losowej Z
%
% Sprawdzone w srodowisku MATLAB 5.x
%

%tic
FILEN=256;
fi=linspace(0,2*pi,FILEN)';
Z=zeros(length(x),1);
p=zeros(length(x),1);

if isempty(fun1) | isempty(fun2)

for i=1:length(x)

Z=(exp(-0.5.*(varargin{3}.^(-2)).*((x(i,1).*cos(fi))-

varargin{1}).^2)./(sqrt(2*pi).*varargin{3})).*(exp(-
0.5.*(varargin{4}.^(-2)).*((x(i,1).*sin(fi))-
varargin{2}).^2)./(sqrt(2*pi).*varargin{4})).*abs(x(i,1));

p(i,1)=trapz(fi,Z);
end

else

for i=1:length(x)

Z=feval(fun1, x(i,1).*cos(fi),

varargin{1:2:length(varargin)}).*feval(fun2, x(i,1).*sin(fi),
varargin{2:2:length(varargin)}).*abs(x(i,1));

p(i,1)=trapz(fi,Z);
end

end
%toc

sumpdf.m

function [py, y]=sumpdf(x, px, N, todo)
%
%[py y]=sumpdf(x, px, N, 'mean') lub [py y]=sumpdf(x, px, N)
%[py y]=sumpdf(x, px, N, 'sum')
%
% Generuje funkcje gestosci prawdopodobienstwa zmiennej losowej
% bedacej suma N tych samych zmiennych losowych X:
% N
% Y = sum(X)
% i=1
%
% lub zmienna bedaca estymatorem wartosci sredniej tzn. jesli
% X - zmienna losowa wejsciowa, to

121

% N
% Y = 1/N * sum(X)
% i=1
%
% Zeby przyspieszyc obliczenia dla duzego N zastosowano przejscie
% przez funkcje tworzaca rozkladu
%
% PARAMETRY WEJSCIOWE:
% x - wektor kolumnowy zmiennej losowej X
% px - wektor kolumnowy funkcji gestosci
% prawdopodobienstwa p(x)
% N - patrz wzor na Y
% 'mean'
% 'sum' - okresla, czy funkcja wygeneruje rozklad sredniej
% ('mean'), czy sumy ('sum') zmiennych losowych.
% Parametr ten jest opcjonalny. Niepodanie go jest
% rownowazne podaniu 'mean'
% WYJSCIE:
% py - wektor kolumnowy opisujacy rozklad zmiennej losowej Y
% y - wektor kolumnowy zmiennej losowej Y
%
% Sprawdzone w srodowisku MATLAB 5.x
%

%tic
if exist('todo')==0 | isempty(todo)
todo='mean';
end

step=x(2)-x(1);
PX=fft(px,length(px).*N);
PY=PX.*step;
PY=PY.^(N-1);
PY=PY.*PX;

if strcmp(todo, 'mean')
y=x;
py=N.*abs(ifft(PY));
else
y=x.*N;
py=abs(ifft(PY));
end

%decimate 'py' N times
py=py(1:N:length(py));
%toc

122

LITERATURA:

1. Brandt S.: "Analiza danych - metody statystyczne i obliczeniowe", PWN,

Warszawa 1998

2. Czekała Z.: "Parada radarów", Dom Wydawniczy BELLONA, Warszawa 1999

3. Hellwig Z.: "Elementy rachunku prawdopodobieństwa i statystyki

matematycznej", PWN, Warszawa 1980, wydanie dziewiąte

4. Helstrom C. W.: "Statystyczna teoria detekcji", WNT, Warszawa 1964

5. Kącki T., Szyszkiewicz J., Majewski T., Słomiński L.: "Statystyczna teoria

radiolokacji", Wydział Wydawniczy WAT, Warszawa 1966

6. Knoch L., Ekiert T.: "Modulacja i detekcja", WKŁ, Warszawa 1979

7. Mrozek B., Mrozek Z.: "MATLAB, uniwersalne środowisko do obliczeń

naukowo-technicznych", Wydawnictwo PLJ, Warszawa 1996

8. Pacut A.: "Prawdopodobieństwo - teoria, modelowanie probabilistyczne w

technice", WNT, Warszawa 1985

9. Papoulis A.: "Prawdopodobieństwo, zmienne losowe i procesy stochastyczne",

WNT, Warszawa 1972

10. Seidler J.: "Statystyczna teoria odbioru sygnałów", PWN, Warszawa - Wrocław

1963

11. Skolnik M. I.: "Radar Handbook", Mc Graw-Hill, New York 1970

12. Smirnow N. W., Dunin-Barkowski I. W.: "Kurs rachunku prawdopodobieństwa i

statystyki matematycznej dla zastosowań technicznych", PWN, Warszawa

1969, wydanie drugie

13. Szabatin J.: "Podstawy teorii sygnałów", WKŁ, Warszawa 1982

14. MATLAB 5.2 User’s Guide

Wyszukiwarka

Podobne podstrony:
wykł-śc, POGODA- to chwilowy stan atmosfery na pewnym obszarze, okre˙lony przez uk˙ad powi˙zanych ze
Analiza podstawowych uk+éad+-w dyskretnych
BUD PRZEGRODA, Sprawdzi˙ pod wzgl˙dem cieplno-wilgotno˙ciowym przegrod˙ budowlan˙ pionow˙ o nast˙puj
uk ad pokarmowy
ANALIZA PRZYCZYN WYBUCHU WYBRANEJ WOJNY NA 3 POZIOMACH
uk-ad krwionoÂny. aq, Biomedyczne podstawy rozwoju i wychowania
UK AD LIMFATYCZNY, rodzaje i zasady masażu
uk+éad kr¦ů+ enia
Generatory drgan sinusoidalnych1, Celem ˙wiczenia jest zapoznanie si˙ z wybranymi podstawowymi uk˙ad
PRZEGR 1, Sprawdzi˙ pod wzgl˙dem cieplno-wilgotno˙ciowym przegrod˙ budowlan˙ pionow˙ o nast˙puj˙cym
TEATR OPRACOWANIA I sem, Raszewski- Uk+éad S, Raszewski - układ S
LABORATORIUM-NAPĘDÓW ELEK, Naped, UK˙AD DO REGULACJI PR˙DKO˙CI OBROTOWEJ
Parkinsonizm - drżączka poraźna, Parkinsonizm - dr˙˙czka pora˙na - jest to zesp˙˙ chorobowy charakte
Anatomia Uk%c5%82ad krwiotw%c3%b3rczy 06 notatki
Patologia Uk%c5%82ad kr%c4%85%c5%bcenia 03 notatki
Uk éad krwiono Ťny i ch éonny
UK AD KRWIONO NY I CH ONNY, rozwiązania z dopytek
8 6, 1. POMIARY W UK˙ADZIE SZEREGOWYM RLC.

więcej podobnych podstron