Nieparametryczne widmo mocy

Widmo mocy

Moc

Moc chwilowa sygnału przez analogię do układów elektrycznych o jednostkowym oporze jest w analizie sygnałów przyjęta jako kwadraty próbek ([math]P = I^2 R = \frac{U^2}{R}[/math]). Oznaczmy sygnał [math]x[n][/math], wówczas jego moc wyraża się wzorem:

[math]P[n]=x[n]^2[/math],

a energia wzorem:

[math]E = \sum _n{x[n]^2}[/math]

Zadanie

Przy pomocy funkcji napisanych na poprzednich zajęciach proszę wygenerować sygnał sinusoidalny o amplitudzie 1, częstości 10 Hz, trwający 0.3 sekundy i próbkowany z częstością 1000 Hz. Proszę narysować ten sygnał przy pomocy funkcji pylab.stem, obliczyć i narysować przebieg mocy w czasie, obliczyć energię tego sygnału.

Widmo mocy: tw. Plancherela i tw. Parsevala

Twierdzenia te omawiane i dowodzone były na wykładzie. Tutaj, tylko krótko przypomnijmy sobie:

Twierdzenie Plancherela

Jeśli [math]X[k][/math] i [math]Y[k][/math] są transformatami [math]x[n][/math] i [math]y[n][/math] odpowiednio to:

[math] \sum _{n=0}^{N-1} x[n]y^*[n] = \frac{1}{N} \sum _{k=0}^{N-1} X[k] Y^*[k] [/math]

gwiazdka oznacza sorzężenie zespolone.

Twierdzenie Parsevala

jest specjalnym przypadkiem twierdzenia Plancherela:

[math] \sum _{n=0}^{N-1} \left|x[n]\right|^2 = \frac{1}{N} \sum _{k=0}^{N-1} \left|X[k]\right|^2. [/math]

Twierdzenie to upoważnia nas do utożsamiania kwadratów wartości bezwzględnej składowych transformaty Fouriera z mocą niesioną przez odpowiadające im składowe.

Zadanie: funkcja do estymacji widma mocy

Proszę napisać i przetestować funkcję realizującą następujący algorytm estymacji widma mocy:

1. Jako parametry funkcja przyjmuje sygnał s i częstość próbkowania tego sygnału
2. Oblicz transformatę Fouriera sygnału przy pomocy fft
3. Unormuj otrzymaną transformatę przez pierwiastek z ilości próbek sygnału
4. Oblicz moc jako iloczyn unormowanej transformaty i jej sprzężenia zespolonego. W pythonie sprzężenie zespolone liczby S otrzymujemy korzystając z metody S.conj(). W wyniku powyższego iloczynu dostaniemy liczby zespolone o zerowej części urojonej (proszę to sprawdzić).
5. Do dalszych operacji wybierz tylko część rzeczywistą mocy. Część rzeczywistą liczby zespolonej M pobieramy w następujący sposób: czesc_rzeczywista_M = M.real
6. Korzystając z funkcji fftfreq obliczamy częstości, dla których policzone są współczynniki Fouriera.
7. Przy pomocy funkcji fftshift porządkujemy kolejność wektorów mocy i częstości, tak aby częstości były reprezentowane od -częstości Nyquista, przez 0, do +częstości Nyquista
8. zwracamy uporządkowane wektory mocy i częstości

Testowanie:

1. przy pomocy funkcji sin napisanej na drugich ćwiczeniach wygeneruj sinusoidę o długości 1 s, o częstości 10Hz, próbkowaną 100 Hz.
2. oblicz moc w czasie (trzeba podnieść wartość każdej próbki do kwadratu)
3. oblicz moc w częstości przy pomocy funkcji estymacji widma mocy
4. oblicz energię sygnału w czasie
5. oblicz energię sygnału w częstości
6. wykreśl:
   • sygnał
   • przebieg jego mocy w czasie
   • przebieg mocy w częstości

*

Okienkowanie a widmo mocy: periodogram

Przypomnijmy wzór na dyskretną transformatę Fouriera DFT zaimplementowaną w FFT:

[math]S[k] = \sum_{n=0}^{n-1} s[n] \exp\left\{-2\pi i{nk \over N}\right\} \qquad k = 0,\ldots,N-1. [/math]

Na podstawie twierdzenia Parsevala możemy policzyć widmo mocy jako: [math] P[k] = \frac{1}{N} \left|S[k]\right|^2 [/math]

Jeśli do liczenia mocy chcielibyśmy posłużyć się techniką okiennkowania sygnału, to powinniśmy używać okienek znormalizowanych, czyli takich których energia jest równa 1, wtedy mnożenie przez okienko nie zaburzy estymaty energii sygnału.

