Nieparametryczne widmo mocy: Różnice pomiędzy wersjami
(Nie pokazano 45 wersji utworzonych przez 2 użytkowników) | |||
Linia 64: | Linia 64: | ||
==== Zadanie 2: Moc i energia sygnału w dziedzinie czasu i częstości ==== | ==== Zadanie 2: Moc i energia sygnału w dziedzinie czasu i częstości ==== | ||
− | Proszę uzupełnić i przetestować funkcję realizującą | + | * 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 | ||
− | + | <source lang = python> | |
− | # | + | #!/usr/bin/env python3 |
− | # Oblicz transformatę Fouriera sygnału przy pomocy funkcji <tt>rfft</tt> | + | # -*- coding: utf-8 -*- |
− | + | ||
− | # Oblicz moc jako iloczyn unormowanej transformaty i jej sprzężenia zespolonego. | + | import pylab as py |
− | # Do dalszych operacji wybierz tylko część rzeczywistą mocy. | + | import numpy as np |
− | # Korzystając z funkcji <tt>rfftfreq</tt> obliczamy częstości, dla których policzone są współczynniki Fouriera. | + | 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 | ||
− | |||
− | <source | + | # 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() | ||
+ | </source> | ||
+ | |||
+ | <!-- | ||
#!/usr/bin/env python3 | #!/usr/bin/env python3 | ||
# -*- coding: utf-8 -*- | # -*- coding: utf-8 -*- | ||
+ | """ | ||
+ | Created on Mon Nov 7 12:21:33 2016 | ||
+ | |||
+ | @author: admin | ||
+ | """ | ||
import pylab as py | import pylab as py | ||
Linia 85: | Linia 137: | ||
from numpy.fft import rfft,rfftfreq | from numpy.fft import rfft,rfftfreq | ||
− | def widmo_mocy(s,Fs): | + | def widmo_mocy(s,Fs,okno): |
− | |||
okno= okno/np.linalg.norm(okno) | okno= okno/np.linalg.norm(okno) | ||
s = s*okno | s = s*okno | ||
S = rfft(s) | S = rfft(s) | ||
− | P = S*S.conj()/Fs | + | P = S*S.conj() |
+ | P = P/Fs | ||
P = P.real | P = P.real | ||
if len(s)%2 ==0: | if len(s)%2 ==0: | ||
Linia 97: | Linia 149: | ||
P[1:] *=2 | P[1:] *=2 | ||
F = rfftfreq(len(s), 1./Fs) | F = rfftfreq(len(s), 1./Fs) | ||
− | return | + | return F,P |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
+ | # część testująca | ||
fs =100 | fs =100 | ||
T = 3 | 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 = s**2 | moc_w_czasie = s**2 | ||
− | (moc_w_czestosci | + | (F, moc_w_czestosci) = widmo_mocy(s, Fs=fs,okno=okno) |
dt = 1/fs | dt = 1/fs | ||
Linia 129: | Linia 181: | ||
py.title('moc w czestosci, energia w czestosci: ' +str(energia_w_czestosci)) | py.title('moc w czestosci, energia w czestosci: ' +str(energia_w_czestosci)) | ||
py.show() | py.show() | ||
− | |||
− | |||
− | |||
− | |||
− | |||
--> | --> | ||
− | === | + | ===Periodogram: widmo mocy okienkowanego sygnału === |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
Aby policzyć widmo mocy sygnału z zastosowaniem okienek wprowadzimy następujące symbole: | Aby policzyć widmo mocy sygnału z zastosowaniem okienek wprowadzimy następujące symbole: | ||
* sygnał: <math>s[n]</math> | * sygnał: <math>s[n]</math> | ||
Linia 155: | Linia 189: | ||
* okienko znormalizowane: <math> \hat w[n] = \frac{1}{\sqrt{\sum_{n=0}^{N-1} (w[n])^2}}w[n]</math> | * okienko znormalizowane: <math> \hat w[n] = \frac{1}{\sqrt{\sum_{n=0}^{N-1} (w[n])^2}}w[n]</math> | ||
<!--(w szczególnym przypadku okienka prostokątnego normalizacja ta daje <math>1/N^2</math> występujące we wzorze na moc)--> | <!--(w szczególnym przypadku okienka prostokątnego normalizacja ta daje <math>1/N^2</math> występujące we wzorze na moc)--> | ||
− | * widmo mocy sygnału okienkowanego: | + | * widmo mocy sygnału okienkowanego, czyli periodogram: |
<math> | <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 | 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 | ||
Linia 161: | Linia 195: | ||
====Zadanie 3: Obliczanie periodogramu ==== | ====Zadanie 3: Obliczanie periodogramu ==== | ||
− | Proszę napisać funkcję obliczającą periodogram. | + | * 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. | + | ** 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. | + | * 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. |
<!-- | <!-- | ||
Linia 199: | Linia 233: | ||
s = s*okno | s = s*okno | ||
N_fft = len(s) | N_fft = len(s) | ||
− | S = | + | S = rfft(s,N_fft) |
P = S*S.conj()/np.sum(okno**2) | P = S*S.conj()/np.sum(okno**2) | ||
− | P = P.real # P i tak ma zerowe | + | P = P.real # P i tak ma zerowe wartości urojone, ale trzeba ykonać konwersję typów |
− | F = | + | F = rfftfreq(N_fft, 1/F_samp) |
− | return | + | if len(s)%2 ==0: # dokładamy moc z ujemnej części widma |
+ | P[1:-1] *=2 | ||
+ | else: | ||
+ | P[1:] *=2 | ||
+ | return P,F | ||
F_samp = 100.0 | F_samp = 100.0 | ||
Linia 245: | Linia 283: | ||
==Sygnały stochastyczne == | ==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 [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Zmienna_losowa|zmienne losowe]]. Można je opisać podając własności [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Rozk.C5.82ad_prawdopodobie.C5.84stwa|rozkładu]], z | + | 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 [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Zmienna_losowa|zmienne losowe]]. Można je opisać podając własności [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Rozk.C5.82ad_prawdopodobie.C5.84stwa|rozkładu]], z którego pochodzą. Często w opisie takich zmiennych posługujemy się [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Momenty|momentami rozkładów]]. |
Jak można sobie wyobrazić rozkłady, z których pochodzą próbki? | 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. | + | 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: | 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: | ||
<equation id="uid23"> | <equation id="uid23"> | ||
Linia 256: | Linia 294: | ||
gdzie <math>x_k(t)</math> są losowymi funkcjami czasu <math>t</math>. | gdzie <math>x_k(t)</math> są losowymi funkcjami czasu <math>t</math>. | ||
− | Procesy stochastyczne można opisywać | + | |
+ | Procesy stochastyczne można opisywać przez wartości oczekiwane liczone po realizacjach. | ||
Dla przypomnienia wartość oczekiwaną liczymy tak: | Dla przypomnienia wartość oczekiwaną liczymy tak: | ||
Linia 264: | Linia 303: | ||
</math> | </math> | ||
</equation> | </equation> | ||
− | średnia <math>\mu _x(t_1)</math> procesu <math>\xi (t)</math> w chwili <math>t_1</math> to suma | + | ś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: | ||
+ | <source lang = python> | ||
+ | #!/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') | ||
+ | </source> | ||
+ | |||
=== Stacjonarność i ergodyczność=== | === Stacjonarność i ergodyczność=== | ||
<dl> | <dl> | ||
Linia 277: | Linia 347: | ||
</dl> | </dl> | ||
− | Założenie o sygnale, że jest stacjonarny i ergodyczny pozwala zamienić sumowanie po realizacjach na sumowanie po czasie w | + | 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: | + | ===Zadanie 4: Estymacja widma sygnału stochastycznego=== |
Bardzo często musimy oszacować widmo mocy sygnału zawierającego znaczny udział szumu. | 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. | + | 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 | * wygeneruj 20 realizacji sygnału będącego sumą sinusoidy (f = 20 Hz, T = 1 s, Fs = 100 Hz) i szumu gaussowskiego | ||
Linia 292: | Linia 362: | ||
* Czy podobny problem występuje dla sygnału bez szumu? | * Czy podobny problem występuje dla sygnału bez szumu? | ||
* Skonstruuj funkcję rysującą średnie widmo wraz z [[WnioskowanieStatystyczne/_Przedzia%C5%82y_ufno%C5%9Bci|przedziałem ufności]]. | * Skonstruuj funkcję rysującą średnie widmo wraz z [[WnioskowanieStatystyczne/_Przedzia%C5%82y_ufno%C5%9Bci|przedziałem ufności]]. | ||
+ | <source lang = python> | ||
+ | #!/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 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 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 | ||
+ | okno = ... | ||
+ | (moc_w_czestosci, F) = periodogram(syg, Fs=FS,okno = okno) | ||
+ | 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) | ||
+ | </source> | ||
<!-- | <!-- | ||
<source lang = python> | <source lang = python> | ||
Linia 327: | Linia 490: | ||
t = np.arange(0,T,dt) | t = np.arange(0,T,dt) | ||
s = np.sin(2*np.pi*f*t + phi) | s = np.sin(2*np.pi*f*t + phi) | ||
− | return (s | + | return (t, s) |
def szum(mu =0 , sigma = 1, T = 1, Fs = 128): | def szum(mu =0 , sigma = 1, T = 1, Fs = 128): | ||
dt = 1.0/Fs | dt = 1.0/Fs | ||
t = np.arange(0,T,dt) | t = np.arange(0,T,dt) | ||
s = np.random.randn(len(t) )*sigma + mu | s = np.random.randn(len(t) )*sigma + mu | ||
− | return (s | + | return (t, s) |
def dwadziescia_realizacji(FS): | def dwadziescia_realizacji(FS): | ||
Linia 341: | Linia 504: | ||
''' | ''' | ||
for i in range(20): | for i in range(20): | ||
− | (s | + | (t, s) = sin(f=10,Fs=FS) |
− | (sz | + | (t, sz) = szum(Fs =FS) |
syg = s + sz | syg = s + sz | ||
(moc_w_czestosci, F) = widmo_mocy(syg, Fs=FS) | (moc_w_czestosci, F) = widmo_mocy(syg, Fs=FS) | ||
Linia 354: | Linia 517: | ||
zbior_widm = np.zeros((20,FS/2+1)) | zbior_widm = np.zeros((20,FS/2+1)) | ||
for i in range(20): | for i in range(20): | ||
− | (s | + | (t, s) = sin(f=10,Fs=FS) |
− | (sz | + | (t, sz) = szum(Fs =FS) |
syg = s + sz | syg = s + sz | ||
(moc_w_czestosci, F) = widmo_mocy(syg, Fs=FS) | (moc_w_czestosci, F) = widmo_mocy(syg, Fs=FS) | ||
Linia 371: | Linia 534: | ||
FS =100.0 | FS =100.0 | ||
− | + | dwadziescia_realizacji(FS) | |
srednie_widmo(FS) | srednie_widmo(FS) | ||
− | |||
− | |||
</source> | </source> | ||
Linia 380: | Linia 541: | ||
=== Oszacowanie błędu transformaty Fouriera dla białego szumu === | === 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), | |
− | <math> | + | * Jego transformata Fouriera <math>X(f)</math> jest liczbą zespoloną |
− | P(f) = |X(f)|^2 = X_R^2(f) + X_I^2(f) | + | * 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. |
− | </math> | + | * 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. |
− | jest sumą kwadratów dwóch niezależnych zmiennych normalnych. Wielkość ta podlega | + | * 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%. | ||
− | :<math>\hat{P}_k = \frac{1}{l}[P_k + P_{k+1} + \dots + P_{k+l-1}]</math> | + | * 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ą [https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.signal.welch.html 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. | Proszę zapoznać się zaimplementowaną w bibliotece scipy.signal funkcją [https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.signal.welch.html 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. | ||
− | + | <source lang = python> | |
− | |||
# -*- coding: utf-8 -*- | # -*- coding: utf-8 -*- | ||
− | |||
import pylab as py | import pylab as py | ||
import numpy as np | import numpy as np | ||
from numpy.fft import rfft, rfftfreq | from numpy.fft import rfft, rfftfreq | ||
− | + | from scipy.signal import welch | |
def sin(f = 1, T = 1, Fs = 128, phi =0 ): | def sin(f = 1, T = 1, Fs = 128, phi =0 ): | ||
Linia 438: | Linia 592: | ||
okno - synał będzie przez nie przemnożony w czasie | okno - synał będzie przez nie przemnożony w czasie | ||
F_samp- częstość próbkowania''' | F_samp- częstość próbkowania''' | ||
− | |||
s = s*okno | s = s*okno | ||
N_fft = len(s) | N_fft = len(s) | ||
− | S = rfft(s, | + | S = rfft(s,N_fft) |
− | P = S*S.conj()/np.sum(okno**2) | + | 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 | |
− | P = P.real # P i tak ma zerowe | ||
F = rfftfreq(N_fft, 1/F_samp) | 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) | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | return ( | ||
Fs = 100.0 | Fs = 100.0 | ||
(x,t) = sin(f = 3.1, T =20, Fs = Fs, phi = 0) | (x,t) = sin(f = 3.1, T =20, Fs = Fs, phi = 0) | ||
− | N = len(x) # | + | N = len(x) # długość sygnału |
− | |||
− | |||
− | |||
− | + | okno = np.hamming(N) | |
+ | okno/=np.linalg.norm(okno) | ||
py.subplot(2,1,1) | py.subplot(2,1,1) | ||
− | (P | + | (F,P) = periodogram(x,okno,Fs) |
py.plot(F,P) | py.plot(F,P) | ||
py.title('periodogram'+' energia: '+ str(np.sum(P))) | py.title('periodogram'+' energia: '+ str(np.sum(P))) | ||
Linia 484: | Linia 620: | ||
#periodogram | #periodogram | ||
py.subplot(2,1,2) | py.subplot(2,1,2) | ||
− | N_s = N/ | + | Nseg =20 |
− | + | N_s = N/Nseg | |
− | + | ||
okno = np.hamming(N_s) | okno = np.hamming(N_s) | ||
− | + | okno/=np.linalg.norm(okno) | |
− | + | (F, P) = welch(...) | |
− | (F, P) = welch( | ||
py.plot(F,P) | py.plot(F,P) | ||
− | py.title('periodogram Welcha'+' energia: '+ str(np.sum(P))) | + | py.title('periodogram Welcha'+' energia: '+ str(Nseg*np.sum(P))) |
py.show() | py.show() | ||
− | |||
</source> | </source> | ||
− | |||
− | |||
===Zadanie 6: Porównanie rozdzielczości i wariancji w periodogramie i w estymatorze Welcha=== | ===Zadanie 6: Porównanie rozdzielczości i wariancji w periodogramie i w estymatorze Welcha=== | ||
Linia 507: | Linia 639: | ||
* Co można powiedzieć o rozdzielczości i względnym błędzie obu metod? | * Co można powiedzieć o rozdzielczości i względnym błędzie obu metod? | ||
<tt>bl_wzg = np.std(S,axis = 0)/np.mean(S,axis = 0)</tt> gdzie S jest tablicą zawierającą widma dla każdej z realizacji. | <tt>bl_wzg = np.std(S,axis = 0)/np.mean(S,axis = 0)</tt> gdzie S jest tablicą zawierającą widma dla każdej z realizacji. | ||
− | |||
<source lang = python> | <source lang = python> | ||
+ | #!/usr/bin/env python3 | ||
# -*- coding: utf-8 -*- | # -*- coding: utf-8 -*- | ||
import pylab as py | import pylab as py | ||
import numpy as np | import numpy as np | ||
− | from numpy.fft import | + | from numpy.fft import rfft, rfftfreq |
− | + | from scipy.signal import welch | |
− | + | import scipy.stats as st | |
− | def | + | from matplotlib.font_manager import FontProperties |
− | ''' | + | |
− | + | font = FontProperties('Arial') | |
− | 1 | + | 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 | dt = 1.0/Fs | ||
t = np.arange(0,T,dt) | t = np.arange(0,T,dt) | ||
− | s = np. | + | # 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() | ||
+ | |||
+ | </source> | ||
+ | <!-- | ||
+ | <source lang = python> | ||
+ | #!/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 | ||
− | def periodogram(s, okno , | + | font = FontProperties('Arial') |
+ | def periodogram(s, okno , F_samp): | ||
'''peiodogram sygnału s | '''peiodogram sygnału s | ||
okno - synał będzie przez nie przemnożony w czasie | okno - synał będzie przez nie przemnożony w czasie | ||
− | + | F_samp- częstość próbkowania''' | |
− | |||
s = s*okno | s = s*okno | ||
N_fft = len(s) | N_fft = len(s) | ||
− | S = | + | S = rfft(s,N_fft) |
− | P = S*S.conj()/np.sum(okno**2) | + | 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 | |
− | P = P.real # P i tak ma zerowe | + | F = rfftfreq(N_fft, 1/F_samp) |
− | F = | + | if len(s)%2 ==0: # dokładamy moc z ujemnej części widma |
− | + | P[1:-1] *=2 | |
− | + | else: | |
− | + | P[1:] *=2 | |
− | + | return (F,P) | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | return ( | ||
def realizacja(T,Fs): | def realizacja(T,Fs): | ||
− | + | dt = 1.0/Fs | |
− | + | t = np.arange(0,T,dt) | |
− | return | + | s = np.sin(2*np.pi*20*t) |
+ | s += 5*np.random.randn(len(s)) | ||
+ | return s | ||
T=10.0 | T=10.0 | ||
Fs = 100.0 | Fs = 100.0 | ||
N = T*Fs | N = T*Fs | ||
− | + | ||
okno = np.blackman(N) | okno = np.blackman(N) | ||
− | + | ile_okien = 10 | |
+ | Nw = N/ile_okien | ||
+ | okno_w = np.blackman(Nw) | ||
N_rep = 100 | N_rep = 100 | ||
− | S_perio = np.zeros((N_rep,N)) | + | S_perio = np.zeros((N_rep,N/2+1)) |
− | S_welch = np.zeros((N_rep,Nw)) | + | S_welch = np.zeros((N_rep,Nw/2+1)) |
− | + | ||
for i in range(N_rep): | for i in range(N_rep): | ||
s = realizacja(T,Fs) | s = realizacja(T,Fs) | ||
− | (P | + | (F_p, P) = periodogram(s,okno,Fs) |
− | S_perio[i] = P | + | S_perio[i,:] = P |
− | + | (F_w, P_w) = welch(s,Fs,window = okno_w, nperseg = Nw, noverlap = Nw/2,scaling = 'density', return_onesided = True) | |
− | + | S_welch[i,:] = P_w * ile_okien | |
− | |||
− | S_welch[i] = | ||
− | |||
− | |||
− | py.figure( | + | py.figure(1) |
− | py.subplot( | + | py.subplot(3,1,1) |
− | py.plot(np. | + | 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( | + | py.subplot(3,1,3) |
− | py.plot(np.std(S_welch,axis = 0)/np.mean(S_welch,axis = 0)) | + | py.plot(F_p, np.std(S_perio,axis = 0)/np.mean(S_perio,axis = 0)*100) |
+ | py.plot(F_w, np.std(S_welch,axis = 0)/np.mean(S_welch,axis = 0)*100) | ||
+ | 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() | py.show() | ||
Linia 680: | Linia 887: | ||
--> | --> | ||
− | + | === Zadanie 7: Estymacja widma mocy metodą wielookienkową Thomsona (multitaper) === | |
− | Proszę napisać funkcję do estymacji mocy metodą | + | Jeśli nie mamy do dyspozycji dostatecznie długiego sygnału stacjonarnego i ergodycznego aby astosować metodę Welcha pomocne może być wykorzystanie zestawów okien ortogonalnych (Discrete Prolate Spheroidal Sequences- DPSS). Ponieważ są ortogonalne więc widma estymowane periodogramem z każdym z tych okienek stanowią niezależne estymaty gęstości mocy. Ich wartość średnia ma mniejszą wariancję niż estymata za pomocą pojedynczego periodogramu. Oczywiście nie ma nic za darmo: za mniejszą wariancję płacimy szerokością piku. |
− | Funkcja powinna pobierać następujące argumenty: sygnał, iloczyn NW, częstość próbkowania sygnału. Funkcja powinna zwracać krotkę <tt>( | + | |
+ | Metoda ta została zaproponowana w pracy: D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proceedings of the IEEE, vol. 70, no. 9, pp. 1055 – 1096, 1982 | ||
+ | |||
+ | ==== Zestawy okien ortogonalnych ==== | ||
+ | Najpierw zobaczmy jak wyglądają sekwencje okien. | ||
+ | * Moduł zawierający funkcję do generowania takiej sekwencji można ściągnąć stąd: http://www.fuw.edu.pl/~jarekz/dpss.py | ||
+ | <source lang = python> | ||
+ | # -*- coding: utf-8 -*- | ||
+ | import pylab as py | ||
+ | import numpy as np | ||
+ | from dpss import dpss_window | ||
+ | |||
+ | NW = 2 # szerokość pasma w którym zlokalizowane są piki główne okien | ||
+ | K = 2*NW-1 # liczba okien | ||
+ | N = 100 # rozmiar okna | ||
+ | py.figure(1) | ||
+ | w, eigen = dpss_window(N, NW, K) # generujemy okna | ||
+ | for i, eig in enumerate(eigen): | ||
+ | py.plot(w[i,:]) # kolejno wykreślamy wszystkie okna | ||
+ | py.legend(range(len(eigen))) | ||
+ | py.show() | ||
+ | |||
+ | print(eigen) | ||
+ | # sprawdzamy czy okna są ortogonalne | ||
+ | print('Iloczyny skalarne sekwencji okien:') | ||
+ | for i in range(len(eigen)): | ||
+ | for j in range(i,len(eigen)): | ||
+ | print('okno '+str(i)+' z oknem '+str(j)+': '+'{:.5f}'.format( np.dot(w[i,:],w[j,:]) ) ) | ||
+ | </source> | ||
+ | |||
+ | <!-- | ||
+ | # -*- coding: utf-8 -*- | ||
+ | import pylab as py | ||
+ | import numpy as np | ||
+ | from dpss import dpss_window | ||
+ | |||
+ | NW = 2 | ||
+ | K = 2*NW-1 | ||
+ | N = 100 | ||
+ | py.figure(1) | ||
+ | w, eigen = dpss_window(N, NW, 2*NW-1) | ||
+ | for i, eig in enumerate(eigen): | ||
+ | py.plot(w[i,:]) | ||
+ | py.legend(range(len(eigen))) | ||
+ | py.show() | ||
+ | |||
+ | print(eigen) | ||
+ | print('Iloczyny skalarne sekwencji okien:') | ||
+ | for i in range(len(eigen)): | ||
+ | for j in range(i,len(eigen)): | ||
+ | print('okno '+str(i)+' z oknem '+str(j)+': '+'{:.5f}'.format( np.dot(w[i,:],w[j,:]) ) ) | ||
+ | --> | ||
+ | |||
+ | ==== Estymacja widma mocy ==== | ||
+ | Proszę napisać funkcję do estymacji mocy metodą wielookienkową. | ||
+ | |||
+ | Funkcja powinna pobierać następujące argumenty: sygnał, iloczyn NW, częstość próbkowania sygnału. Funkcja powinna zwracać krotkę <tt>(F,P)</tt> gdzie <tt>P</tt> widmo mocy, <tt>F</tt> skala częstości. | ||
Przykładowe wywołanie takiej funkcji powinno wyglądać tak: | Przykładowe wywołanie takiej funkcji powinno wyglądać tak: | ||
<tt> (S,F) = mtm(s, NW = 3, Fs = 128)</tt> | <tt> (S,F) = mtm(s, NW = 3, Fs = 128)</tt> | ||
+ | |||
+ | 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 [[Nieparametryczne_widmo_mocy#Okienkowanie_a_widmo_mocy:_periodogram|periodogram]] z oknem prostokątnym. | ||
+ | |||
+ | |||
Algorytm do zastosowania wewnątrz funkcji: | Algorytm do zastosowania wewnątrz funkcji: | ||
# Oblicz maksymalną liczbę okienek <tt> K = 2*NW-1</tt> | # Oblicz maksymalną liczbę okienek <tt> K = 2*NW-1</tt> | ||
# Oblicz długość sygnału | # Oblicz długość sygnału | ||
# wygeneruj serię okienek dpss | # 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: <tt>Si = np.abs(fft(s*w | + | # dla każdego z otrzymanych okienek oblicz widmo mocy iloczynu tego okienka i sygnału. Dla i-tego okienka będzie to: <tt>Si = np.abs(fft(s*w[i]))**2</tt> |
# uśrednij widma otrzymane dla wszystkich okienek | # uśrednij widma otrzymane dla wszystkich okienek | ||
# wygeneruj oś częstości (<tt>fftfreq</tt>) | # wygeneruj oś częstości (<tt>fftfreq</tt>) | ||
− | + | Uzupełnij poniższy kod: | |
− | + | ||
+ | <source lang = python> | ||
+ | # -*- coding: utf-8 -*- | ||
+ | import pylab as py | ||
+ | import numpy as np | ||
+ | from dpss import dpss_window | ||
+ | from numpy.fft import rfft,rfftfreq | ||
+ | |||
+ | from matplotlib.font_manager import FontProperties | ||
+ | font = FontProperties('Arial') | ||
+ | |||
+ | |||
+ | 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 | ||
+ | ''' | ||
+ | ... | ||
+ | return (s,t) | ||
+ | |||
+ | def periodogram(s, okno , Fs): | ||
+ | '''peiodogram sygnału s | ||
+ | okno - synał będzie przez nie przemnożony w czasie | ||
+ | F_samp- częstość próbkowania''' | ||
+ | ... | ||
+ | return (F,P) | ||
+ | |||
+ | def mtm(s, NW , Fs): | ||
+ | '''estymacja widma w oparciu o metodę Multiteper | ||
+ | D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proceedings of the | ||
+ | IEEE, vol. 70, no. 9, pp. 1055 – 1096, 1982. | ||
+ | x - sygnał | ||
+ | N -ilość punktów okna | ||
+ | NW - iloczyn długości okna w czasie i szerokości w częstości | ||
+ | K - ilość okien | ||
+ | |||
+ | funkcja zwraca estymatę mocy widmowej | ||
+ | ''' | ||
+ | K = 2*NW-1 | ||
+ | N = len(s) | ||
+ | w, eigen = ...# wygeneruj sekwencję okien DPSS | ||
+ | P_tmp =0 | ||
+ | for i in range(K): #dla każdego okna | ||
+ | (F,Pi)= ...# oblicz periodogram | ||
+ | P_tmp+= ...# dodaj do zmiennej tymczasowej | ||
+ | P = ...# moc jest średnią z periodogramów dla poszczególnych okien | ||
+ | F = rfftfreq(N,1.0/Fs) | ||
+ | return (F, P) | ||
+ | |||
+ | #prezentacja widma | ||
+ | Fs = 200.0 # częstość próbkowania | ||
+ | |||
+ | # tworzymy sygnał testowy | ||
+ | (s1,t) = sin(f=10.2,Fs=Fs) | ||
+ | (s2,t) = sin(f=17.2,Fs=Fs) | ||
+ | s = s1+s2+np.random.randn(len(t)) | ||
+ | |||
+ | py.figure(1) | ||
+ | NW = 2 # ustalamy szerokość pasma | ||
+ | (F_m,P_m) = ... # estymujemy widmo metodą mtm | ||
+ | (F_p,P_p) = ... # estymujemy widmo metodą periodogram z oknem prostokątnym | ||
+ | # wykreślamy wyniki | ||
+ | py.plot(F_m,P_m) | ||
+ | py.plot(F_p,P_p ,'g') | ||
+ | |||
+ | # opisy wykresu | ||
+ | py.xlabel('Częstość [Hz]', fontproperties = font) | ||
+ | py.ylabel('Gestość mocy V^2/Hz', fontproperties = font) | ||
+ | py.title('Porównanie estymat gęstości mocy: wielokoienkowej i periodogramu', fontproperties = font) | ||
+ | py.legend(('wielookienkowa','periodogram')) | ||
+ | |||
+ | # test zacowania energii | ||
+ | print('Test zachowania energii:') | ||
+ | print( 'energia w czasie = \t\t'+ '{:.5f}'.format( ... )) | ||
+ | print( 'energia w mtm = \t\t'+ '{:.5f}'.format( ... )) | ||
+ | print( 'energia w periodogramie = \t'+ '{:.5f}'.format( ... )) | ||
+ | py.show() | ||
+ | |||
+ | </source> | ||
+ | <!-- | ||
+ | # -*- coding: utf-8 -*- | ||
+ | import pylab as py | ||
+ | import numpy as np | ||
+ | from dpss import dpss_window | ||
+ | from numpy.fft import rfft,rfftfreq | ||
+ | |||
+ | from matplotlib.font_manager import FontProperties | ||
+ | font = FontProperties('Arial') | ||
+ | |||
+ | |||
+ | 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 , Fs): | ||
+ | '''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/Fs) | ||
+ | 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 mtm(s, NW , Fs): | ||
+ | '''estymacja widma w oparciu o metodę Multiteper | ||
+ | D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proceedings of the | ||
+ | IEEE, vol. 70, no. 9, pp. 1055 – 1096, 1982. | ||
+ | x - sygnał | ||
+ | N -ilość punktów okna | ||
+ | NW - iloczyn długości okna w czasie i szerokości w częstości | ||
+ | K - ilość okien | ||
+ | |||
+ | funkcja zwraca estymatę mocy widmowej | ||
+ | ''' | ||
+ | K = 2*NW-1 | ||
+ | N = len(s) | ||
+ | w, eigen = dpss_window(N, NW, K) | ||
+ | P=0 | ||
+ | for i in range(K): | ||
+ | (F,Pi)=periodogram(s, w[i,:] , Fs) | ||
+ | P+=Pi | ||
+ | P = P/K | ||
+ | F = rfftfreq(N,1.0/Fs) | ||
+ | return (F, P) | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | #prezentacja widma | ||
+ | Fs = 200.0 | ||
+ | |||
+ | (s1,t) = sin(f=10.2,Fs=Fs) | ||
+ | (s2,t) = sin(f=17.2,Fs=Fs) | ||
+ | s = s1+s2+np.random.randn(len(t)) | ||
+ | |||
+ | py.figure(1) | ||
+ | NW = 2 | ||
+ | (F_m,P_m) = mtm(s, NW, Fs) | ||
+ | (F_p,P_p) = periodogram(s, np.ones(len(t)), Fs) | ||
+ | py.plot(F_m,P_m) | ||
+ | py.plot(F_p,P_p ,'g') | ||
+ | py.xlabel('Częstość [Hz]', fontproperties = font) | ||
+ | py.ylabel('Gestość mocy V^2/Hz', fontproperties = font) | ||
+ | py.title('Porównanie estymat gęstości mocy: wielokoienkowej i periodogramu', fontproperties = font) | ||
+ | py.legend(('wielookienkowa','periodogram')) | ||
+ | # test zacowania energii | ||
+ | print('Test zachowania energii:') | ||
+ | print( 'energia w czasie = \t\t'+ '{:.5f}'.format(np.sum(s**2*1/Fs))) | ||
+ | print( 'energia w mtm = \t\t'+ '{:.5f}'.format(np.sum(P_m))) | ||
+ | print( 'energia w periodogramie = \t'+ '{:.5f}'.format(np.sum(P_p))) | ||
+ | py.show() | ||
+ | --> | ||
+ | |||
+ | ===Zadanie 8 === | ||
+ | Proszę wykonać ilustrację średniej wraz z przedziałami ufności 95% oraz błędu względnego w estymatorze wielookienkowym (dla porównania periodogram), analogicznie jak w zadaniu 6. | ||
+ | |||
<!-- | <!-- | ||
− | + | # -*- coding: utf-8 -*- | |
import pylab as py | import pylab as py | ||
import numpy as np | import numpy as np | ||
− | import | + | from dpss import dpss_window |
− | from numpy.fft import | + | from numpy.fft import rfft,rfftfreq |
+ | import scipy.stats as st | ||
+ | from matplotlib.font_manager import FontProperties | ||
+ | font = FontProperties('Arial') | ||
+ | |||
def sin(f = 1, T = 1, Fs = 128, phi =0 ): | 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 mtm(s, NW | + | def mtm(s, NW , Fs = 128): |
'''estymacja widma w oparciu o metodę Multiteper | '''estymacja widma w oparciu o metodę Multiteper | ||
D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proceedings of the | D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proceedings of the | ||
Linia 727: | Linia 1166: | ||
K = 2*NW-1 | K = 2*NW-1 | ||
N = len(s) | N = len(s) | ||
− | w = | + | w, eigen = dpss_window(N, NW, K) |
− | + | P=0 | |
for i in range(K): | for i in range(K): | ||
− | + | (F,Pi)=periodogram(s, w[i,:] , Fs) | |
− | + | P+=Pi | |
− | + | P = P/K | |
− | F = | + | F = rfftfreq(N,1.0/Fs) |
− | return ( | + | return (F, P) |
+ | |||
+ | |||
+ | 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) | ||
+ | s = np.sin(2*np.pi*20*t) | ||
+ | s += 5*np.random.randn(len(s)) | ||
+ | return s | ||
+ | |||
+ | T=2.0 | ||
Fs = 100.0 | Fs = 100.0 | ||
− | + | N = T*Fs | |
− | + | NW = 2 | |
− | + | N_rep = 1000 | |
− | ( | + | S_perio = np.zeros((N_rep, N//2+1)) |
− | + | S_mtm = np.zeros((N_rep, N//2+1)) | |
+ | for i in range(N_rep): | ||
+ | s = realizacja(T,Fs) | ||
+ | (F_p, P) = periodogram(s,np.hamming(len(s)),Fs) | ||
+ | S_perio[i,:] = P | ||
+ | (F_m, P_w) = mtm(s, NW = NW, Fs = Fs) | ||
+ | S_mtm[i,:] = P_w | ||
− | py.plot( | + | 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_m,np.mean(S_mtm,axis = 0),'r', | ||
+ | F_m, st.scoreatpercentile(S_mtm, 2.5,axis =0),'b', | ||
+ | F_m, st.scoreatpercentile(S_mtm,97.5,axis =0),'b' ) | ||
+ | py.title('MTM: średnie widmo realizacji wraz z 95% CI', fontproperties = font) | ||
+ | py.ylabel('V^2/Hz') | ||
− | + | py.subplot(3,1,3) | |
− | + | py.plot(F_p, np.std(S_perio,axis = 0)/np.mean(S_perio,axis = 0)*100) | |
+ | py.plot(F_m, np.std(S_mtm,axis = 0)/np.mean(S_mtm,axis = 0)*100) | ||
+ | py.ylim([0,150]) | ||
+ | py.xlabel('Częstość Hz', fontproperties = font) | ||
+ | py.ylabel('%') | ||
+ | py.legend(('periodogram','MTM')) | ||
+ | py.title('Błąd względny', fontproperties = font) | ||
py.show() | py.show() | ||
− | + | --> | |
− | + | ==Co musimy z tego zapamiętać == | |
− | + | * Jak definiujemy moc sygnału i energię w dziedzinie czasu w analizie sygnałów? | |
− | - | + | * Jak definiujemy gęstość energii i energię sygnału w dziedzinie częstości? |
+ | * Jak estymować periodogram? | ||
+ | * Co to znaczy że sygnał jest stochastyczny? | ||
+ | * Co to znaczy że sygnał jest stacjonarny i ergodyczny? | ||
+ | * Jaki jest błąd względny widma białego szumu estymowanego za pomocą periodogramu? | ||
+ | * Metody zmniejszenia błędu względnego: metoda Welcha i metoda wielookienkowa Thomsona - na czym polegają, jakie są podobieństwa i różniece w stosowaniu tych metod? | ||
+ | * W jakich sytuacjach wybrać metodę Welcha a w jakich Thomsona? | ||
[[Analiza_sygnałów_-_ćwiczenia]]/Fourier_4 | [[Analiza_sygnałów_-_ćwiczenia]]/Fourier_4 |
Aktualna wersja na dzień 14:14, 16 lis 2016
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ą wielookienkową Thomsona (multitaper)
- 2.7 Zadanie 8
- 3 Co musimy z tego zapamiętać
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:
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:
ś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 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 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
okno = ...
(moc_w_czestosci, F) = periodogram(syg, Fs=FS,okno = okno)
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ą wielookienkową Thomsona (multitaper)
Jeśli nie mamy do dyspozycji dostatecznie długiego sygnału stacjonarnego i ergodycznego aby astosować metodę Welcha pomocne może być wykorzystanie zestawów okien ortogonalnych (Discrete Prolate Spheroidal Sequences- DPSS). Ponieważ są ortogonalne więc widma estymowane periodogramem z każdym z tych okienek stanowią niezależne estymaty gęstości mocy. Ich wartość średnia ma mniejszą wariancję niż estymata za pomocą pojedynczego periodogramu. Oczywiście nie ma nic za darmo: za mniejszą wariancję płacimy szerokością piku.
Metoda ta została zaproponowana w pracy: D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proceedings of the IEEE, vol. 70, no. 9, pp. 1055 – 1096, 1982
Zestawy okien ortogonalnych
Najpierw zobaczmy jak wyglądają sekwencje okien.
- Moduł zawierający funkcję do generowania takiej sekwencji można ściągnąć stąd: http://www.fuw.edu.pl/~jarekz/dpss.py
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
from dpss import dpss_window
NW = 2 # szerokość pasma w którym zlokalizowane są piki główne okien
K = 2*NW-1 # liczba okien
N = 100 # rozmiar okna
py.figure(1)
w, eigen = dpss_window(N, NW, K) # generujemy okna
for i, eig in enumerate(eigen):
py.plot(w[i,:]) # kolejno wykreślamy wszystkie okna
py.legend(range(len(eigen)))
py.show()
print(eigen)
# sprawdzamy czy okna są ortogonalne
print('Iloczyny skalarne sekwencji okien:')
for i in range(len(eigen)):
for j in range(i,len(eigen)):
print('okno '+str(i)+' z oknem '+str(j)+': '+'{:.5f}'.format( np.dot(w[i,:],w[j,:]) ) )
Estymacja widma mocy
Proszę napisać funkcję do estymacji mocy metodą wielookienkową.
Funkcja powinna pobierać następujące argumenty: sygnał, iloczyn NW, częstość próbkowania sygnału. Funkcja powinna zwracać krotkę (F,P) gdzie P 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)
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.
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[i]))**2
- uśrednij widma otrzymane dla wszystkich okienek
- wygeneruj oś częstości (fftfreq)
Uzupełnij poniższy kod:
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
from dpss import dpss_window
from numpy.fft import rfft,rfftfreq
from matplotlib.font_manager import FontProperties
font = FontProperties('Arial')
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
'''
...
return (s,t)
def periodogram(s, okno , Fs):
'''peiodogram sygnału s
okno - synał będzie przez nie przemnożony w czasie
F_samp- częstość próbkowania'''
...
return (F,P)
def mtm(s, NW , Fs):
'''estymacja widma w oparciu o metodę Multiteper
D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proceedings of the
IEEE, vol. 70, no. 9, pp. 1055 – 1096, 1982.
x - sygnał
N -ilość punktów okna
NW - iloczyn długości okna w czasie i szerokości w częstości
K - ilość okien
funkcja zwraca estymatę mocy widmowej
'''
K = 2*NW-1
N = len(s)
w, eigen = ...# wygeneruj sekwencję okien DPSS
P_tmp =0
for i in range(K): #dla każdego okna
(F,Pi)= ...# oblicz periodogram
P_tmp+= ...# dodaj do zmiennej tymczasowej
P = ...# moc jest średnią z periodogramów dla poszczególnych okien
F = rfftfreq(N,1.0/Fs)
return (F, P)
#prezentacja widma
Fs = 200.0 # częstość próbkowania
# tworzymy sygnał testowy
(s1,t) = sin(f=10.2,Fs=Fs)
(s2,t) = sin(f=17.2,Fs=Fs)
s = s1+s2+np.random.randn(len(t))
py.figure(1)
NW = 2 # ustalamy szerokość pasma
(F_m,P_m) = ... # estymujemy widmo metodą mtm
(F_p,P_p) = ... # estymujemy widmo metodą periodogram z oknem prostokątnym
# wykreślamy wyniki
py.plot(F_m,P_m)
py.plot(F_p,P_p ,'g')
# opisy wykresu
py.xlabel('Częstość [Hz]', fontproperties = font)
py.ylabel('Gestość mocy V^2/Hz', fontproperties = font)
py.title('Porównanie estymat gęstości mocy: wielokoienkowej i periodogramu', fontproperties = font)
py.legend(('wielookienkowa','periodogram'))
# test zacowania energii
print('Test zachowania energii:')
print( 'energia w czasie = \t\t'+ '{:.5f}'.format( ... ))
print( 'energia w mtm = \t\t'+ '{:.5f}'.format( ... ))
print( 'energia w periodogramie = \t'+ '{:.5f}'.format( ... ))
py.show()
Zadanie 8
Proszę wykonać ilustrację średniej wraz z przedziałami ufności 95% oraz błędu względnego w estymatorze wielookienkowym (dla porównania periodogram), analogicznie jak w zadaniu 6.
Co musimy z tego zapamiętać
- Jak definiujemy moc sygnału i energię w dziedzinie czasu w analizie sygnałów?
- Jak definiujemy gęstość energii i energię sygnału w dziedzinie częstości?
- Jak estymować periodogram?
- Co to znaczy że sygnał jest stochastyczny?
- Co to znaczy że sygnał jest stacjonarny i ergodyczny?
- Jaki jest błąd względny widma białego szumu estymowanego za pomocą periodogramu?
- Metody zmniejszenia błędu względnego: metoda Welcha i metoda wielookienkowa Thomsona - na czym polegają, jakie są podobieństwa i różniece w stosowaniu tych metod?
- W jakich sytuacjach wybrać metodę Welcha a w jakich Thomsona?
Analiza_sygnałów_-_ćwiczenia/Fourier_4