Nieparametryczne widmo mocy: Różnice pomiędzy wersjami

Z Brain-wiki
m
 
(Nie pokazano 64 wersji utworzonych przez 2 użytkowników)
Linia 2: Linia 2:
  
 
==Widmo mocy==
 
==Widmo mocy==
 +
<!--
 
=== Moc===
 
=== Moc===
 
Moc chwilowa sygnału przez analogię do układów elektrycznych o jednostkowym oporze jest w analizie sygnałów przyjęta jako kwadraty próbek (<math>P = I^2 R = \frac{U^2}{R}</math>).
 
Moc chwilowa sygnału przez analogię do układów elektrycznych o jednostkowym oporze jest w analizie sygnałów przyjęta jako kwadraty próbek (<math>P = I^2 R = \frac{U^2}{R}</math>).
Linia 9: Linia 10:
 
a energia wzorem:
 
a energia wzorem:
 
:<math>E = \sum _n{x[n]^2}</math>
 
:<math>E = \sum _n{x[n]^2}</math>
==== Zadanie ====
 
Przy pomocy funkcji napisanych na poprzednich zajęciach proszę wygenerować sygnał sinusoidalny o amplitudzie 1, częstości 10 Hz, trwający 0.3 sekundy i próbkowany z częstością 1000 Hz. Proszę narysować ten sygnał przy pomocy funkcji <tt>pylab.stem</tt>, obliczyć i narysować przebieg mocy w czasie, obliczyć energię tego sygnału.
 
  
 
===Widmo mocy: tw. Plancherela i tw. Parsevala ===
 
===Widmo mocy: tw. Plancherela i tw. Parsevala ===
Linia 54: Linia 53:
 
czyli
 
czyli
 
:<math>\sum_{k=0}^{N - 1} |x[k]|^2  = \frac{1}{N} \sum_{r=0}^{N - 1} |X[r]|^2</math>
 
:<math>\sum_{k=0}^{N - 1} |x[k]|^2  = \frac{1}{N} \sum_{r=0}^{N - 1} |X[r]|^2</math>
 +
-->
 +
===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 <tt>pylab.stem</tt>,
 +
* 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: funkcja do estymacji widma mocy ===
+
==== Zadanie 2: Moc i energia sygnału w dziedzinie czasu i częstości ====
Proszę napisać i przetestować funkcję realizującą następujący algorytm estymacji widma mocy:
 
# Jako parametry funkcja przyjmuje sygnał <tt>s</tt> i częstość próbkowania tego sygnału
 
# Oblicz transformatę Fouriera sygnału przy pomocy <tt>fft</tt>
 
# Unormuj otrzymaną transformatę przez pierwiastek z ilości próbek sygnału
 
# Oblicz moc jako iloczyn unormowanej transformaty i jej sprzężenia zespolonego. W pythonie sprzężenie zespolone liczby <tt>S</tt> otrzymujemy korzystając z metody <tt>S.conj()</tt>. W wyniku powyższego iloczynu dostaniemy liczby zespolone o zerowej części urojonej (proszę to sprawdzić).
 
# Do dalszych operacji wybierz tylko część rzeczywistą mocy. Część rzeczywistą liczby zespolonej M pobieramy w następujący sposób: <tt>czesc_rzeczywista_M = M.real</tt>
 
# Korzystając z funkcji <tt>fftfreq</tt> obliczamy częstości, dla których policzone są współczynniki Fouriera.
 
# Przy pomocy funkcji <tt>fftshift</tt> porządkujemy kolejność wektorów mocy i częstości, tak aby częstości były reprezentowane od -częstości Nyquista, przez 0, do +częstości Nyquista
 
# zwracamy uporządkowane wektory mocy i częstości
 
  
Testowanie:
+
* 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
  