Aby policzyć widmo mocy sygnału z zastosowaniem okienek wprowadzimy następujące symbole:

• sygnał: [math]s[n][/math]
• okienko: [math] w[n][/math]
• okienko znormalizowane: [math] \hat w[n] = \frac{1}{\sqrt{\sum_{n=0}^{N-1} (w[n])^2}}w[n][/math]
• widmo mocy sygnału okienkowanego:

[math] P[k] = \frac{1}{\sum_{n=0}^{N-1} (w[n])^2} \left|\sum_{n=0}^{N-1} s[n]w[n] e^{i\frac{2 \pi }{N} k n}\right|^2 [/math]

Zadanie

Proszę napisać funkcję obliczającą periodogram. Funkcja jako argumenty powinna przyjmować sygnał, okno (podane jako sekwencja próbek), i częstość próbkowania. Zwracać powinna widmo mocy i skalę osi częstości. Wewnątrz funkcja powinna implementować liczenie widma z sygnału okienkowanego znormalizowanym oknem.

Funkcję proszę przetestować obliczając dla funkcji sinus energię sygnału w dziedzinie czasu i w dziedzinie częstości. Testy proszę wykonać dla okna prostokątnego, Blackmana i Haminga.

*

Sygnały stochastyczne

Sygnał stochastyczny to taki sygnał, dla którego ciągu próbek nie da się opisać funkcją czasu. Kolejne próbki w takim sygnale to zmienne losowe. Można je opisać podając własności rozkładu, z k†órego pochodzą. Często w opisie takich zmiennych posługujemy się momentami rozkładów. Jak można sobie wyobrazić rozkłady, z których pochodzą próbki? Można sobie wyobrazić,że obserwowany przez nas sygnał stochastyczny to jedna z możliwych realizacji procesu stochastycznego. Jeśli [math]K[/math] jest zbiorem [math]k[/math] zdarzeń ([math]k \in K[/math]) i każde z tych zdarzeń ma przypisaną funkcję [math]x_k(t)[/math] zwaną realizacją procesu [math]\xi (t)[/math], to proces stochastyczny może być zdefiniowany jako zbiór funkcji:

[math] \xi (t) = \left\lbrace x_1(t),x_2(t),\dots , x_N(t) \right\rbrace [/math]

gdzie [math]x_k(t)[/math] są losowymi funkcjami czasu [math]t[/math].

Procesy stochastyczne można opisywać prze wartości oczekiwane liczone po realizacjach.

Dla przypomnienia wartość oczekiwaną liczymy tak:

[math] {\mu _x(t_1) = E\left[\xi (t_1) \right]= \lim _{N \rightarrow \infty }\sum _{k=1}^{N}{x_k(t_1)} p(x_k,t_1)} [/math]

średnia [math]\mu _x(t_1)[/math] procesu [math]\xi (t)[/math] w chwili [math]t_1[/math] to suma wartośći zaobserwowanych w chwili we wszystkich realizacjach [math]t_1[/math] ważona prawdopodobieństwem wystąpienia tej realizacji:

Stacjonarność i ergodyczność

Stacjonarność:

Jeśli dla procesu stochastycznego [math]\xi (t)[/math] wszystkie momenty są niezależne od czasu to jest on stajonarny w ścisłym sensie. Jeśli tylko średnia [math]\mu _x[/math] i autokorelacja [math]R_x(\tau )[/math] nie zależą od czasu to proces jest stacjonarny w słabym sensie, co dla wielu zastosowań jest wystarczające.

Ergodyczność:

Proces jest ergodyczny jeśli jego średnie po czasie i po realizacjach są sobie równe. Oznacza to, że dla takiego procesu jedna realizacja jest reprezentatywna i zawiera całą informację o tym procesie.

Założenie o sygnale, że jest stacjonarny i ergodyczny pozwala zamienić sumowanie po realizacjach na sumowanie po czasie w estymatory momentów statystycznych.

Transformata Fouriera sygnału stochastycznego

Bardzo często musimy oszacować widmo mocy sygnału zawierającego znaczny udział szumu.

Poniższe ćwiczenie ilustruje niepewność szacowania pików w widmie otrzymanym z transformaty Fouriera dla sygnału zawierającego szum.

• wygeneruj 20 realizacji sygnału będącego sumą sinusoidy (f=20Hz, T=1s, Fs =100Hz) i szumu gaussowskiego
• dla każdej realizacji oblicz widmo mocy
• wykreśl wszystkie otrzymane widma na wspólnym wykresie

