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

Z Brain-wiki
 
(Nie pokazano 33 wersji utworzonych przez 3 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 38: Linia 78:
 
</math>
 
</math>
 
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 ) - \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.
+
===Rekonstrukcja sygnałów: ===
 +
Na podstawie wyników poprzedniego zadania proszę napisać program, który demonstruje jaki jest wynik składania coraz większej liczby 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ą 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ą 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ń 07:11, 20 paź 2023

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 ) - \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 liczby 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