# przy pomocy funkcji sin napisanej na drugich ćwiczeniach  wygeneruj sinusoidę o długości 1 s, o częstości 10Hz, próbkowaną 100 Hz.
+
<source lang = python>
# oblicz moc w czasie (trzeba podnieść wartość każdej próbki do kwadratu)
+
#!/usr/bin/env python3
# oblicz moc w częstości przy pomocy funkcji estymacji widma mocy
+
# -*- coding: utf-8 -*-
# oblicz energię sygnału w czasie
 
# oblicz energię sygnału w częstości
 
# wykreśl:  
 
#* sygnał
 
#* przebieg jego mocy w czasie
 
#* przebieg mocy w częstości
 
  
*
 
<!--
 
<source lang =python>
 
# -*- coding: utf-8 -*-
 
 
 
import pylab as py
 
import pylab as py
 
import numpy as np
 
import numpy as np
from numpy.fft import rfft,fft,fftfreq,fftshift
+
from numpy.fft import rfft,rfftfreq
 
   
 
   
def widmo_mocy(s, Fs):
+
def widmo_mocy(s,Fs,okno):
     S = fft(s)/np.sqrt(len(s))
+
    okno= # znormalizuj okno
     S_moc = S*S.conj()
+
    s = # zokienkuj sygnał
     S_moc = S_moc.real
+
     S = # Oblicz transformatę Fouriera sygnału przy pomocy funkcji <tt>rfft</tt>
     F = fftfreq(len(s), 1/Fs)
+
    P =  # Oblicz moc jako iloczyn unormowanej transformaty i jej sprzężenia zespolonego.
    return (fftshift(S_moc),fftshift(F))
+
    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 )
  
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
+
# okno prostokątne
    '''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania
+
okno = np.ones(len(s))
    Domyślnie wytwarzany jest sygnał reprezentujący
+
moc_w_czasie = ...
    1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
+
(F, moc_w_czestosci) = widmo_mocy(s, Fs=fs, okno = okno)
    '''
 
    dt = 1.0/Fs
 
    t = np.arange(0,T,dt)
 
    s = np.sin(2*np.pi*f*t + phi)
 
    return (s,t)
 
 
 
(s,t) = sin(f=10,Fs=100.0)
 
moc_w_czasie = s**2
 
(moc_w_czestosci, F) = widmo_mocy(s, Fs=100.0)
 
  
energia_w_czasie = np.sum(moc_w_czasie)
+
dt = 1/fs
energia_w_czestosci = np.sum(moc_w_czestosci)
+
energia_w_czasie = ...
print 'energia w czasie: ', energia_w_czasie
+
energia_w_czestosci = ...
print 'energia w czestosci: ', energia_w_czestosci
 
  
 
py.subplot(3,1,1)
 
py.subplot(3,1,1)
 
py.plot(t,s)
 
py.plot(t,s)
py.title('Sygnał')
+
py.title('Sygnal')
 
py.subplot(3,1,2)
 
py.subplot(3,1,2)
 
py.plot(t,moc_w_czasie)
 
py.plot(t,moc_w_czasie)
py.title('moc w czasie')
+
py.title('moc w czasie, energia w czasie: ' +str(energia_w_czasie))
 
py.subplot(3,1,3)
 
py.subplot(3,1,3)
 
py.plot(F,moc_w_czestosci)
 
py.plot(F,moc_w_czestosci)
py.title('moc w częstości')
+
py.title('moc w czestosci, energia w czestosci: ' +str(energia_w_czestosci))
 
py.show()
 
py.show()
 
</source>
 
</source>
-->
 
  
 
<!--
 
<!--
===Rozdzielczość widma===
+
#!/usr/bin/env python3
Proszę zbadać rozdzielczość częstotliwościową widma w zależności od długości sygnału. Używamy funkcji <tt>m_sin</tt>. Zwrócić uwagę na to co dzieje się na końcach odcinka.
+
# -*- coding: utf-8 -*-
-->
+
"""
 +
Created on Mon Nov  7 12:21:33 2016
 +
 
 +
@author: admin
 +
"""
 +
 
 +
import pylab as py
 +
