Ćwiczenia 2

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

• Dokumentacja numpy: http://docs.scipy.org/doc/numpy/reference/routines.fft.html
• numpy implementuje algorytm: "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 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

• 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:

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