Ćwiczenia 2

Transformata Fouriera

Szereg Fouriera

Ciągłą funkcję okresową (sygnał okresowy) można przedstawić w postaci szeregu Fouriera — 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:

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

szereg można zwinąć do szeregu funkcji sinusoidalnych.

Zadanie

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

• 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

Przykład

Oblicz analitycznie, 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:

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

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

# -*- coding: utf-8 -*-

import pylab as py
import numpy as np

def skladnik_modul_t(n ):
	'''n - który element szeregu'''
	if n==0:
		c_n = 0.5
	else:
		c_n = 1/(np.pi**2 * n**2) * (-1  + np.cos(np.pi*n) )
	return c_n
	
t = np.arange(-4,4,0.1)
f = np.zeros(len(t))

N = 10
c = []
for n in range(-N,N+1):
	cn = skladnik_modul_t(n)
	c.append(cn)
	f += cn * np.exp(-1j*2*np.pi*t*n/2.0)
py.plot(t,f)
for n in range(-N,N+1):
	print '%(c).4f'%{'c':c[n+N]}
py.show()

Zadanie

Analogicznie do powyższego przykładu proszę znaleźć i zbadać szereg Fouriera dla funkcji:

[math] f(t) = |cos(t)|[/math]

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

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 1s, częstości 2Hz próbkowaną 10Hz. 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?

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.