import numpy as np
 +
from numpy.fft import rfft,rfftfreq
 +
 +
def widmo_mocy(s,Fs,okno):
 +
    okno= okno/np.linalg.norm(okno)
 +
    s = s*okno
 +
    S = rfft(s)
 +
    P = S*S.conj()
 +
    P = P/Fs
 +
    P = P.real
 +
    if len(s)%2 ==0:
 +
        P[1:-1] *=2
 +
    else:
 +
        P[1:] *=2
 +
    F = rfftfreq(len(s), 1./Fs)
 +
    return F,P
  
===Okienkowanie a widmo mocy: periodogram===
+
# część testująca
Przypomnijmy wzór na dyskretną transformatę Fouriera [http://haar.zfb.fuw.edu.pl/edu/index.php/%C4%86wiczenia_2 DFT] zaimplementowaną w FFT:
+
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 )
  
:<math>S[k] = \sum_{n=0}^{n-1} s[n] \exp\left\{-2\pi i{nk \over N}\right\}      \qquad k = 0,\ldots,N-1.
+
# okno prostokątne
</math>
+
okno = np.ones(len(s))
  
 +
moc_w_czasie = s**2
 +
(F, moc_w_czestosci) = widmo_mocy(s, Fs=fs,okno=okno)
 +
dt = 1/fs
  
