Nieparametryczne widmo mocy
Analiza_sygnałów_-_ćwiczenia/Fourier_4
Spis treści
- 1 Widmo mocy
- 2 Sygnały stochastyczne
Widmo mocy
Obliczanie mocy sygnału
Zadanie 1: Moc i energia sygnału w dziedzinie czasu
Proszę:
- wygenerować sygnał sinusoidalny [math]s[/math] o amplitudzie 1, częstości 10 Hz, trwający 0.3 sekundy i próbkowany z częstością 1000 Hz.
- narysować ten sygnał przy pomocy funkcji pylab.stem,
- obliczyć i narysować przebieg mocy w czasie [math]P_t = s_t^2[/math]: moc w danej chwili to kwadrat wartości próbki sygnału
- obliczyć energię tego sygnału [math]E = \sum_t P_t \Delta t [/math]: energia to suma mocy mnożonej przez przyrosty czasu między próbkami
Zadanie 2: Moc i energia sygnału w dziedzinie czasu i częstości
- Proszę uzupełnić i przetestować funkcję realizującą poniższy algorytm estymacji widma mocy.
- Następnie proszę obliczyć energię oraz wyświetlić przebieg widma mocy dla sygnału z Zadania 1.
- Sprawdzić czy energia zależy od częstości próbkowania i od długości sygnału
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
from numpy.fft import rfft,rfftfreq
def widmo_mocy(s,Fs,okno):
okno= # znormalizuj okno
s = # zokienkuj sygnał
S = # Oblicz transformatę Fouriera sygnału przy pomocy funkcji <tt>rfft</tt>
P = # Oblicz moc jako iloczyn unormowanej transformaty i jej sprzężenia zespolonego.
P = # Unormuj widmo dzieląc przez częstość próbkowania
P = # Do dalszych operacji wybierz tylko część rzeczywistą mocy.
if len(s)%2 ==0: # dokładamy moc z ujemnej części widma
P[1:-1] *=2
else:
P[1:] *=2
F = # Korzystając z funkcji <tt>rfftfreq</tt> obliczamy częstości, dla których policzone są współczynniki Fouriera.
return P,F
# część testująca
fs =100
T = 3
f = 10.3
# okno prostokątne
okno = np.ones(len(s))
# sygnał testowy
dt = 1.0/Fs
t = np.arange(0,T,dt)
s = np.sin(2*np.pi*f*t )
moc_w_czasie = ...
(moc_w_czestosci, F) = widmo_mocy(s, Fs=fs, okno = okno)
dt = 1/fs
energia_w_czasie = ...
energia_w_czestosci = ...
py.subplot(3,1,1)
py.plot(t,s)
py.title('Sygnal')
py.subplot(3,1,2)
py.plot(t,moc_w_czasie)
py.title('moc w czasie, energia w czasie: ' +str(energia_w_czasie))
py.subplot(3,1,3)
py.plot(F,moc_w_czestosci)
py.title('moc w czestosci, energia w czestosci: ' +str(energia_w_czestosci))
py.show()
Okienkowanie a widmo mocy: periodogram
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 3: Obliczanie periodogramu
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:
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:
ś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.
Zadanie 4: 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 = 20 Hz, T = 1 s, Fs = 100 Hz) 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.
Zadanie 5: Metoda Welcha
Proszę zapoznać się zaimplementowaną w bibliotece scipy.signal funkcją welch. 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.
Zadanie 6: Porównanie rozdzielczości i wariancji w periodogramie i w estymatorze Welcha
- wygeneruj 100 realizacji sygnału będącego sumą sinusoidy (f = 20 Hz, T = 10 s, Fs = 100 Hz) i szumu gaussowskiego
- dla każdej realizacji oblicz widmo mocy za pomocą periodogramu okienkowanego oknem Blackmana
- wykreśl wszystkie otrzymane widma na wspólnym wykresie (subplot(2,1,1))
- 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.
Zadanie 7: 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:
- Oblicz maksymalną liczbę okienek K = 2*NW-1
- Oblicz długość sygnału
- wygeneruj serię okienek dpss
- 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
- uśrednij widma otrzymane dla wszystkich okienek
- 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.
*
Analiza_sygnałów_-_ćwiczenia/Fourier_4