Ćwiczenia 2: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 30 wersji utworzonych przez 2 użytkowników)
Linia 1: Linia 1:
 +
[[Analiza_sygnałów_-_ćwiczenia]]/Fourier
 +
 
==Transformata Fouriera==
 
==Transformata Fouriera==
=== Szereg Fouriera ===
+
===Analiza Fouriera jako rzutowanie wektorów ===
Ciągłą funkcję  okresową (sygnał okresowy) można przedstawić w postaci  [[Szereg_Fouriera|szeregu Fouriera]] &mdash; czyli sumy zespolonych funkcji exponencjalnych. Współczynniki w szeregu Fouriera to liczby, w ogólności zespolone, które mówią z jaką amplitudą i z jaką fazą wchodzi odpowiadająca mu  funkcja, współczynniki te numerowane są od <math>-\infty</math> do <math>\infty</math>. Dla sygnałów rzeczywistych współczynniki występują parami co powoduje, że zgodnie ze wzorami Eulera:
+
W tym zadaniu chciałbym abyśmy wyrobili sobie intuicję, że analiza Fouriera moze być rozumiana jako badanie sygnału (niech nazywa się syg) za pomocą pewnych znanych funkcji (będą to funkcje sin i cos). Jak już mówiliśmy o sygnałach można myśleć, że są to wektory. Jeśli chcemy zbadać, czy dwa wektory są podobne do siebie to liczymy ich iloczyn skalarny. Zobaczmy jak to działą w poniższym zadaniu:
 +
 
 +
<source lang = python>
 +
import numpy as np
 +
 
 +
T = 1 # ustawiamy długość czasu
 +
Fs = 100 # ustawiamy częstość próbkowania
 +
t = #wytwarzamy wektor czasu
 +
 
 +
syg = # wytwarzamy sygnał będący sumą sinusa o cz. 3Hz i cos o częstości 7Hz
 +
syg = syg/np.linalg.norm(syg) # normalizujemy badany sygnał
 +
 
 +
print("częstość\tsin\tcos" ) # szykujemy nagłówek wydruku
 +
for n in # dla liczb n z zakresu 1 do 10
 +
    test_sin = #tworzymy sygnał testowy sinus o częstości n
 +
    test_sin = #normalizujemy ten sygnał testowy
 +
    test_cos = # tworzymy testowego cosinusa o częstości n
 +
    test_cos = # normalizujemy
 +
   
 +
    rzut_s =#obliczamy iloczyn skalarny sygnału z testowym sinusem
 +
    rzut_p = #obliczamy iloczyn skalarny z testowym cosinusem
 +
    print("%d\t%.3f\t%.3f" %(n,rzut_s, rzut_p)) # wypisujemy wyniki
 +
</source>
 +
 
 +
 
 +
=== Zespolone eksponensy ===
 +
<!--Ciągłą funkcję  okresową (sygnał okresowy) można przedstawić w postaci  [[Szereg_Fouriera|szeregu Fouriera]] &mdash; czyli sumy zespolonych funkcji exponencjalnych. Współczynniki w szeregu Fouriera to liczby, w ogólności zespolone, które mówią z jaką amplitudą i z jaką fazą wchodzi odpowiadająca mu  funkcja, współczynniki te numerowane są od <math>-\infty</math> do <math>\infty</math>. Dla sygnałów rzeczywistych współczynniki występują parami co powoduje, że zgodnie ze wzorami -->
 +
Dla przypomnienia wzory Eulera:
  
 
<math>
 
<math>
Linia 7: Linia 36:
 
</math>
 
</math>
  
szereg można zwinąć do szeregu funkcji sinusoidalnych.
+
<!--szereg można zwinąć do szeregu funkcji sinusoidalnych.-->
  
 
=== Zadanie ===
 
=== Zadanie ===
 +
Dla przypomnienia:
 +
* uzyskanie liczby zespolonej <math>i</math> w pythonie mamy jako <tt> 1j</tt>.
 +
* jeśli prowadzimy obliczenia na liczbach zespolonych to w wyniku dostajemy liczby zespolone (nawet jeśli ich część urojona jest 0). Aby rzutować wynik na liczby rzeczywiste musimy wykonać funkcję <tt>np.real(.)</tt>
 
