Ćwiczenia 2 2

Zadania

Współczynniki obliczane przez FFT

Najprostsza sytuacja: Badamy współczynniki zwracane przez fft dla sinusoid o różnych częstościach. Wykorzystujemy funkcje do generacji sygnałów napisane na poprzednich zajęciach.

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

• Proszę kolejno wygenerować sinusoidy o długości 1s próbkowaną 128Hz i częstościach 1,10,30, 64 i 0 Hz. Dla tych sinusoid proszę policzyć transformaty Fouriera i wykreślić zarówno sygnały jak i wartość bezwzględne otrzymanych współczynników.
  • Jak wyglądają otrzymane wykresy?
  • Czy coś szczególnego dzieje się dla częstości 0 i 64Hz? Czy w tych skrajnych przypadkach faza sygnału ma wpływ na wynik transformaty?

• Proszę wygenerować sygnał delta położony w sekundzie 0,5 na odcinku czasu o długości 1s próbkowany 128Hz. Dla takiego sygnału proszę policzyć transformatę Fouriera i wykreślić zarówno sygnały jak i wartość bezwzględne otrzymanych współczynników.
  • Jak wygląda transformata funkcji delta? Jakie częstości w sobie zawiera?

*

Rekonstrukcja sygnału ze współczynników

Reprezentacja sygnałów w dziedzinie częstości jest dualna do reprezentacji w dziedzinie czasu. To znaczy, że jedną reprezentację można przekształcić w drugą. Do przejścia z dziedziny czasu do częstości używaliśmy transformaty Fouriera (zaimplemantowanej w fft). Przejścia z dziedziny częstości do czasu dokonujemy przy pomocy odwrotnej transformaty Fouriera (zaimplementowanej jako ifft. Mając (zespolone) współczynniki w dziedzinie częstości dla pewnego sygnału,możemy odzyskać jego przebieg czasowy.

• Proszę wygenerować sygnał [math]s(t) = \sin(2\pi t \cdot 1)+\sin\left(2 \pi t \cdot 3+\frac{\pi}{5}\right) [/math] o długości 2,5s próbkowany 10Hz, obliczyć jego transformatę Fouriera, a następnie zrekonstruować przebieg czasowy. Sygnał oryginalny i zrekonstruowany wykreślić na jednym rysunku. Korzystając z naszych funkcji generację sygnału można zapisać tak:

(s1,t) = sin(f=1,T=2.5, Fs= 10, phi =0) #generuję sygnał
(s2,t) = sin(f=3,T=2.5, Fs= 10, phi =np.pi/5)
s = s1 + s2

• Dla porównania proszę zrekonstruować sygnał korzystając jawnie z postaci odwrotnej transformaty Fouriera danej wzorem:

[math]x(n)=\frac{1}{N} \sum_{k=0}^{N-1} S[k] \cdot \exp\left(2j \pi n\frac{ k}{N}\right)[/math]

Zwróćmy uwagę, że n w powyższym wzorze to indeksy punktów czasu, zatem wiążą się one z czasem zwracanym przez nasze funkcje następująco:

n=np.arange(0,len(t),1).

zatem kod rekonstruujący może być np. taki:

n=np.arange(0,len(t),1)
s_rekonstrukcja = np.zeros(len(t))
for k in range(0,N):
	s_rekonstrukcja += 1.0/N * S[k]*np.exp((2j*np.pi*n*k)/N)

*

Rekonstrukcja na dłuższym odcinku czasu:

Korzystając z wyników poprzedniego polecenia proszę porównać wykres funkcji zadanych wzorami dla odcinak czasu o długości 5s i rekonstrukcję sygnału o długości 5s otrzymanego ze współczynników Fouriera obliczonych dla 2,5s odcinak czasu.

• Na jakim odcinku czasu sygnały są zgodne?
• Co dzieje się na dalszym odcinku czasu?

*

Efekt nieciągłości funkcji

• Wygenerować sinusoidę o następujących własnościach: f=10 Hz, T=1, Fs=100 Hz, i fazie = 1;
• Przy pomocy subplotów proszę sporządzić rysunek zgodnie z ponższym opisem:
  • subplot(2,2,1): przebieg sygnału w czasie
  • subplot(2,2,2): moduł jego transformaty Fouriera (narysować za pomocą funkcji py.stem, skalę na osi częstości można uzyskać wywołując funkcję F = FFT.fftfreq(len(s), 1/Fs) gdzie s to nasz sygnał),
  • subplot(2,2,3): Proszę wykreślić trzykrotnie periodycznie powielony oryginalny sygnał. Można go skonstruować wywołując funkcję: s_period = np.concatenate((s,s,s)).
  • subplot(2,2,4): moduł transformaty Fouriera s_period (narysować za pomocą funkcji py.stem, skalę na osi częstości można uzyskać wywołując funkcję F = FFT.fftfreq(3*len(s), 1/Fs)

• Powtórz te same kroki dla sinusa o częstości 10.3 Hz.

Pytania:

1. Czym różnią się przedłużenia sinusoidy 10Hz od sinusoidy 10.3Hz? Proszę zwrócić uwagę na miejsca sklejania sygnałów.Porównaj z wynikami otrzymanymi w zagadnieniu Rekonstrukcja na dłuższym odcinku czasu.
2. Skąd bierze się widoczna różnica w widmie sinusoidy 10Hz i 10.3Hz?

*

Długość sygnału a rozdzielczość widma FFT

Z dotychczasowych rozważań o transformacie Fouriera ograniczonych w czasie sygnałów dyskretnych wynika, że w widmie reprezentowane są częstości od [math]-F_N[/math] do [math]F_N[/math] gdzie [math]F_N[/math] to częstości Nyquista. Dostępnych binów częstości jest N - tyle samo ile obserwowanych punktów sygnału. Zatem zwiększenie długości sygnału w czasie poprawia "rozdzielczość" reprezentacji częstotliwościowej sygnału.

Załóżmy, że dysponujemy jedynie sekwencją N próbek pewnego sygnału. Rozważymy teraz jakie można przyjąć strategie przedłużania tego sygnału w celu zwiększenia gęstości binów częstotliwościowych i jakie te strategie mają konsekwencje.

Przedłużenie przez periodyczną replikację

Rozważając poprzednie przykłady zauważyliśmy, że FFT "widzi" sygnał tak jakby to była nieskończona periodyczna replikacja fragmentu sygnału podanego na wejście. Zatem najbardziej naturalną formą przedłużenia sygnału może wydawać się postępowanie zgodnie z tym sposobem widzenia i dołożenie do sygnału kolejnych segmentów zawierających kopie analizowanego fragmentu sygnału. Zbadajmy empirycznie efekty takiego podejścia.

• Proszę wytworzyć i wykreślić sinusoidy o częstościach 15 i 20 Hz, trwające 0.1s i próbkowane 100Hz.
  • Ile próbek sygnału dostajemy?
  • Proszę obliczyć transformatę Fouriera. Ile punktów w widmie amplitudowym otrzymujemy? Dla jakich cżęstości mamy biny?
  • Proszę skonstruować sygnał będący złożeniem 10 kopii oryginalnego sygnału.

Jest to 10-krotnie dłuższy fragment sygnału, którego wersję nieskończoną opisują współczynniki obliczone przez transformatę Fouriera. Transformata policzona z takiego wydłużonego odcinak ma 10-krotnie więcej binów. W szczególności zawiera on bin 15Hz. Proszę porównać wykresy widma amplitudowego otrzymanego dla wersji krótkiej i przedłużonej sygnału. Czy wynik jest zaskakujący? Jak go można zrozumieć?

* Sygnał przedłużony periodycznie nie zawiera żadnej dodatkowej informacji względem pojedynczego okresu!

Przedłużanie sygnału zerami

Inną popularną metodą na zwiększanie ilości binów w transformacie Fouriera jest przedłużanie sygnału zerami (zero-padding). Jest to szczególny przypadek następującego podejścia: Nasz "prawdziwy" sygnał jest długi. Oglądamy go przez prostokątne okno, które ma wartość 1 na odcinku czasu dla którego próbki mamy dostępne i 0 dla pozostałego czasu (więcej o różnych oknach będzie na kolejnych zajęciach). W efekcie możemy myśleć, że oglądany przez nas sygnał to efekt przemnożenia "prawdziwego" sygnału przez okno. Efekty takiego przedłużania proszę zbadać przy użyciu poniższego kodu.

• Jak można zinterpretować wyniki tego eksperymentu w świetle twierdzenia o splocie?

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

import pylab as py
import numpy as np
import numpy.fft as FFT

 
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
	'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania
	Domyślnie wytwarzany jest sygnał reprezentujący 
	1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
	'''
 
	dt = 1.0/Fs
	t = np.arange(0,T,dt)
	s = np.sin(2*np.pi*f*t + phi)
	return (s,t)	
	
(s1,t) = sin(f = 15.0, T =0.1, Fs = 100, phi = 0)
(s2,t)= sin(f = 20.0, T =0.1, Fs = 100, phi = 0)
s=s1+s2
py.clf()
py.subplot(2,2,1)
py.plot(t,s)
py.subplot(2,2,2)
S = FFT.fft(s)
F = FFT.fftfreq(len(s),0.01)
py.stem(F,np.abs(S)/len(S))
py.xlim((-50,50))
py.ylim((0,0.7))
z= np.zeros(len(s))
py.subplot(2,2,3)
s_period = np.concatenate((s,z,z,z,z,z,z,z,z,z))
t_period = np.arange(0,len(s_period)/100.0,0.01)
py.plot(t_period,s_period)

py.subplot(2,2,4)
S_period = FFT.fft(s_period)
F_period = FFT.fftfreq(len(s_period),0.01)
py.stem(F_period,np.abs(S_period)/len(S))
py.stem(F,np.abs(S)/len(S),linefmt='r-', markerfmt='ro')
py.xlim((-50,50))
py.ylim((0,0.7))
py.show()