Nieparametryczne widmo mocy

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ą następujący algorytm estymacji widma mocy:

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

Następnie proszę obliczyć energię oraz wyświetlić przebieg widma mocy dla sygnału z Zadania 1.

#!/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 = np.ones(len(s))
    okno= okno/np.linalg.norm(okno)
    s = s*okno
    S = rfft(s)
    P = S*S.conj()/Fs
    P = P.real
    if len(s)%2 ==0:
        P[1:-1] *=2
    else:
        P[1:] *=2
    F = rfftfreq(len(s), 1./Fs)
    return P,F

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)

fs =100
T = 3
(s,t) = sin(f=10.3,T= T,Fs=fs)
moc_w_czasie = s**2
(moc_w_czestosci, F) = widmo_mocy(s, Fs=fs)
dt = 1/fs

energia_w_czasie = np.sum(moc_w_czasie*dt)
energia_w_czestosci = np.sum(moc_w_czestosci)

py.subplot(3,1,1)
py.plot(t,s)
py.title('Sygnal')
py.subplot(3,1,2)
py.plot(t,moc_w_czasie)
py.title('moc w czasie, energia w czasie: ' +str(energia_w_czasie))
py.subplot(3,1,3)
py.plot(F,moc_w_czestosci)
py.title('moc w czestosci, energia w czestosci: ' +str(energia_w_czestosci))
py.show()

Okienkowanie a widmo mocy: periodogram

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

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

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

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

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

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

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

Zadanie 3: Obliczanie periodogramu

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

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

Sygnały stochastyczne

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

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

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

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

Dla przypomnienia wartość oczekiwaną liczymy tak:

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

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

Stacjonarność i ergodyczność

Stacjonarność:

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

Ergodyczność:

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

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

Zadanie 4: Transformata Fouriera sygnału stochastycznego

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

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

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

Proszę obejrzeć otrzymane widma.

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

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

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

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

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

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

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

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

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

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

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

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

Zadanie 5: Metoda Welcha

Proszę zapoznać się zaimplementowaną w bibliotece scipy.signal funkcją welch. Funkcję proszę przetestować obliczając dla funkcji sinus energię sygnału w dziedzinie czasu i w dziedzinie częstości. Testy proszę wykonać dla okna prostokątnego, Blackmana i Haminga.

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

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.

Zadanie 7: Estymacja widma mocy metodą multitaper

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

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

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

*

Analiza_sygnałów_-_ćwiczenia/Fourier_4