Proszę narysować wykres funkcji <math>f(t) = c \left(e^{i2 \pi \frac{  n  }{T}t } + e^{-i 2 \pi \frac{ n }{T} t  }\right) </math> dla <math>t \in (-\pi,\, \pi); \quad T = 2 \pi</math>.  
 
Proszę narysować wykres funkcji <math>f(t) = c \left(e^{i2 \pi \frac{  n  }{T}t } + e^{-i 2 \pi \frac{ n }{T} t  }\right) </math> dla <math>t \in (-\pi,\, \pi); \quad T = 2 \pi</math>.  
 +
Odstęp między próbkami proszę ustawić na 0.001. 
 +
 
* Obejrzeć wykresy dla <math>n \in \{0, 1, 2, 5\}</math>.  
 
* Obejrzeć wykresy dla <math>n \in \{0, 1, 2, 5\}</math>.  
 
* Jaki sens ma liczba <math>n</math>?  
 
* Jaki sens ma liczba <math>n</math>?  
Linia 16: Linia 50:
  
 
===  Obliczanie transformaty ===
 
===  Obliczanie transformaty ===
====Przykład ====
 
'''Oblicz analitycznie''',  [[Szereg_Fouriera|korzystając ze wzorów]] współczynniki Fouriera dla funkcji ''f''. Za pomocą wykresu coraz to większej ilości składników otrzymanego szeregu sprawdź zbieżność szeregu.
 
 
Funkcja  o okresie 2, w podstawowym okresie dana jest ona wzorem:
 
Funkcja  o okresie 2, w podstawowym okresie dana jest ona wzorem:
 
:<math> f(t)=|t| </math> dla <math>|t|<1  </math>,
 
:<math> f(t)=|t| </math> dla <math>|t|<1  </math>,
 +
Za pomocą wykresu coraz to większej ilości składników jej szeregu Fouriera sprawdź jego  zbieżność.
 +
===Obliczenia analityczne===
 +
<!--
 +
Za pomocą wykresu coraz to większej ilości składników otrzymanego szeregu sprawdź zbieżność szeregu.
 +
[[Szereg_Fouriera|korzystając ze wzorów]] współczynniki Fouriera dla funkcji ''f''.
 +
-->
 +
 +
  
 
Musimy obliczyć następującą całkę:
 
Musimy obliczyć następującą całkę:
 
:<math> c_n = \frac{1}{2} \int_{-1}^1 {|t| e^{\frac{2 \pi j n t }{2}} dt} =  
 
:<math> c_n = \frac{1}{2} \int_{-1}^1 {|t| e^{\frac{2 \pi j n t }{2}} dt} =  
 
\frac{1}{2} \left[ \int_{-1}^0{ -t e^{\pi j n t }dt} + \int_{0}^1 {t e^{\pi j n t }dt} \right]</math>
 
