Nieparametryczne widmo mocy
Z Brain-wiki
Wersja z dnia 16:09, 14 lis 2016 autorstwa Jarekz (dyskusja | edycje) (→Zadanie 6: Porównanie rozdzielczości i wariancji w periodogramie i w estymatorze Welcha)
Analiza_sygnałów_-_ćwiczenia/Fourier_4
Spis treści
- 1 Widmo mocy
- 2 Sygnały stochastyczne
- 2.1 Stacjonarność i ergodyczność
- 2.2 Zadanie 4: Estymacja widma sygnału stochastycznego
- 2.3 Oszacowanie błędu transformaty Fouriera dla białego szumu
- 2.4 Zadanie 5: Metoda Welcha
- 2.5 Zadanie 6: Porównanie rozdzielczości i wariancji w periodogramie i w estymatorze Welcha
- 2.6 Zadanie 7: Estymacja widma mocy metodą multitaper
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 F,P
# część testująca
fs =100
T = 3
f = 10.3
# sygnał testowy
dt = 1.0/fs
t = np.arange(0,T,dt)
s = np.sin(2*np.pi*f*t )
# okno prostokątne
okno = np.ones(len(s))
moc_w_czasie = ...
(F, moc_w_czestosci) = 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()
Periodogram: widmo mocy okienkowanego 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, czyli periodogram:
[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 któ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ć przez 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ści zaobserwowanych w chwili we wszystkich realizacjach [math]t_1[/math] ważona prawdopodobieństwem wystąpienia tej realizacji.
Poniżej mamy przykład wytwarzania procesu złożonego z dwóch realizacji po 50 próbek oraz estymowania jego wartości średniej. Każda próbka jest niezależną zmienną losową z rozkładu normalnego o średniej 0 i wariancji 1:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import pylab as py
t = np.arange(0,50,1)
# realizacja 1
x1 = np.random.randn(t.size)
# realizacja 2
x2 = np.random.randn(t.size)
# średnia procesu
xm = 0.5*(x1+x2)
# ilustracja
py.subplot(3,1,1)
py.stem(t,x1)
py.title('realizacja 1')
py.subplot(3,1,2)
py.stem(t,x2)
py.title('realizacja 2')
py.subplot(3,1,3)
py.stem(t,xm,'r')
py.title('średnia procesu')
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 estymatorach momentów statystycznych.
Zadanie 4: Estymacja widma 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 (stochastycznego).
- 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.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 8 11:45:15 2016
@author: admin
"""
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
import scipy.stats as st
from numpy.fft import rfft,rfftfreq
def widmo_mocy(s,Fs):
'''funkcja oblicza widmo mocy sygnału rzeczywistego i oś częstości
s - sygnał
Fs - częstość próbkowania
zwraca dodatnią część widma
'''
...
return (F,P)
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania
Domyślnie wytwarzany jest sygnał reprezentujący
1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
'''
dt = 1.0/Fs
t = np.arange(0,T,dt)
s = np.sin(2*np.pi*f*t + phi)
return (s,t)
def szum(mu =0 , sigma = 1, T = 1, Fs = 128):
'''szum gaussowski o zadanej:
średniej mu
wariancji sigma**2
długości T,
częstości próbkowania Fs
'''
dt = 1.0/Fs
t = np.arange(0,T,dt)
s = ...
return (s,t)
def dwadziescia_realizacji(FS):
'''
* wygeneruj 20 realizacji sygnału będącego sumą sinusoidy (f=20Hz, T=1s, Fs =100Hz) i szumu gassowskiego
* dla każdej realizacji oblicz widmo mocy
* wykreśl wszystkie otrzymane widma na wspólnym wykresie
'''
for i in range(20):
(s,t) = ... #realizacja sinusa
(sz,t) = ...#realizacja szumu
syg = ...# sygnał będący sumą powyższych
(F, moc_w_czestosci) = widmo_mocy(syg, Fs=FS)
py.plot(F,moc_w_czestosci)
py.show()
def srednie_widmo(FS):
'''
# Skonstruuj funkcję rysującą średnie widmo wraz z 95% przedziałem ufności.
'''
zbior_widm = np.zeros((20,FS/2+1))
for i in range(20):
(s,t) = ... #realizacja sinusa
(sz,t) = ...#realizacja szumu
syg = ...# sygnał będący sumą powyższych
(moc_w_czestosci, F) = widmo_mocy(syg, Fs=FS)
zbior_widm[i][:] = ...#zapamiętaj widmo i-tej realizacji
srednie_w = ...# usrednij widma po realizacjach
przedzial_d = np.zeros(len(F)) # tablice na dolną i górną granicę przedziału ufności
przedzial_g = np.zeros(len(F))
for f in F: # dla każdej częstości znajdujemy granice przedziałów ufności
przedzial_d[f] = st.scoreatpercentile(..., 2.5)
przedzial_g[f] = st.scoreatpercentile(..., 97.5)
py.plot(F,srednie_w,'r') # rysujemy średnią
py.plot(F,przedzial_d,'b')# rysujemy granicę dolną
py.plot(F,przedzial_g,'b')# rysujemy granicę górną
py.show()
FS =100.0
dwadziescia_realizacji(FS)
srednie_widmo(FS)
Oszacowanie błędu transformaty Fouriera dla białego szumu
- Niech [math]x(t)[/math] - sygnał stochastyczny, 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ą
- Wówczas, część rzeczywista [math]X_R(f)[/math] i urojona [math]X_I(f)[/math] są nieskorelowanymi zmiennymi losowymi 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.
- 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 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],
- Oznacza to, ż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].
- Druga
- to 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.
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
from numpy.fft import rfft, rfftfreq
from scipy.signal import welch
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania
Domyślnie wytwarzany jest sygnał reprezentujący
1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
'''
dt = 1.0/Fs
t = np.arange(0,T,dt)
s = np.sin(2*np.pi*f*t + phi)
return (s,t)
def periodogram(s, okno , F_samp):
'''peiodogram sygnału s
okno - synał będzie przez nie przemnożony w czasie
F_samp- częstość próbkowania'''
s = s*okno
N_fft = len(s)
S = rfft(s,N_fft)
P = S*S.conj()/np.sum(okno**2)
P = P.real/Fs # P i tak ma zerowe wartości urojone, ale trzeba ykonać konwersję typów
F = rfftfreq(N_fft, 1/F_samp)
if len(s)%2 ==0: # dokładamy moc z ujemnej części widma
P[1:-1] *=2
else:
P[1:] *=2
return (F,P)
Fs = 100.0
(x,t) = sin(f = 3.1, T =20, Fs = Fs, phi = 0)
N = len(x) # długość sygnału
okno = np.hamming(N)
okno/=np.linalg.norm(okno)
py.subplot(2,1,1)
(F,P) = periodogram(x,okno,Fs)
py.plot(F,P)
py.title('periodogram'+' energia: '+ str(np.sum(P)))
#periodogram
py.subplot(2,1,2)
Nseg =20
N_s = N/Nseg
okno = np.hamming(N_s)
okno/=np.linalg.norm(okno)
(F, P) = welch(...)
py.plot(F,P)
py.title('periodogram Welcha'+' energia: '+ str(Nseg*np.sum(P)))
py.show()
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.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
from numpy.fft import rfft, rfftfreq
from scipy.signal import welch
import scipy.stats as st
from matplotlib.font_manager import FontProperties
font = FontProperties('Arial')
def periodogram(s, okno , F_samp):
'''peiodogram sygnału s
okno - synał będzie przez nie przemnożony w czasie
F_samp- częstość próbkowania'''
s = s*okno
N_fft = len(s)
S = rfft(s,N_fft)
P = S*S.conj()/np.sum(okno**2)
P = P.real/Fs # P i tak ma zerowe wartości urojone, ale trzeba ykonać konwersję typów
F = rfftfreq(N_fft, 1/F_samp)
if len(s)%2 ==0: # dokładamy moc z ujemnej części widma
P[1:-1] *=2
else:
P[1:] *=2
return (F,P)
def realizacja(T,Fs):
dt = 1.0/Fs
t = np.arange(0,T,dt)
# jedna realizacja sygnału będącego sumą sinusoidy (f = 20 Hz, T = 10 s, Fs = 100 Hz) i szumu gaussowskiego o std 5-krotnie większym niż amplituda sinusoidy
s =...
return s
T=10.0
Fs = 100.0
N = T*Fs
okno = np.blackman(N) # okno blakmana dla periodogramu
ile_okien = 10
Nw = N/ile_okien
okno_w = ...#okno blackmana dla welcha
N_rep = 100
S_perio = np.zeros((N_rep,N/2+1)) # uwaga, to jest dobrze tylko dla Fs parzystych
S_welch = np.zeros((N_rep,Nw/2+1)) # uwaga, to jest dobrze tylko dla Fs parzystych
for i in range(N_rep):#wygeneruj 100 realizacji sygnału
s = ...#sygnał
(F_p, P) = ...# jego periodogram
S_perio[i,:] = ...# zachowaj periodogram tej realizacji na później
(F_w, P_w) = ... # widmo Welcha
S_welch[i,:] = ...# zachowaj widmo Welcha tej realizacji na później
# poniżej robimy wykresy
py.figure(1)
py.subplot(3,1,1)
py.plot(F_p,np.mean(S_perio,axis = 0),'r',
F_p, st.scoreatpercentile(S_perio, 2.5,axis =0),'b',
F_p, st.scoreatpercentile(S_perio,97.5,axis =0),'b' )
py.title(u'Periodogram: średnie widmo realizacji wraz z 95% CI', fontproperties = font)
py.ylabel('V^2/Hz')
py.subplot(3,1,2)
py.plot(F_w,np.mean(S_welch,axis = 0),'r',
F_w, st.scoreatpercentile(S_welch, 2.5,axis =0),'b',
F_w, st.scoreatpercentile(S_welch,97.5,axis =0),'b' )
py.title('Welch: średnie widmo realizacji wraz z 95% CI', fontproperties = font)
py.ylabel('V^2/Hz')
py.subplot(3,1,3)
py.plot(F_p, ...) # błąd względny ( std/mean *100)w % dla periodogramu
py.plot(F_w, ...) # błąd względny w % dla Welcha
py.ylim([0,250])
py.xlabel('Częstość Hz', fontproperties = font)
py.ylabel('%')
py.legend(('periodogram','Welch'))
py.title('Błąd względny', fontproperties = font)
py.show()
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