Proszę obejrzeć otrzymane widma.

• Zaobserwuj jakiego rzędu jest niepewność wyniku.
• Czy podobny problem występuje dla sygnału bez szumu?
• Skonstruuj funkcję rysującą średnie widmo wraz z przedziałem ufności.

*

Oszacowanie błędu transformaty Fouriera dla białego szumu

Dla sygnału stochastycznego [math]x(t)[/math], którego kolejne próbki pochodzą z niezależnych rozkładów normalnych (biały szum), jego transformata Fouriera [math]X(f)[/math] jest liczbą zespoloną, której część rzeczywista [math]X_R(f)[/math] i urojona [math]X_I(f)[/math] mogą być uznane za nieskorelowane zmienne losowe o średniej zero i równych wariancjach. Ponieważ transformata Fouriera jest operacją liniową więc składowe [math]X_R(f)[/math] i [math]X_I(f)[/math] mają rozkłady normalne. Zatem wielkość: [math] P(f) = |X(f)|^2 = X_R^2(f) + X_I^2(f) [/math] jest sumą kwadratów dwóch niezależnych zmiennych normalnych. Wielkość ta podlega zatem rozkładowi [math]\chi^2[/math] o dwóch stopniach swobody. Możemy oszacować względny błąd [math]P(f_1) [/math] dla danej częstości [math]f_1[/math]:

[math]\epsilon_r= \sigma_{P_{f_1}}/\mu_{P_{f_1}}[/math]

Dla rozkładu [math]\chi_2^2[/math]: [math]\sigma^2 = 2n[/math] i [math]\mu = n[/math], gdzie [math]n[/math] jest ilością stopni swobody. W naszym przypadku [math]n =2[/math] więc mamy [math]\epsilon_f = 1[/math], co oznacza, że dla pojedynczego binu częstości w widmie [math]P(f)[/math] względny błąd wynosi 100%.

Aby zmniejszyć ten błąd trzeba zwiększyć ilość stopni swobody. Są generalnie stosowane dwie techniki. Pierwsza to uśrednianie sąsiednich binów częstości. Otrzymujemy wówczas wygładzony estymator mocy [math]\hat{P}_k[/math]:

[math]\hat{P}_k = \frac{1}{l}[P_k + P_{k+1} + \dots + P_{k+l-1}][/math]

Zakładając, że biny częstości [math]P_i[/math] są niezależne estymator [math]P_k[/math] ma rozkład [math]\chi^2[/math] o ilości stopni swobody równej [math]n= 2l[/math]. Względny błąd takiego estymatora to: [math]\epsilon_r= \sqrt{\frac{1}{l}}[/math].

Innym sposobem poprawy estymatora mocy jest podzielenie sygnału na fragmenty, obliczenie periodogramu dla każdego fragmentu, a następnie zsumowanie otrzymanych wartości:

[math]\hat{P}_k=[P_{k,1}+P_{k,2}+\dots+P_{k,j}+\dots+P_{k,q}][/math]

gdzie [math]S_{k,j}[/math] jest estymatą składowej o częstości [math]k[/math] w oparciu o [math]j-ty[/math] fragment sygnału. Ilość stopni swobody wynosi w tym przypadku [math]q[/math] zatem względny błąd wynosi: [math]\epsilon_r = \sqrt{\frac{1}{q}}[/math].

Zauważmy, że w obu metodach zmniejszamy wariancję estymatora kosztem rozdzielczości w częstości.

Metoda Welcha

Proszę napisać i przetestować funkcję implementującą metodę Welcha estymacji widma mocy. Algorytm Welcha:

1. sygnał [math]x[/math] o długości N jest dzielony na segmenty, każdy o długości [math]N_s[/math]. Odcinki mogą na siebie zachodzić na [math]N_z[/math] punktów. Czyli są względem siebie przesunięte o [math]N_p = N_s-N_z[/math].
2. z każdego segmentu liczony jest okienkowany periodogram
3. periodogramy są sumowane
4. wynik dzielony jest przez efektywne wykorzystanie każdego kawałka sygnału w estymacie: K_eff = dlogosc_okna * ilosc_okien / dlugosc_sygnalu

Funkcję proszę przetestować obliczając dla funkcji sinus energię sygnału w dziedzinie czasu i w dziedzinie częstości. Testy proszę wykonać dla okna prostokątnego, Blackmana i Haminga.

* 

Porównanie rozdzielczości i wariancji w periodogramie i w estymatorze Welcha