\frac{1}{2} \left[ \int_{-1}^0{ -t e^{\pi j n t }dt} + \int_{0}^1 {t e^{\pi j n t }dt} \right]</math>
Przypomnijmy sobie, że [http://www.wolframalpha.com/index.html np. korzystając ze strony WolframAlpha]:
+
Przypomnijmy sobie, że [http://www.wolframalpha.com/input/?i=integrate+x+exp(ax)+dx np. korzystając ze strony WolframAlpha]:
 
:<math>\int x e^{ax}dx = \frac{1}{a^2} e^{ax} (ax-1) </math>
 
:<math>\int x e^{ax}dx = \frac{1}{a^2} e^{ax} (ax-1) </math>
 
Zatem nasze <math>c_n</math>:
 
Zatem nasze <math>c_n</math>:
Linia 39: Linia 79:
 
Korzystając ze wzorów Eulera możemy zwinąć to wyrażenie do funkcji trygonometrycznych:
 
Korzystając ze wzorów Eulera możemy zwinąć to wyrażenie do funkcji trygonometrycznych:
 
:<math>c_n =  \frac{1}{ \pi^2 n^2} (\cos(\pi n ) - \sin(\pi n)  -1)</math>
 
:<math>c_n =  \frac{1}{ \pi^2 n^2} (\cos(\pi n ) - \sin(\pi n)  -1)</math>
'''Rekonstrukcja sygnałów: '''
+
 
 +
===Rekonstrukcja sygnałów: ===
 
Na podstawie wyników poprzedniego zadania proszę napisać program, który demonstruje jaki jest wynik składania coraz większej ilości czynników.
 
Na podstawie wyników poprzedniego zadania proszę napisać program, który demonstruje jaki jest wynik składania coraz większej ilości czynników.
 
Poniższy kod implementuje sumowanie zadanej ilości składników szeregu i ilustruje wynik.
 
Poniższy kod implementuje sumowanie zadanej ilości składników szeregu i ilustruje wynik.
 
Proszę uruchomić go dla <tt>N={2, 4, 6, 20}</tt>
 
Proszę uruchomić go dla <tt>N={2, 4, 6, 20}</tt>
  
W implementacji musimy zwrócić uwagę na fakt, że <math>\lim_{n \rightarrow 0} \frac{1}{ \pi^2 n^2} \left(\cos(\pi n )-\sin(\pi n)  -1\right) = \frac{1}{2}  </math>
+
W implementacji musimy zwrócić uwagę na fakt, że:
zaś dla pozostałych ''n'' całkowitych możemy uwzględnić, że <math> \sin(\pi n) = 0</math>.
+
* <math>\lim_{n \rightarrow 0} \frac{1}{ \pi^2 n^2} \left(\cos(\pi n )-\sin(\pi n)  -1\right) = \frac{1}{2}  </math>
 +
* zaś dla pozostałych ''n'' całkowitych możemy uwzględnić, że <math> \sin(\pi n) = 0</math>.
 +
* proszę zwrócić uwagę, że otrzymywane sumy są równe swojej części rzeczywistej
 +
* jak wygląda rekonstruowany sygnał poza przedziałem t  =(-1,1) ?
 +
 
 
<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
 
+
 
def skladnik_modul_t(n ):
 
def skladnik_modul_t(n ):
'''n - który element szeregu'''
+
'''Funkcja obliczająca n-ty element szeregu Fouriera dla funkcji moduł t
 +
          n - który element szeregu'''
 +
   
 
if n==0:
 
if n==0:
c_n = 0.5
+
c_n = ...
 
else:
 
else:
c_n = 1/(np.pi**2 * n**2) * (-1  + np.cos(np.pi*n) )
+
c_n =...
 
return c_n
 
return c_n
+
t = np.arange(-4,4,0.1)
+
t = np.arange(-1,1,0.1)
f = np.zeros(len(t))
+
syg = np.zeros(len(t),dtype='complex')
 
+
N = 10
+
N = 1 # ustalamy ile par zespolonych eksponensów sumujemy
c = []
 
 
for n in range(-N,N+1):
 
for n in range(-N,N+1):
cn = skladnik_modul_t(n)
+
c = skladnik_modul_t(n)
c.append(cn)
+
print ('n= %d, c= %.2f'%(n,c))
f += cn * np.exp(-1j*2*np.pi*t*n/2.0)
+
syg += c * ...
py.plot(t,f)
+
py.plot(t,syg,'b', t, syg.imag,'r')
for n in range(-N,N+1):
+
py.ylim((-0.1, 1.1))
print '%(c).4f'%{'c':c[n+N]}
 
py.show()
 
 
</source>
 
</source>
  
====Zadanie ====
+
 
Analogicznie do powyższego przykładu proszę znaleźć i zbadać szereg Fouriera dla funkcji:
+
====Zadanie domowe====
 +
* Analogicznie do powyższego przykładu proszę znaleźć i zbadać szereg Fouriera dla funkcji:
 
: <math> f(t) = |cos(t)|</math>
 
: <math> f(t) = |cos(t)|</math>
 +
* zaimplementować ilustrację sumowania szeregu Fouriera dla sygnału prostokątnego (przykład z wykładu): [[Szereg_Fouriera]]
  
 
===Dyskretna Transformata Fouriera (DFT) ===
 
===Dyskretna Transformata Fouriera (DFT) ===
Linia 108: Linia 153:
  
 
Zwróćmy uwagę, że różni się ona do transformaty ''wprost'' jedynie znakiem w exponencie i normalizacją <math>1/n</math>.
 
Zwróćmy uwagę, że różni się ona do transformaty ''wprost'' jedynie znakiem w exponencie i normalizacją <math>1/n</math>.
   
+
 +
===FFT===
 +
* Dokumentacja numpy: http://docs.scipy.org/doc/numpy/reference/routines.fft.html
 +
* numpy implementuje algorytm: [https://www.jstor.org/stable/2003354?seq=1#page_scan_tab_contents  "An Algorithm for the Machine Calculation of Complex Fourier Series"] James W. Cooley and John W. Tukey Mathematics of Computation Vol. 19, No. 90 (Apr., 1965), pp. 297-301.
 +
 
 
Wartości zwracane przez <tt>fft(a,n)</tt> (<tt>a</tt> sygnał, <tt>n</tt> ilość punktów transformaty) mają następujący ''standardowy'' porządek:     
 
Wartości zwracane przez <tt>fft(a,n)</tt> (<tt>a</tt> sygnał, <tt>n</tt> ilość punktów transformaty) mają następujący ''standardowy'' porządek:     
 
Jeśli <tt>A = fft(a, n)</tt>, to  
 
Jeśli <tt>A = fft(a, n)</tt>, to  
Linia 131: Linia 180:
 
Jeśli potraktujemy wejście <tt>a</tt> jako sygnał w dziedzinie czasu i policzymy <tt>A = fft(a)</tt>, wówczas <tt>np.abs(A)</tt> jest widmem amplitudowym, zaś <tt>np.abs(A)**2</tt> jest widmem mocy. Można obliczyć także widmo fazowe za pomocą funkcji  <tt>np.angle(A)</tt>.
 
Jeśli potraktujemy wejście <tt>a</tt> jako sygnał w dziedzinie czasu i policzymy <tt>A = fft(a)</tt>, wówczas <tt>np.abs(A)</tt> jest widmem amplitudowym, zaś <tt>np.abs(A)**2</tt> jest widmem mocy. Można obliczyć także widmo fazowe za pomocą funkcji  <tt>np.angle(A)</tt>.
  
=== Transformaty rzeczywiste i Hermitowskie ===
+
==== Zadanie ====
Jeśli sygnał wejściowy jest rzeczywisty to jego transformata jest hermitowska, tzn. współczynnik przy częstości <math>f_k</math> jest sprzężony ze współczynnikiem przy częstości <math>-f_k</math>. Oznacza to, że dla sygnałów rzeczywistych współczynniki przy ujemnych częstościach nie wnoszą żadnej dodatkowej informacji. Rodzina funkcji  <tt>rfft</tt> wykorzystuje tą symetrię i zwracają tylko dodatnią część widma włącznie z częstością Nyquista. Tak więc, <tt>n</tt> punktów rzeczywistych na wejściu daje na wyjściu <math>\frac{n}{2}+1</math> punktów zespolonych. Funkcje odwrotne w tej rodzinie zakładają tą samą symetrię i aby na wyjściu uzyskać <tt>n</tt>  punktów rzeczywistych na wejściu trzeba podać <math>\frac{n}{2}+1</math> wartości zespolonych.
+
 
 +
* Proszę wygenerować cosinusoidę oraz sinusoidę o długości 1 s, częstości 2 Hz próbkowaną 10 Hz. Proszę obliczyć współczynniki i wypisać je.
 +
** Jakiego typu liczby otrzymaliśmy?
 +
** Czy istnieją jakieś związki między współczynnikami?
 +
** Jaki jest związek między długością wejściowego sygnału i wyjściowych współczynników?
 +
 
 +
*Proszę wygenerować sygnał będący sumą sinusoid o częstościach f = 30 Hz i f = 21 Hz, amplitudach A = 2 i A = 5, fazach pi/3 i pi/4 oraz składowej stałej 0.5, o długości T = 5 s próbkowany 128 Hz.
 +
**Narysuj wygenerowany sygnał.
 +
**Oblicz współczynniki i wypisz je. Jakiego rodzaju liczby otrzymaliśmy? Jaki jest związek między długością wejściowego sygnału i wyjściowych współczynników?
 +
** Jaki jest związek pomiędzy argumentami zwróconych współczynników a fazą rzeczywistych sygnałów wejściowych?
 +
 
 +
=== Widmo mocy dla sygnałów rzeczywistych: Transformaty rzeczywiste i Hermitowskie ===
 +
Jeśli sygnał wejściowy jest rzeczywisty to jego transformata jest hermitowska, tzn. współczynnik przy częstości <math>f_k</math> jest sprzężony ze współczynnikiem przy częstości <math>-f_k</math>. Oznacza to, że dla sygnałów rzeczywistych współczynniki przy ujemnych częstościach nie wnoszą żadnej dodatkowej informacji.  
 +
 
 +
'''Rodzina funkcji''' <tt>rfft</tt> wykorzystuje tą symetrię i zwracają tylko dodatnią część widma włącznie z częstością Nyquista. Tak więc, <tt>n</tt> punktów rzeczywistych na wejściu daje na wyjściu <math>\frac{n}{2}+1</math> punktów zespolonych. Funkcje odwrotne w tej rodzinie zakładają tą samą symetrię i aby na wyjściu uzyskać <tt>n</tt>  punktów rzeczywistych na wejściu trzeba podać <math>\frac{n}{2}+1</math> wartości zespolonych.
 
      
 
      
 
Dla kompletności powiedzmy jeszcze, że możliwy jest przypadek odwrotny, tzn. widmo jest czysto rzeczywiste i odpowiada mu hermitowski sygnał zespolony. Tę symetrię wykorzystują funkcje <tt>hfft</tt>, które zakładają, że operujemy w dziedzinie czasu <math>\frac{n}{2}+1</math> punktami zespolonymi i odpowiadającymi im w dziedzinie częstości <math>n</math> punktami rzeczywistymi.
 
Dla kompletności powiedzmy jeszcze, że możliwy jest przypadek odwrotny, tzn. widmo jest czysto rzeczywiste i odpowiada mu hermitowski sygnał zespolony. Tę symetrię wykorzystują funkcje <tt>hfft</tt>, które zakładają, że operujemy w dziedzinie czasu <math>\frac{n}{2}+1</math> punktami zespolonymi i odpowiadającymi im w dziedzinie częstości <math>n</math> punktami rzeczywistymi.
 +
===Zadanie===
 +
* Pobierz sygnał http://www.fuw.edu.pl/~jarekz/SYGNALY/klawisz.b.
 +
* wczytaj ten sygnał jako jeden kanał typu int16
 +
* wykreśl ten sygnał
 +
* Jest to dźwięk próbkowany częstościa Fs = 44100
 +
* można go posłuchać za omocą kodu:
 +
 +
<source lang = python>
 +
import sounddevice as sd
 +
sd.play(syg, Fs)
 +
</source>
 +
 +
* wykreśl '''widmo mocy''' tego sygnału tj. '''kwadrat modułu transformaty Fouriera''' (jest to syg. rzeczywisty więc używamy funkcji z rodziny rfft)
 +
 +
 +
;przy okazji: wczytywanie pliku wav, tak aby był on tablicą numpy można zrobić jak w poniższym kodzie, wówczas w Fs mamy częstośc próbkowania sygnału audio, a w macierzy syg mamy  sygnał (może być stereo: lewy i prawy kanał)
 +
<source lang = python>
 +
from scipy.io.wavfile import read
 +
import sounddevice as sd
 +
(Fs, syg) = read('plik.wav')
 +
L = syg[:,0]
 +
sd.play(syg[:,0], Fs)
 +
</source>
 +
 +
==Co musimy z tego zapamiętać?==
 +
* Analiza Fourierowska polega na rzutowaniu sygnału na sygnały próbne: sinusoidy o znanych częstościach
 +
* Ze współczynników transformaty można odtworzyć oryginalny sygnał: mamy więc dwie reprezentacje sygnału:
 +
** w dziedzinie czasu
 +
** w dziedzinie częstości
 +
** między reprezentacjami przechodzimy za pomocą transformaty Fouriera i odwrotnej transformaty Fouriera
 +
* reprezentacja częstościowa jest na ogół zespolona: widmo mocy to kwadrat modułu transformaty Fouriera
 +
[[Analiza_sygnałów_-_ćwiczenia]]/Fourier

Aktualna wersja na dzień 09:14, 2 paź 2017

Analiza_sygnałów_-_ćwiczenia/Fourier

Transformata Fouriera

Analiza Fouriera jako rzutowanie wektorów

W tym zadaniu chciałbym abyśmy wyrobili sobie intuicję, że analiza Fouriera moze być rozumiana jako badanie sygnału (niech nazywa się syg) za pomocą pewnych znanych funkcji (będą to funkcje sin i cos). Jak już mówiliśmy o sygnałach można myśleć, że są to wektory. Jeśli chcemy zbadać, czy dwa wektory są podobne do siebie to liczymy ich iloczyn skalarny. Zobaczmy jak to działą w poniższym zadaniu:

import numpy as np

T = 1 # ustawiamy długość czasu
Fs = 100 # ustawiamy częstość próbkowania
t = #wytwarzamy wektor czasu

syg = # wytwarzamy sygnał będący sumą sinusa o cz. 3Hz i cos o częstości 7Hz
syg = syg/np.linalg.norm(syg) # normalizujemy badany sygnał

print("częstość\tsin\tcos" ) # szykujemy nagłówek wydruku
for n in # dla liczb n z zakresu 1 do 10
    test_sin = #tworzymy sygnał testowy sinus o częstości n
    test_sin = #normalizujemy ten sygnał testowy
    test_cos = # tworzymy testowego cosinusa o częstości n
    test_cos = # normalizujemy 
    
    rzut_s =#obliczamy iloczyn skalarny sygnału z testowym sinusem
    rzut_p = #obliczamy iloczyn skalarny z testowym cosinusem
    print("%d\t%.3f\t%.3f" %(n,rzut_s, rzut_p)) # wypisujemy wyniki


Zespolone eksponensy

Dla przypomnienia wzory Eulera:

[math] e^{ix}=\cos x + i\sin x \quad \Rightarrow \begin{cases}\cos x = \frac12(e^{ix}+e^{-ix})\\ \sin x = \frac{1}{2i}(e^{ix}-e^{-ix}) \end{cases} [/math]


Zadanie

Dla przypomnienia:

  • uzyskanie liczby zespolonej [math]i[/math] w pythonie mamy jako 1j.
  • jeśli prowadzimy obliczenia na liczbach zespolonych to w wyniku dostajemy liczby zespolone (nawet jeśli ich część urojona jest 0). Aby rzutować wynik na liczby rzeczywiste musimy wykonać funkcję np.real(.)

Proszę narysować wykres funkcji [math]f(t) = c \left(e^{i2 \pi \frac{ n }{T}t } + e^{-i 2 \pi \frac{ n }{T} t }\right) [/math] dla [math]t \in (-\pi,\, \pi); \quad T = 2 \pi[/math]. Odstęp między próbkami proszę ustawić na 0.001.

  • Obejrzeć wykresy dla [math]n \in \{0, 1, 2, 5\}[/math].
  • Jaki sens ma liczba [math]n[/math]?
  • Jaki sens ma liczba [math]c[/math]?

Obliczanie transformaty

Funkcja o okresie 2, w podstawowym okresie dana jest ona wzorem:

[math] f(t)=|t| [/math] dla [math]|t|\lt 1 [/math],

Za pomocą wykresu coraz to większej ilości składników jej szeregu Fouriera sprawdź jego zbieżność.

Obliczenia analityczne

Musimy obliczyć następującą całkę:

[math] c_n = \frac{1}{2} \int_{-1}^1 {|t| e^{\frac{2 \pi j n t }{2}} dt} = \frac{1}{2} \left[ \int_{-1}^0{ -t e^{\pi j n t }dt} + \int_{0}^1 {t e^{\pi j n t }dt} \right][/math]

Przypomnijmy sobie, że np. korzystając ze strony WolframAlpha:

[math]\int x e^{ax}dx = \frac{1}{a^2} e^{ax} (ax-1) [/math]

Zatem nasze [math]c_n[/math]:

[math]c_n = \frac{1}{\pi^2 n^2} \left[ - \left[ e^{\pi j n t} \frac{\pi j n t -1}{(\pi j n)^2} \right]_{-1}^0 + \left[ e^{\pi j n t} \frac{\pi j n t -1 }{(\pi j n)^2} \right]_{0}^1 \right] = [/math]
[math] = \frac{1}{2 \pi^2 n^2} \left[ -2+ \pi j n \left( e^{-\pi j n} - e^{\pi j n} \right) + \left( e^{-\pi j n} + e^{\pi j n} \right) \right] [/math]

Korzystając ze wzorów Eulera możemy zwinąć to wyrażenie do funkcji trygonometrycznych:

[math]c_n = \frac{1}{ \pi^2 n^2} (\cos(\pi n ) - \sin(\pi n) -1)[/math]

Rekonstrukcja sygnałów:

Na podstawie wyników poprzedniego zadania proszę napisać program, który demonstruje jaki jest wynik składania coraz większej ilości czynników. Poniższy kod implementuje sumowanie zadanej ilości składników szeregu i ilustruje wynik. Proszę uruchomić go dla N={2, 4, 6, 20}

W implementacji musimy zwrócić uwagę na fakt, że:

  • [math]\lim_{n \rightarrow 0} \frac{1}{ \pi^2 n^2} \left(\cos(\pi n )-\sin(\pi n) -1\right) = \frac{1}{2} [/math]
  • zaś dla pozostałych n całkowitych możemy uwzględnić, że [math] \sin(\pi n) = 0[/math].
  • proszę zwrócić uwagę, że otrzymywane sumy są równe swojej części rzeczywistej
  • jak wygląda rekonstruowany sygnał poza przedziałem t =(-1,1) ?
# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
 
def skladnik_modul_t(n ):
	'''Funkcja obliczająca n-ty element szeregu Fouriera dla funkcji moduł t
           n - który element szeregu'''
     
	if n==0:
		c_n = ...
	else:
		c_n =...
	return c_n
 
t = np.arange(-1,1,0.1)
syg = np.zeros(len(t),dtype='complex')
 
N = 1 # ustalamy ile par zespolonych eksponensów sumujemy
for n in range(-N,N+1):
	c = skladnik_modul_t(n)
	print ('n= %d, c= %.2f'%(n,c))
	syg += c * ...
py.plot(t,syg,'b', t, syg.imag,'r')
py.ylim((-0.1, 1.1))


Zadanie domowe

  • Analogicznie do powyższego przykładu proszę znaleźć i zbadać szereg Fouriera dla funkcji:
[math] f(t) = |cos(t)|[/math]
  • zaimplementować ilustrację sumowania szeregu Fouriera dla sygnału prostokątnego (przykład z wykładu): Szereg_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łą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 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(nT_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] 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]


DFT zaimplementowana w numpy.fft 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. 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].

FFT

Wartości zwracane przez fft(a,n) (a sygnał, n ilość punktów transformaty) mają następujący standardowy porządek: Jeśli A = fft(a, n), to

  • A[0] zawiera składową stałą (średnią sygnału)
  • A[1:n/2] zawiera współczynniki odpowiadające dodatnim częstościom
  • A[n/2+1:] zawiera współczynniki odpowiadające ujemnym częstościom w kolejności od bardziej do mniej ujemnych.
  • Dla parzystego n A[n/2] reprezentuje dodatnia i ujemną częstość Nyquista i dla sygnałów rzeczywistych jest liczbą rzeczywistą.
  • Dla nieparzystego n, element A[(n-1)/2] zawiera współczynnik dla największej częstości dodatniej a element A[(n+1)/2] zawiera współczynnik dla największej częstości ujemnej.

Funkcja numpy.fft.fftfreq(len(A),1.0/Fs) zwraca macierz częstości odpowiadających poszczególnym elementom wyjściowym.

Składnia: numpy.fft.fftfreq(n, d=1.0) Parametry:

n : int — długość okna.
d : skalar — okres próbkowania (odwrotność częstości próbkowania).
Zwracane częstości są obliczane w następujący sposób:
f = [0,1,...,n/2-1,-n/2,...,-1]/(d*n) jeśli n jest przyste
f = [0,1,...,(n-1)/2,-(n-1)/2,...,-1]/(d*n) jeśli n jest nieparzyste.

Funkcja numpy.fft.fftshift(A) przestawia wektor wyjściowy fft i wektor częstości, tak aby częstość zero wypadała w środku. Zastosowanie funkcji numpy.fft.ifftshift(A) odwraca działanie numpy.fft.fftshift(.).

Jeśli potraktujemy wejście a jako sygnał w dziedzinie czasu i policzymy A = fft(a), wówczas np.abs(A) jest widmem amplitudowym, zaś np.abs(A)**2 jest widmem mocy. Można obliczyć także widmo fazowe za pomocą funkcji np.angle(A).

Zadanie

  • Proszę wygenerować cosinusoidę oraz sinusoidę o długości 1 s, częstości 2 Hz próbkowaną 10 Hz. Proszę obliczyć współczynniki i wypisać je.
    • Jakiego typu liczby otrzymaliśmy?
    • Czy istnieją jakieś związki między współczynnikami?
    • Jaki jest związek między długością wejściowego sygnału i wyjściowych współczynników?
  • Proszę wygenerować sygnał będący sumą sinusoid o częstościach f = 30 Hz i f = 21 Hz, amplitudach A = 2 i A = 5, fazach pi/3 i pi/4 oraz składowej stałej 0.5, o długości T = 5 s próbkowany 128 Hz.
    • Narysuj wygenerowany sygnał.
    • Oblicz współczynniki i wypisz je. Jakiego rodzaju liczby otrzymaliśmy? Jaki jest związek między długością wejściowego sygnału i wyjściowych współczynników?
    • Jaki jest związek pomiędzy argumentami zwróconych współczynników a fazą rzeczywistych sygnałów wejściowych?

Widmo mocy dla sygnałów rzeczywistych: Transformaty rzeczywiste i Hermitowskie

Jeśli sygnał wejściowy jest rzeczywisty to jego transformata jest hermitowska, tzn. współczynnik przy częstości [math]f_k[/math] jest sprzężony ze współczynnikiem przy częstości [math]-f_k[/math]. Oznacza to, że dla sygnałów rzeczywistych współczynniki przy ujemnych częstościach nie wnoszą żadnej dodatkowej informacji.

Rodzina funkcji rfft wykorzystuje tą symetrię i zwracają tylko dodatnią część widma włącznie z częstością Nyquista. Tak więc, n punktów rzeczywistych na wejściu daje na wyjściu [math]\frac{n}{2}+1[/math] punktów zespolonych. Funkcje odwrotne w tej rodzinie zakładają tą samą symetrię i aby na wyjściu uzyskać n punktów rzeczywistych na wejściu trzeba podać [math]\frac{n}{2}+1[/math] wartości zespolonych.

Dla kompletności powiedzmy jeszcze, że możliwy jest przypadek odwrotny, tzn. widmo jest czysto rzeczywiste i odpowiada mu hermitowski sygnał zespolony. Tę symetrię wykorzystują funkcje hfft, które zakładają, że operujemy w dziedzinie czasu [math]\frac{n}{2}+1[/math] punktami zespolonymi i odpowiadającymi im w dziedzinie częstości [math]n[/math] punktami rzeczywistymi.

Zadanie

import sounddevice as sd
sd.play(syg, Fs)
  • wykreśl widmo mocy tego sygnału tj. kwadrat modułu transformaty Fouriera (jest to syg. rzeczywisty więc używamy funkcji z rodziny rfft)


przy okazji
wczytywanie pliku wav, tak aby był on tablicą numpy można zrobić jak w poniższym kodzie, wówczas w Fs mamy częstośc próbkowania sygnału audio, a w macierzy syg mamy sygnał (może być stereo: lewy i prawy kanał)
from scipy.io.wavfile import read
import sounddevice as sd
(Fs, syg) = read('plik.wav')
L = syg[:,0]
sd.play(syg[:,0], Fs)

Co musimy z tego zapamiętać?

  • Analiza Fourierowska polega na rzutowaniu sygnału na sygnały próbne: sinusoidy o znanych częstościach
  • Ze współczynników transformaty można odtworzyć oryginalny sygnał: mamy więc dwie reprezentacje sygnału:
    • w dziedzinie czasu
    • w dziedzinie częstości
    • między reprezentacjami przechodzimy za pomocą transformaty Fouriera i odwrotnej transformaty Fouriera
  • reprezentacja częstościowa jest na ogół zespolona: widmo mocy to kwadrat modułu transformaty Fouriera

Analiza_sygnałów_-_ćwiczenia/Fourier