Estymacja widma na podstawie FT: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 28 pośrednich wersji utworzonych przez tego samego użytkownika)
Linia 1: Linia 1:
−
=Praktyczna estymacja widma Fourierowskiego sygnałów=
 
  
−
Dla sygnałów dyskretnych obliczamy Dyskretną Transformatę Fouriera (omawianą szerzej na [[Ćwiczenia_2#Dyskretna_Transformata_Fouriera_.28DFT.29|ćwiczeniach]]) oraz poniżej. Kwadrat jej modułu to inaczej periodogram, czyli estymata geśtości widmowej mocy dla sygnałów dyskretnych.
 
  
−
Sygnały z którymi mamy do czynienia w praktyce są nie tylko dyskretne, ale też skończone.
+
=[[Analiza_sygnałów_-_lecture|AS/]] Estymacja widma na podstawie Transformaty Fouriera=
−
Obliczanie transformaty Fouriera dla skończonego odcinka niesie ze
+
 
−
sobą dodatkowe komplikacje.  Znamy wartości sygnału <math>x[n]</math>
+
==Dyskretna Transformata Fouriera (DFT)==
−
dla <math>i=1\ldots N</math>. Odpowiada to iloczynowi sygnału
 
−
<math>\left\{s[n]\right\}_{n\in\mathbb{Z}}</math> z oknem prostokątnym
 
−
<math>w_p[k]</math>:
 
  
−
<math> w_p[k]=\left\{\begin{array}{rl}
+
W praktycznych zastosowaniach mamy do czynienia z sygnałami próbkowanymi o skończonej długości. Transformata Fouriera działająca na takich sygnałach nazywana jest Dyskretną Transformatą Fouriera, a algorytm najczęściej wykorzystywany do jej obliczania to szybka trasnsformata Fouriera (fast Fourier transform FFT).
−
1 & \mathrm{dla} \;k=1 .. N\\
+
Formułę na współczynniki transformaty Fouriera można otrzymać z [[Szereg_Fouriera|szeregu Fouriera]]. Załóżmy, że sygnał który chcemy przetransformować składa się z <math>N</math> próbek.
−
0 & \mathrm{dla} \;k<0 \vee k>N\\
 
−
\end{array}
 
−
\right.
 
−
</math>
 
  
−
W efekcie [[Twierdzenia_o_splocie_i_o_próbkowaniu_(aliasing)|(patrz twierdzenie o splocie)]] otrzymujemy splot transformaty Fouriera
+
<math> s =\{ s[0],\dots,s[n],\dots s[N-1]\}</math>  
−
sygnału (nieskończonego) z transformatą Fouriera okna
 
−
<math>\hat{w}_p[k]</math>. Na przykład dla okna prostokątnego będzie to funkcja postaci <math>sin(x)/x</math>, która może wprowadzić w widmie sztuczne oscylacje, które mylnie możemy zidentyfikować z pikami widma. Dlatego w praktyce stosujemy okna o
 
−
łagodniejszym przebiegu transformaty Fouriera. Czyli:
 
−
#Obliczamy iloczyn sygnału <math>s[n]</math> z wybranym oknem <math>w[n]</math>, dopasowanym do jego rozmiaru
 
−
#Obliczamy periodogram sygnału <math>s[n] w[n]</math>
 
  
−
Przy założeniu stacjonarności sygnału możemy obliczyć widmo [[Nieparametryczne_widmo_mocy#Metoda_Welcha|omawianą na ćwiczeniach]] metodą Welcha, według której dzielimy sygnał na zachodzące na siebie odcinki, każdy odcinek mnożymy przez okno <math>w[n]</math> po czy otrzymane widma uśredniamy. W ten sposób dla każdej częstości mamy po kilka estymat mocy widmowej, wyliczonych z kolejnych odcinków, co pozwala na oszacowanie błędu estymaty.
+
i próbki pobierane były co <math>T_s</math> sekund. Zakładamy, że analizowany sygnał <math>s</math>  to jeden okres nieskończonego sygnału o  okresie  <math>T=N\cdot T_s</math>. Wprowadźmy oznaczenie:
  
 +
<math>s[n]=s(n T_s)</math>.
  
−
===Dyskretna Transformata Fouriera (DFT) ===
+
Przepiszmy wzór na współczynniki [[Szereg_Fouriera|szeregu Fouriera]]
  
−
W praktycznych zastosowaniach mamy do czynienia z sygnałami próbkowanymi o skończonej długości. Transformata Fouriera działąjąca na takich sygnałach nazywana jest Dyskretną Transformatą Fouriera, a algorytm najczęściej wykorzystywany do jej obliczania to szybka trasnsformata Fouriera (fast Fourier transform FFT).
 
−
Formułę na współczynniki FFT można otrzymać z [[Szereg_Fouriera|szeregu Fouriera]]. Załóżmy, że sygnał który chcemy przetransformować składa się z <math>N</math> próbek.
 
  
−
<math> s =\{ s[0],\dots,s[n],\dots s[N-1]\}</math>  
+
<math>
 +
c_{k} = \frac{1}{T}\int_{0}^{T} s(t) e^\frac{2\pi i k t}{T} d t
 +
</math>
  
−
i próbki pobierane były co <math>T_s</math> sekund. Zakładamy, że analizowany sygnał <math>s</math>  to jeden okres nieskończonego sygnału o  okresie  <math>T=N\cdot T_s</math>. Wprowadźmy oznaczenie:
 
  
−
<math>s[n]=s(nT_s)</math>.
+
Ponieważ sygnał jest teraz dyskretny, całka zamieni się na sumę pól prostokątów o bokach równych wartości funkcji podcałkowej w zadanych punktach <math>x(nT_s)e^{(2i{\pi}knT_s/T)}</math> i odległości między punktami <math>T_s</math>:
  
−
Przepiszmy wzór na współczynniki szeregu Fouriera. Ponieważ sygnał jest teraz dyskretny, całka zamieni się na sumę Riemanna: pole będzie sumą pól prostokątów o bokach równych wartości funkcji podcałkowej w zadanych punktach <math>x(nT_s)exp(2i{\pi}knT_s/T)</math> i odległości między punktami <math>T_s</math>:
 
  
 
<math>
 
<math>
−
  S[k] = \frac{1}{NT_s}\sum_{n=0}^{N-1}s(nT_s)e^{2i\pi\frac{knT_s}{NT_s}}T_s = \frac{1}{N}\sum_{n=0}^{N-1}s[n]e^{2i{\pi}\frac{kn}{N}}
+
  \hat{s}[k] = \frac{1}{NT_s}\sum_{n=0}^{N-1}s(nT_s)e^{2i\pi\frac{knT_s}{NT_s}} \; T_s = \frac{1}{N}\sum_{n=0}^{N-1}s[n]e^{2i{\pi}\frac{kn}{N}}
 
</math>
 
</math>
  
  
−
DFT zaimplementowana w <tt>numpy.fft</tt> jest określona jako:
 
−
:<math>A[k] =  \sum_{m=0}^{n-1} a[m] \exp\left\{-2\pi i{mk \over n}\right\}      \qquad k = 0,\ldots,n-1.
 
−
</math>
 
  
−
DFT jest w ogólności zdefiniowane dla zespolonych argumentów i zwraca zespolone współczynniki.
+
==Praktyczna estymacja widma Fourierowskiego sygnałów==
−
Odwrotna dyskretna transformata Fouriera jest zdefiniowana jako:
 
−
:<math> a[m] = \frac{1}{n}\sum_{k=0}^{n-1}A[k]\exp\left\{2\pi i{mk\over n}\right\}
 
−
      \qquad m = 0,\ldots,n-1.
 
−
</math>
 
  
−
Zwróćmy uwagę, że różni się ona do transformaty ''wprost'' jedynie znakiem w exponencie i normalizacją <math>1/n</math>.
+
Dla sygnałów dyskretnych obliczamy Dyskretną Transformatę Fouriera (omawianą też szerzej na [[Ćwiczenia_2#Dyskretna_Transformata_Fouriera_.28DFT.29|ćwiczeniach]]). Kwadrat jej modułu to inaczej periodogram, czyli estymata gęstości widmowej mocy dla sygnałów dyskretnych.  
  
 +
Sygnały z którymi mamy do czynienia w praktyce są nie tylko dyskretne, ale też skończone.
 +
Obliczanie transformaty Fouriera dla skończonego odcinka niesie ze
 +
sobą dodatkowe komplikacje.  Znamy wartości sygnału <math>s[n]</math>
 +
dla <math>i=1\ldots N</math>. Odpowiada to iloczynowi sygnału
 +
<math>\left\{s[n]\right\}_{n\in\mathbb{Z}}</math> z oknem prostokątnym
 +
<math>w_p[k]</math>:
  
−
===Okienkowanie a widmo mocy: periodogram===
+
<math> w_p[k]=\left\{\begin{array}{rl}
−
Przypomnijmy wzór na dyskretną transformatę Fouriera [http://haar.zfb.fuw.edu.pl/edu/index.php/%C4%86wiczenia_2 DFT] zaimplementowaną w FFT:
+
1 & \mathrm{dla} \;k=1 .. N\\
 +
0 & \mathrm{dla} \;k<0 \vee k>N\\
 +
\end{array}
 +
\right.
 +
</math>
  
−
:<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.
+
W efekcie [[Twierdzenia_o_splocie_i_o_próbkowaniu_(aliasing)|(patrz twierdzenie o splocie)]] otrzymujemy splot transformaty Fouriera
−
</math>
+
sygnału (nieskończonego) z transformatą Fouriera okna
 +
<math>\hat{w}_p[k]</math>. Na przykład dla okna prostokątnego będzie to funkcja postaci <math>sin(x)/x</math>, która może wprowadzić w widmie sztuczne oscylacje, które mylnie możemy zidentyfikować z pikami widma. Dlatego w praktyce stosujemy okna o łagodniejszym przebiegu transformaty Fouriera. Czyli:
 +
#Obliczamy iloczyn sygnału <math>s[n]</math> z wybranym oknem <math>w[n]</math>
 +
#Obliczamy periodogram sygnału <math>s[n] w[n]</math>
  
  
−
Na podstawie twierdzenia [[Nieparametryczne_widmo_mocy#Twierdzenie_Parsevala|Parsevala]] możemy policzyć widmo mocy jako:
+
W ogólnym przypadku, biorąc pod uwagę normalizację okna
−
<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.
+
<math>w[n] = \frac{1}{\sqrt{\sum_{n=0}^{N-1} (w[n])^2}}w[n]</math>
  
−
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>
 
 
<!--(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:
 
−
<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>
 
−
 
  
 +
dostajemy widmo mocy sygnału okienkowanego:
  
−
=== 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>
 
<math>
−
P(f) = |X(f)|^2 = X_R^2(f) + X_I^2(f)
+
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>
−
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:
+
Przy założeniu stacjonarności sygnału możemy obliczyć widmo [[Nieparametryczne_widmo_mocy#Metoda_Welcha|testowaną na ćwiczeniach]] metodą Welcha, według której dzielimy sygnał na zachodzące na siebie odcinki, każdy odcinek mnożymy przez okno <math>w[n]</math> po czy otrzymane widma uśredniamy. W ten sposób dla każdej częstości mamy po kilka estymat mocy widmowej, wyliczonych z kolejnych odcinków, co pozwala na oszacowanie błędu estymaty.
  
−
:<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.
+
<div align="right">
 +
[[Twierdzenia_o_splocie_i_o_próbkowaniu_(aliasing)|⬅]] [[Analiza_sygnałów_-_wykład|⬆]] [[Model_autoregresyjny_(AR)|⮕]]
 +
</div>

Aktualna wersja na dzień 10:14, 8 lis 2024


AS/ Estymacja widma na podstawie Transformaty Fouriera

Dyskretna Transformata Fouriera (DFT)

W praktycznych zastosowaniach mamy do czynienia z sygnałami próbkowanymi o skończonej długości. Transformata Fouriera działająca na takich sygnałach nazywana jest Dyskretną Transformatą Fouriera, a algorytm najczęściej wykorzystywany do jej obliczania to szybka trasnsformata Fouriera (fast Fourier transform FFT). Formułę na współczynniki transformaty Fouriera można otrzymać z szeregu Fouriera. Załóżmy, że sygnał który chcemy przetransformować składa się z [math]N[/math] próbek.

[math] s =\{ s[0],\dots,s[n],\dots s[N-1]\}[/math]

i próbki pobierane były co [math]T_s[/math] sekund. Zakładamy, że analizowany sygnał [math]s[/math] to jeden okres nieskończonego sygnału o okresie [math]T=N\cdot T_s[/math]. Wprowadźmy oznaczenie:

[math]s[n]=s(n T_s)[/math].

Przepiszmy wzór na współczynniki szeregu Fouriera


[math] c_{k} = \frac{1}{T}\int_{0}^{T} s(t) e^\frac{2\pi i k t}{T} d t [/math]


Ponieważ sygnał jest teraz dyskretny, całka zamieni się na sumę pól prostokątów o bokach równych wartości funkcji podcałkowej w zadanych punktach [math]x(nT_s)e^{(2i{\pi}knT_s/T)}[/math] i odległości między punktami [math]T_s[/math]:


[math] \hat{s}[k] = \frac{1}{NT_s}\sum_{n=0}^{N-1}s(nT_s)e^{2i\pi\frac{knT_s}{NT_s}} \; T_s = \frac{1}{N}\sum_{n=0}^{N-1}s[n]e^{2i{\pi}\frac{kn}{N}} [/math]


Praktyczna estymacja widma Fourierowskiego sygnałów

Dla sygnałów dyskretnych obliczamy Dyskretną Transformatę Fouriera (omawianą też szerzej na ćwiczeniach). Kwadrat jej modułu to inaczej periodogram, czyli estymata gęstości widmowej mocy dla sygnałów dyskretnych.

Sygnały z którymi mamy do czynienia w praktyce są nie tylko dyskretne, ale też skończone. Obliczanie transformaty Fouriera dla skończonego odcinka niesie ze sobą dodatkowe komplikacje. Znamy wartości sygnału [math]s[n][/math] dla [math]i=1\ldots N[/math]. Odpowiada to iloczynowi sygnału [math]\left\{s[n]\right\}_{n\in\mathbb{Z}}[/math] z oknem prostokątnym [math]w_p[k][/math]:

[math] w_p[k]=\left\{\begin{array}{rl} 1 & \mathrm{dla} \;k=1 .. N\\ 0 & \mathrm{dla} \;k\lt 0 \vee k\gt N\\ \end{array} \right. [/math]

W efekcie (patrz twierdzenie o splocie) otrzymujemy splot transformaty Fouriera sygnału (nieskończonego) z transformatą Fouriera okna [math]\hat{w}_p[k][/math]. Na przykład dla okna prostokątnego będzie to funkcja postaci [math]sin(x)/x[/math], która może wprowadzić w widmie sztuczne oscylacje, które mylnie możemy zidentyfikować z pikami widma. Dlatego w praktyce stosujemy okna o łagodniejszym przebiegu transformaty Fouriera. Czyli:

  1. Obliczamy iloczyn sygnału [math]s[n][/math] z wybranym oknem [math]w[n][/math]
  2. Obliczamy periodogram sygnału [math]s[n] w[n][/math]


W ogólnym przypadku, biorąc pod uwagę normalizację okna

[math]w[n] = \frac{1}{\sqrt{\sum_{n=0}^{N-1} (w[n])^2}}w[n][/math]


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

Przy założeniu stacjonarności sygnału możemy obliczyć widmo testowaną na ćwiczeniach metodą Welcha, według której dzielimy sygnał na zachodzące na siebie odcinki, każdy odcinek mnożymy przez okno [math]w[n][/math] po czy otrzymane widma uśredniamy. W ten sposób dla każdej częstości mamy po kilka estymat mocy widmowej, wyliczonych z kolejnych odcinków, co pozwala na oszacowanie błędu estymaty.


Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=Estymacja_widma_na_podstawie_FT&oldid=11133”