Na podstawie twierdzenia [[Nieparametryczne_widmo_mocy#Twierdzenie_Parsevala|Parsevala]] możemy policzyć widmo mocy jako:
+
energia_w_czasie = np.sum(moc_w_czasie*dt)
<math>
+
energia_w_czestosci = np.sum(moc_w_czestosci)
P[k] = \frac{1}{N}  \left|S[k]\right|^2
 
</math>
 
  
Jeśli do liczenia mocy chcielibyśmy posłużyć się techniką okiennkowania sygnału, to powinniśmy używać okienek znormalizowanych, czyli takich których energia jest równa 1, wtedy mnożenie przez okienko nie zaburzy estymaty energii sygnału.  
+
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:
 
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 151: 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  
 
</math>
 
</math>
  
====Zadanie ====
+
====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.
 
*
 
 
<!--
 
<!--
 
<source lang = python>
 
<source lang = python>
Linia 168: Linia 206:
 
import pylab as py
 
import pylab as py
 
import numpy as np
 
import numpy as np
from numpy.fft import  fft, fftfreq, fftshift
+
from numpy.fft import  rfft, rfftfreq
 
   
 
   
 
   
 
   
Linia 183: Linia 221:
  
 
def widmo_dB(s, N_fft , F_samp):
 
def widmo_dB(s, N_fft , F_samp):
     S = fft(s,N_fft)/np.sqrt(N_fft)
+
     S = rfft(s,N_fft)/np.sqrt(N_fft)
 
     S_dB = 20*np.log10(np.abs(S))
 
     S_dB = 20*np.log10(np.abs(S))
     F = fftfreq(N_fft, 1/F_samp)
+
     F = rfftfreq(N_fft, 1/F_samp)
     return (fftshift(S_dB),fftshift(F))
+
     return S_dB,F
  
 
def periodogram(s, okno , F_samp):
 
def periodogram(s, okno , F_samp):
Linia 195: Linia 233:
 
     s = s*okno
 
     s = s*okno
 
     N_fft = len(s)
 
     N_fft = len(s)
     S = fft(s,N_fft)#/np.sqrt(N_fft)
+
     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 wartośći urojone, ale trzeba ykonać konwersję typów
+
     P = P.real # P i tak ma zerowe wartości urojone, ale trzeba ykonać konwersję typów
     F = fftfreq(N_fft, 1/F_samp)
+
     F = rfftfreq(N_fft, 1/F_samp)
     return (fftshift(P),fftshift(F))
+
    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 241: 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 k†órego pochodzą. Często w opisie takich zmiennych posługujemy się [[WnioskowanieStatystyczne/Zmienne_losowe_i_generatory_liczb_pseudolosowych#Momenty|momentami rozkładów]].
+
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 252: 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ć prze wartości oczekiwane liczone po realizacjach.
+
 
 +
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 260: 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 wartośći zaobserwowanych w chwili we wszystkich realizacjach <math>t_1</math> ważona prawdopodobieństwem wystąpienia tej realizacji:
+
ś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 273: Linia 347:
 
</dl>
 
</dl>
  
Założenie o sygnale, że jest stacjonarny i ergodyczny pozwala zamienić sumowanie po realizacjach na sumowanie po czasie w estymatory momentów statystycznych.
+
Założenie o sygnale, że jest stacjonarny i ergodyczny pozwala zamienić sumowanie po realizacjach na sumowanie po czasie w estymatorach momentów statystycznych.
  
===Transformata Fouriera sygnału stochastycznego===
+
===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 288: 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 298: Linia 465:
 
import scipy.stats as st
 
import scipy.stats as st
  
from numpy.fft import rfft,fft,fftfreq,fftshift
+
from numpy.fft import rfft,rfftfreq
 
   
 
   
def widmo_mocy(s, Fs):
 
    '''funkcja oblicza widmo mocy sygnału i oś częstości
 
    s - sygnał
 
    Fs - częstość próbkowania
 
    '''
 
    S = fft(s)/np.sqrt(len(s))
 
    S_moc = S*S.conj()
 
    S_moc = S_moc.real
 
    F = fftfreq(len(s), 1/Fs)
 
    return (fftshift(S_moc),fftshift(F))
 
  
def widmo_mocy_rzeczywistego(s,FS):
+
 
 +
def widmo_mocy(s,Fs):
 
     '''funkcja oblicza widmo mocy sygnału rzeczywistego i oś częstości
 
     '''funkcja oblicza widmo mocy sygnału rzeczywistego i oś częstości
 
     s - sygnał
 
     s - sygnał
Linia 321: Linia 479:
 
     S_moc = np.abs(S)**2
 
     S_moc = np.abs(S)**2
 
     S_moc = S_moc.real
 
     S_moc = S_moc.real
     F = fftfreq(len(s), 1/Fs)
+
     F = rfftfreq(len(s), 1/Fs)
 
     return (S_moc,F)
 
     return (S_moc,F)
  
Linia 332: 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,t)
+
     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,t)
+
     return (t, s)
  
 
def dwadziescia_realizacji(FS):
 
def dwadziescia_realizacji(FS):
Linia 346: Linia 504:
 
     '''
 
     '''
 
     for i in range(20):
 
     for i in range(20):
         (s,t) = sin(f=10,Fs=FS)
+
         (t, s) = sin(f=10,Fs=FS)
         (sz,t) = szum(Fs =FS)
+
         (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 357: Linia 515:
 
     #  Skonstruuj funkcję rysującą średnie widmo wraz z 95% przedziałem ufności.  
 
     #  Skonstruuj funkcję rysującą średnie widmo wraz z 95% przedziałem ufności.  
 
     '''
 
     '''
     zbior_widm = np.zeros((20,FS))
+
     zbior_widm = np.zeros((20,FS/2+1))
 
     for i in range(20):
 
     for i in range(20):
         (s,t) = sin(f=10,Fs=FS)
+
         (t, s) = sin(f=10,Fs=FS)
         (sz,t) = szum(Fs =FS)
+
         (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 376: Linia 534:
  
 
FS =100.0     
 
FS =100.0     
#dwadziescia_realizacji(FS)
+
dwadziescia_realizacji(FS)
 
srednie_widmo(FS)
 
srednie_widmo(FS)
 
 
 
</source>
 
</source>
  
Linia 385: Linia 541:
  
 
=== Oszacowanie błędu transformaty Fouriera dla białego szumu ===
 
=== 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ść:
+
* 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 zatem rozkładowi <math>\chi^2</math> o dwóch stopniach swobody.
+
* Wielkość:
Możemy oszacować względny błąd <math>P(f_1) </math> dla danej częstości <math>f_1</math>:
+
: <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.  
:<math>\epsilon_r= \sigma_{P_{f_1}}/\mu_{P_{f_1}}</math>
+
* Wielkość ta podlega rozkładowi <math>\chi^2</math> o dwóch stopniach swobody.
 
 
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>
+
* 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%.
  
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>.
+
* 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>.  
  
Zauważmy, że w obu metodach zmniejszamy wariancję estymatora kosztem rozdzielczości w częstości.
+
:*;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>.
  
===Metoda Welcha===
+
'''Zauważmy, że w obu metodach zmniejszamy wariancję estymatora kosztem rozdzielczości w częstości.'''
Proszę napisać i przetestować funkcję implementującą metodę Welcha estymacji widma mocy. Algorytm Welcha:
 
# sygnał <math>x</math> o długości ''N'' jest dzielony na segmenty, każdy o długości <math>N_s</math>. Odcinki mogą na siebie zachodzić na <math>N_z</math> punktów. Czyli są względem siebie przesunięte o <math>N_p = N_s-N_z</math>.
 
# z każdego segmentu liczony jest okienkowany periodogram
 
# periodogramy są sumowane
 
# wynik dzielony jest przez efektywne wykorzystanie każdego kawałka sygnału w estymacie: <tt>K_eff = dlogosc_okna * ilosc_okien / dlugosc_sygnalu</tt>
 
  
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 5: Metoda Welcha===
 
 
* <!--
 
<source lang = python>
 
def pwelch(s,okienko, przesuniencie, Fs):
 
    '''s - sygnał
 
    okienko - przebieg czasowy okienka
 
    przesuniecie - o ile punktów okienka są przesówane względem siebie
 
    Fs - częstość próbkowania'''
 
    N = len(s)
 
    N_s = len(okienko)
 
   
 
    start_fragmentow = np.arange(0,N-N_s+1,przesuniencie)
 
    ile_fragmentow = len(start_fragmentow)
 
    ile_przekrycia = N_s*ile_fragmentow/float(N)
 
  
    P_sredni = np.zeros(N_s)
+
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.
    for i in range(ile_fragmentow):
 
        s_fragment = s[start_fragmentow[i]:start_fragmentow[i]+N_s]
 
        (P, F) = periodogram(s_fragment,okno,F_samp)
 
        P_sredni += P
 
    return (P_sredni/ile_przekrycia,F)
 
</source>
 
  
Całość razem z kodem testującym:
 
 
<source lang = python>
 
<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  fft, fftfreq, fftshift
+
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 461: Linia 588:
 
     return (s,t)
 
     return (s,t)
  
def periodogram(s, okno , Fs):
+
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
     Fs- 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 = fft(s,N_fft)#/np.sqrt(N_fft)
+
     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 wartośći urojone, ale trzeba ykonać konwersję typów
+
     F = rfftfreq(N_fft, 1/F_samp)
     F = fftfreq(N_fft, 1/F_samp)
+
     if len(s)%2 ==0: # dokładamy moc z ujemnej części widma
     return (fftshift(P),fftshift(F))
+
        P[1:-1] *=2
def pwelch(s,okienko, przesuniencie, Fs):
+
     else:
    '''s - sygnał
+
         P[1:] *=2
    okienko - przebieg czasowy okienka
+
     return (F,P)
    przesuniecie - o ile punktów okienka są przesówane względem siebie
 
    Fs - częstość próbkowania'''
 
    N = len(s)
 
    N_s = len(okienko)
 
   
 
    start_fragmentow = np.arange(0,N-N_s+1,przesuniencie)
 
    ile_fragmentow = len(start_fragmentow)
 
    ile_przekrycia = N_s*ile_fragmentow/float(N)
 
    P_sredni = np.zeros(N_s)
 
     for i in range(ile_fragmentow):
 
         s_fragment = s[start_fragmentow[i]:start_fragmentow[i]+N_s]
 
        (P, F) = periodogram(s_fragment,okno,Fs)
 
        P_sredni += P
 
     return (P_sredni/ile_przekrycia,F)
 
  
 
      
 
      
F_samp = 100.0
+
Fs = 100.0
(x,t) = sin(f = 3.1, T =20, Fs = F_samp, phi = 0)
+
(x,t) = sin(f = 3.1, T =20, Fs = Fs, phi = 0)
N = len(x) # dłogość sygnału
+
N = len(x) # długość sygnału
  
okno = np.ones(N)#/np.sqrt(N)
 
#okno = np.blackman(N)
 
#okno = np.hamming(N)
 
  
s = x*okno # sygnał okienkowany
+
okno = np.hamming(N)
 +
okno/=np.linalg.norm(okno)
  
 
py.subplot(2,1,1)
 
py.subplot(2,1,1)
(P, F) = periodogram(x,okno,F_samp)
+
(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 510: Linia 620:
 
#periodogram
 
#periodogram
 
py.subplot(2,1,2)
 
py.subplot(2,1,2)
N_s = N/8
+
Nseg =20
#okno = np.ones(N_s)#/np.sqrt(N)
+
N_s = N/Nseg
#okno = np.blackman(N_s)
+
 
 
okno = np.hamming(N_s)
 
okno = np.hamming(N_s)
 
+
okno/=np.linalg.norm(okno)
       
+
(F, P) = welch(...)
(P, F) = pwelch(x,okno,N_s/10,F_samp)
 
 
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===
 
 
=== 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
 
#  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
 
# dla każdej realizacji oblicz widmo mocy za pomocą periodogramu okienkowanego oknem Blackmana
Linia 533: 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  fft, fftfreq, fftshift
+
from numpy.fft import  rfft, rfftfreq
+
from scipy.signal import welch
+
import scipy.stats as st
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
+
from matplotlib.font_manager import FontProperties
     '''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
+
font = FontProperties('Arial')
     1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
+
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.sin(2*np.pi*f*t + phi)
+
     # 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
     return (s,t)
+
    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 , Fs):
+
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
     Fs- 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 = fft(s,N_fft)#/np.sqrt(N_fft)
+
     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 wartośći urojone, ale trzeba ykonać konwersję typów
+
     F = rfftfreq(N_fft, 1/F_samp)
     F = fftfreq(N_fft, 1/Fs)
+
     if len(s)%2 ==0: # dokładamy moc z ujemnej części widma
     return (fftshift(P),fftshift(F))
+
        P[1:-1] *=2
def pwelch(s,okno, przesuniencie, Fs):
+
     else:
    '''s - sygnał
+
         P[1:] *=2
    okienko - przebieg czasowy okienka
+
     return (F,P)
    przesuniecie - o ile punktów okienka są przesówane względem siebie
 
    Fs - częstość próbkowania'''
 
    N = len(s)
 
    N_s = len(okno)
 
   
 
    start_fragmentow = np.arange(0,N-N_s+1,przesuniencie)
 
    ile_fragmentow = len(start_fragmentow)
 
    ile_przekrycia = N_s*ile_fragmentow/float(N)
 
    P_sredni = np.zeros(N_s)
 
     for i in range(ile_fragmentow):
 
         s_fragment = s[start_fragmentow[i]:start_fragmentow[i]+N_s]
 
        (P, F) = periodogram(s_fragment,okno,Fs)
 
        P_sredni += P
 
     return (P_sredni/ile_przekrycia,F)
 
  
 
def realizacja(T,Fs):
 
def realizacja(T,Fs):
     (x,t) = sin(f = 20.0, T = T, Fs = Fs, phi = 0)
+
     dt = 1.0/Fs
     x += 2*np.random.randn(len(x))
+
    t = np.arange(0,T,dt)
     return x
+
    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
Nw = N/10.0
+
 
 
okno = np.blackman(N)
 
okno = np.blackman(N)
okno_welch = np.blackman(Nw)
+
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))
py.figure(1)
+
 
 
for i in range(N_rep):
 
for i in range(N_rep):
 
     s = realizacja(T,Fs)
 
     s = realizacja(T,Fs)
     (P, F) = periodogram(s,okno,Fs)  
+
     (F_p, P) = periodogram(s,okno,Fs)  
     S_perio[i] = P
+
     S_perio[i,:] = P
     py.subplot(2,1,1)
+
     (F_w, P_w) = welch(s,Fs,window = okno_w, nperseg = Nw, noverlap = Nw/2,scaling = 'density', return_onesided = True)
    py.plot(F,P)
+
     S_welch[i,:] = P_w * ile_okien
    (P, F) = pwelch(s,okno_welch,Nw/10,Fs)
 
     S_welch[i] = P
 
    py.subplot(2,1,2)
 
    py.plot(F,P)
 
  
py.figure(2)
+
py.figure(1)  
py.subplot(2,1,1)
+
py.subplot(3,1,1)
py.plot(np.std(S_perio,axis = 0)/np.mean(S_perio,axis = 0))
+
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(2,1,2)
+
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 707: Linia 887:
 
-->
 
-->
  
====Estymacja widma mocy metodą multitaper ====
+
=== Zadanie 7: Estymacja widma mocy metodą wielookienkową Thomsona (multitaper) ===
Proszę napisać funkcję do estymacji mocy metodą 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.
Funkcja powinna pobierać następujące argumenty: sygnał, iloczyn NW, częstość próbkowania sygnału. Funkcja powinna zwracać krotkę <tt>(S,F)</tt> gdzie <tt>S</tt> widmo mocy, <tt>F</tt> skala częstości.
+
 
 +
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.dpssarray[i]))**2</tt>
+
# 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>)
  
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.
+
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.
 +
 
 
<!--
 
<!--
<source lang = python>
+
# -*- coding: utf-8 -*-
 
import pylab as py
 
import pylab as py
 
import numpy as np
 
import numpy as np
import gendpss as dpss
+
from dpss import dpss_window
from numpy.fft import fft,fftshift,fftfreq
+
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
+
    '''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  
+
    Domyślnie wytwarzany jest sygnał reprezentujący  
1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
+
    1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
'''
+
    '''
 
   
 
   
dt = 1.0/Fs
+
    dt = 1.0/Fs
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,t)
+
    return (s,t)
  
  
def mtm(s, NW = 3, Fs = 128):
+
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 754: Linia 1166:
 
     K = 2*NW-1
 
     K = 2*NW-1
 
     N = len(s)
 
     N = len(s)
     w = dpss.gendpss(N,NW,K)
+
     w, eigen = dpss_window(N, NW, K)
     S=np.zeros(N)
+
     P=0
 
     for i in range(K):
 
     for i in range(K):
         Si = np.abs(fft(s*w.dpssarray[i]))**2
+
         (F,Pi)=periodogram(s, w[i,:] , Fs)
         S[:] += Si.real
+
         P+=Pi
     S = S/K
+
     P = P/K
     F = fftfreq(N,1.0/Fs)
+
     F = rfftfreq(N,1.0/Fs)
     return (fftshift(S),fftshift(F))
+
     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
NW = 3
+
N = T*Fs
(s,t) = sin(f=10.2,Fs=Fs)
+
NW = 2
s = s+np.random.randn(len(s))
+
N_rep = 1000
(S,F) = mtm(s, NW = NW,  Fs = Fs)
+
S_perio = np.zeros((N_rep, N//2+1))
py.plot(F,S)
+
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( F,fftshift(np.abs(fft(s))**2/len(s) ) ,'g')
+
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')
  
print np.sum(S)
+
py.subplot(3,1,3)
print np.sum(s**2)
+
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()
 
+
-->
</source>
+
==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

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 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

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

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

#!/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.

# -*- 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:

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