1. wygeneruj 100 realizacji sygnału będącego sumą sinusoidy (f=20Hz, T=10s, Fs =100Hz) i szumu gaussowskiego
2. dla każdej realizacji oblicz widmo mocy za pomocą periodogramu okienkowanego oknem Blackmana
3. wykreśl wszystkie otrzymane widma na wspólnym wykresie (subplot(2,1,1))
4. Powtórz krok 2) dla estymatora Welcha z oknem Blackmana o długości 1/10 długości sygnału przesuwanym co 2 punkty, otrzymane widma wykreśl na wspólnym wykresie (subplot(2,1,2))

• Co można powiedzieć o rozdzielczości i względnym błędzie obu metod?

bl_wzg = np.std(S,axis = 0)/np.mean(S,axis = 0) gdzie S jest tablicą zawierającą widma dla każdej z realizacji.

*

Wielookienkowa metoda Thomsona

Metoda ta Spectrum estimation and harmonic analysis znana jest pod anglojęzyczną nazwą multitaper.

Można ją opisać poniższym algorytmem:

• wygeneruj sekwencję ortogonalnych okienek charakteryzujących się minimalnymi wyciekami widma (stosunek energii w piku centralnym do energii w listkach bocznych jest wysoki). Sekwencja takich okien nazywana jest discrete prolate spheroidal sequences (DPSS) lub sekwencją Slepiana.
• oblicz widmo sygnału okienkowanego każdym z okien w sekwencji
• uśrednij otrzymane widma

Kolejne dwa zadania służą zapoznaniu się z tą metodą.

Własności okienek DPSS

Do generacji sekwencji okienek DPSS wykorzystamy moduł Gendpss.py. Proszę go zapisać w swoim katalogu roboczym. Importujemy go do naszych programów tak jak każdy inny moduł np.:

import Gendpss as dpss

Funkcja potrzebna nam z tego modułu to gendpss(). Funkcja ta wytwarza obiekt reprezentujący konkretną sekwencję DPSS. Wywołujemy ją następująco:

w = dpss.gendpss(N,NW,K)

gdzie: N: długość okna,

NW: iloczyn czas-szerokość pasma

K: ile okien w sekwencji

Po powyższym wywołaniu obiekt w posiada dwie interesujące nas tablice:

• w.lambdas - to wartości własne okienek. Są one miarą koncentracji energii w piku głównym, jest to zatem miara jakości okienka (dobre okienka mają wartości własne bliskie 1). Zgodnie z teorią takich wartości powinno być 2*NW-1.
• w.dpssarray[i] - i-te okienko.

Polecenia:

Proszę:

• wygenerować okienka o długości 256, NW = 2. Ilość okien K raz ustalić na 3 (2*NW-1) a drugi raz na 5. Dla ilu okienek ich wartości własne są bliskie 1?
• narysować przebieg czasowy okienek
• sprawdzić czy energia okienek jest znormalizowana do 1.
• sprawdzić czy kolejne okienka są do siebie ortogonalne. W tym celu należy obliczyć iloczyn skalarny pomiędzy kolejnymi okienkami (np.sum(w.dpssarray[i]*w.dpssarray[j])).
• wywrysować widma okienek analogicznie jak w tym ćwiczeniu

*

Estymacja widma mocy metodą multitaper

Proszę napisać funkcję do estymacji mocy metodą multitaper. Funkcja powinna pobierać następujące argumenty: sygnał, iloczyn NW, częstość próbkowania sygnału. Funkcja powinna zwracać krotkę (S,F) gdzie S widmo mocy, F skala częstości. Przykładowe wywołanie takiej funkcji powinno wyglądać tak: (S,F) = mtm(s, NW = 3, Fs = 128) Algorytm do zastosowania wewnątrz funkcji:

1. Oblicz maksymalną liczbę okienek K = 2*NW-1
2. Oblicz długość sygnału
3. wygeneruj serię okienek dpss
4. dla każdego z otrzymanych okienek oblicz widmo mocy iloczynu tego okienka i sygnału. Dla i-tego okienka będzie to: Si = np.abs(fft(s*w.dpssarray[i]))**2
5. uśrednij widma otrzymane dla wszystkich okienek
6. wygeneruj oś częstości (fftfreq)

Działanie funkcji sprawdź estymując i wykreślając widmo sinusoidy np. o częstości 10 Hz, czasie trwania 1s, próbkowanej 100Hz z dodanym szumem gaussowskim o średniej 0 i wariancji 1. Sprawdź także zachowanie energii przez tą estymatę. Dla porównania na tym samym wykresie dorysuj widmo otrzymane przez periodogram z oknem prostokątnym.

*