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

Z Brain-wiki
 
(Nie pokazano 29 wersji utworzonych przez 4 użytkowników)
Linia 1: Linia 1:
−
=== Zadania ===
+
[[Analiza_sygnałów_-_ćwiczenia]]/Fourier_2
  
−
==== Współczynniki obliczane przez FFT ====
 
−
Najprostsza sytuacja: Badamy współczynniki zwracane przez <tt>fft</tt> dla sinusoid o różnych częstościach. Wykorzystujemy funkcje do generacji sygnałów napisane na poprzednich zajęciach.
 
  
−
* 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.
+
=Odwracalność transformaty=
−
** Jak wygląda transformata funkcji delta? Jakie częstości w sobie zawiera?
+
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 <tt>fft</tt>). Przejścia z dziedziny częstości do czasu dokonujemy przy pomocy odwrotnej transformaty Fouriera (zaimplementowanej jako <tt>ifft</tt>. Mając (zespolone) współczynniki w dziedzinie częstości dla pewnego sygnału, możemy odzyskać jego przebieg czasowy.
−
*
+
===Zadanie 1===
 +
* 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,5 s próbkowany 100 Hz, obliczyć jego transformatę Fouriera za pomocą <tt>fft</tt>, a następnie zrekonstruować przebieg czasowy za pomocą <tt>ifft</tt>. Sygnał oryginalny i zrekonstruowany wykreślić na jednym rysunku.  ''Uwaga: funkcja ifft zwraca wektor liczb zespolonych. Sprawdź jaka jest jegeo część urojona. Na wykresie rekonstrukcji przedstaw jego część rzeczywistą.''
 +
 
 
<!--
 
<!--
−
<source lang=python>
 
−
# -*- coding: utf-8 -*-
 
−
# skrypt implemetujący polecenia w sekcji
 
−
# Współczynniki obliczane przez FFT
 
−
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)
 
−
 
−
def delta(t0=0.5, T=1  ,Fs = 128):
 
−
dt = 1.0/Fs
 
−
t = np.arange(0,T,dt)
 
−
d = np.zeros(len(t))
 
−
d[np.ceil(t0*Fs)]=1
 
−
return (d,t)
 
−
 
−
# do pytania 1
 
−
(s,t) = sin(f=2,T=1,Fs=10,phi=np.pi/2)
 
−
S = FFT.fft(s)
 
−
print S
 
−
# na wydruku należy zauważyć, że wsółczynniki są zespolone i parami sprzężone
 
−
 
−
# do pytania 2
 
−
(s,t) = sin(f=1,T=1,Fs=128,phi=np.pi/3) # f=0,10,30,64
 
−
S = FFT.fft(s)
 
−
py.figure(1)
 
−
py.subplot(2,1,1)
 
−
py.plot(t,s,'o')
 
−
py.subplot(2,1,2)
 
−
py.plot(np.abs(S),'o')
 
−
 
−
# do pytania 3
 
−
(s,t) = delta(t0=0.5, T=1  ,Fs = 128)
 
−
 
−
S = FFT.fft(s)
 
−
py.figure(2)
 
−
py.subplot(2,1,1)
 
−
py.plot(t,s,'o')
 
−
py.subplot(2,1,2)
 
−
py.plot(np.abs(S),'o')
 
−
</source>
 
−
-->
 
−
 
−
=====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 <tt>fft</tt>). Przejścia z dziedziny częstości do czasu dokonujemy przy pomocy odwrotnej transformaty Fouriera (zaimplementowanej jako <tt>ifft</tt>. 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:
 
−
<source lang = python>
 
−
(s1,t) = sin(f=1, T=2.5, Fs=10, phi=0)
 
−
(s2,t) = sin(f=3, T=2.5, Fs=10, phi=np.pi/5)
 
−
s = s1 + s2
 
−
</source>
 
 
* Dla porównania proszę zrekonstruować sygnał korzystając jawnie z postaci odwrotnej transformaty Fouriera danej wzorem:
 
* 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>
 
:<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:
 
: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:
−
:<tt>n=np.arange(0,len(t),1)</tt>.
+
:<tt>n = np.arange(0,len(t),1)</tt>.
 +
 
 
zatem kod rekonstruujący może być np. taki:
 
zatem kod rekonstruujący może być np. taki:
 
<source lang = python>
 
<source lang = python>
−
n=np.arange(0,len(t),1)
+
n = np.arange(0,len(t),1)
 
s_rekonstrukcja = np.zeros(len(t))
 
s_rekonstrukcja = np.zeros(len(t))
 
for k in range(0,N):
 
for k in range(0,N):
 
s_rekonstrukcja += 1.0/N * S[k]*np.exp((2j*np.pi*n*k)/N)
 
s_rekonstrukcja += 1.0/N * S[k]*np.exp((2j*np.pi*n*k)/N)
 
</source>
 
</source>
−
*
+
-->
 +
 
 
<!--
 
<!--
 
<source lang = python>
 
<source lang = python>
Linia 196: Linia 133:
  
 
-->
 
-->
−
==== Efekt nieciągłości funkcji ====
+
 
 +
=Badanie rozdzielczości sygnałami testowymi=
 +
* Poniżej będziemy zajmować się sygnałami rzeczywistymi, więc stosujemy funkcje z rodziny Real FFT:
 +
https://docs.scipy.org/doc/numpy/reference/routines.fft.html
 +
 
 +
* W poniższych przykładach jako widmo będziemy rozumieli widmo amplitudowe, tzn wartość bezwzględną ze współczynników szeregu Fouriera.
 +
 
 +
==Widmo sinusoidy i delty==
 +
Najprostsza sytuacja: Badamy współczynniki zwracane przez <tt>fft</tt> dla sinusoid o różnych częstościach.
 +
===Zadanie 2===
 +
* Proszę kolejno wygenerować sinusoidy o długości 1s próbkowaną 32Hz i częstościach 1,10, 16  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 16Hz? Czy w tych skrajnych przypadkach faza sygnału ma wpływ na wynik transformaty?
 +
===Zadanie 3===
 +
 
 +
* 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?
 +
 
 +
<!--
 +
<source lang=python>
 +
# -*- coding: utf-8 -*-
 +
# skrypt implemetujący polecenia w sekcji
 +
# ==Badanie rozdzielczości sygnałami testowymi==
 +
import pylab as py
 +
import numpy as np
 +
from numpy.fft import rfft, rfftfreq
 +
 +
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)
 +
 +
def delta(t0=0.5, T=1  ,Fs = 128):
 +
dt = 1.0/Fs
 +
t = np.arange(0,T,dt)
 +
d = np.zeros(len(t))
 +
d[np.ceil(t0*Fs)]=1
 +
return (d,t)
 +
 +
# do pytania 1
 +
(s,t) = sin(f=2,T=1,Fs=10,phi=np.pi/2)
 +
S = rfft(s)
 +
print( S)
 +
# na wydruku należy zauważyć, że wsółczynniki są zespolone i parami sprzężone
 +
 
 +
# do pytania 2
 +
Fs = 30
 +
for f in (0,1,10,16):
 +
    py.figure()
 +
    (s,t) = sin(f=f,T=1,Fs=Fs,phi=np.pi/3) # f=0,10,30,64
 +
    S = rfft(s)
 +
    F = rfftfreq(s.size, 1/Fs)
 +
 
 +
    py.subplot(2,1,1)
 +
    py.plot(t,s,'o-')
 +
    py.title(str(f))
 +
    py.subplot(2,1,2)
 +
    py.plot(F, np.abs(S),'o')
 +
    py.show()
 +
 
 +
# do pytania 3
 +
py.figure()
 +
 
 +
(s,t) = delta(t0=0.5, T=1  ,Fs = 128)
 +
 
 +
S = rfft(s)
 +
F = rfftfreq(s.size, 1/Fs)
 +
py.figure(2)
 +
py.subplot(2,1,1)
 +
py.plot(t,s,'o')
 +
py.subplot(2,1,2)
 +
py.plot(F, np.abs(S),'o')
 +
</source>
 +
-->
 +
 
 +
=== Efekt nieciągłości funkcji ===
 +
====Zadanie 4====
 
* Wygenerować sinusoidę o następujących własnościach: f=10 Hz, T=1, Fs=100 Hz, i fazie = 1;
 
* 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:
 
* 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,1): przebieg sygnału w czasie  
−
** subplot(2,2,2): moduł jego transformaty Fouriera (narysować za pomocą funkcji <tt>py.stem</tt>, skalę na osi częstości można uzyskać wywołując funkcję <tt>F = FFT.fftfreq(len(s), 1/Fs)</tt> gdzie <tt>s</tt> to nasz sygnał),
+
** subplot(2,2,2): moduł jego transformaty Fouriera (narysować za pomocą funkcji <tt>py.stem</tt> wraz zprawidłową osią częstości,
 
** subplot(2,2,3): Proszę wykreślić trzykrotnie periodycznie powielony oryginalny sygnał. Można go skonstruować wywołując funkcję: <tt>s_period = np.concatenate((s,s,s))</tt>.
 
** subplot(2,2,3): Proszę wykreślić trzykrotnie periodycznie powielony oryginalny sygnał. Można go skonstruować wywołując funkcję: <tt>s_period = np.concatenate((s,s,s))</tt>.
−
** subplot(2,2,4): moduł transformaty Fouriera <tt>s_period</tt> (narysować za pomocą funkcji <tt>py.stem</tt>, skalę na osi częstości można uzyskać wywołując funkcję <tt>F = FFT.fftfreq(3*len(s), 1/Fs)</tt>
+
** subplot(2,2,4): moduł transformaty Fouriera <tt>s_period</tt> (narysować za pomocą funkcji <tt>py.stem</tt> wraz zprawidłową osią częstości
  
 
* Powtórz te same kroki dla sinusa o częstości 10.3 Hz.
 
* Powtórz te same kroki dla sinusa o częstości 10.3 Hz.
 
Pytania:
 
Pytania:
−
# Czym różnią się przedłużenia sinusoidy 10 Hz od sinusoidy 10.3 Hz? Proszę zwrócić uwagę na miejsca sklejania sygnałów. <!--Porównaj z wynikami otrzymanymi w zagadnieniu '''Rekonstrukcja na dłuższym odcinku czasu'''.-->
+
# Czym różnią się przedłużenia sinusoidy 10 Hz od sinusoidy 10.3 Hz? Proszę zwrócić uwagę na miejsca sklejania sygnałów.
 +
<!--Porównaj z wynikami otrzymanymi w zagadnieniu '''Rekonstrukcja na dłuższym odcinku czasu'''.-->
 
# Skąd bierze się widoczna różnica w widmie sinusoidy 10 Hz i 10.3 Hz?  
 
# Skąd bierze się widoczna różnica w widmie sinusoidy 10 Hz i 10.3 Hz?  
−
*
+
 
 
<!--
 
<!--
 
<source lang = python>
 
<source lang = python>
Linia 215: Linia 235:
 
import pylab as py
 
import pylab as py
 
import numpy as np
 
import numpy as np
−
import numpy.fft as FFT
+
from numpy.fft import rfft, rfftfreq
  
 
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
 
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
+
    '''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  
+
    Domyślnie wytwarzany jest sygnał reprezentujący  
−
1 sekundę sinusa o częstości 1Hz i zerowej fazie próbkowanego 128 Hz
+
    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)
 
−
 
−
(s,t) = sin(f = 10.3, T =1, Fs = 100, phi = 1)
 
  
−
py.subplot(2,2,1)
+
    dt = 1.0/Fs
−
py.plot(t,s)
+
    t = np.arange(0,T,dt)
−
py.subplot(2,2,2)
+
    s = np.sin(2*np.pi*f*t + phi)
−
S = FFT.fft(s)
+
    return (s,t)   
−
F = FFT.fftfreq(len(s),0.01)
+
   
−
py.stem(F,np.abs(S)/len(S))
+
for f in (10, 10.3):
 +
    py.figure()
 +
    (s,t) = sin(f = f, T =1, Fs = 100, phi = 1)
 +
   
 +
    py.subplot(2,2,1)
 +
    py.plot(t,s)
 +
    py.title(f)
 +
    py.subplot(2,2,2)
 +
    S = rfft(s)
 +
    F = rfftfreq(len(s),0.01)
 +
    py.stem(F,np.abs(S)/len(S))
 +
   
 +
    py.subplot(2,2,3)
 +
    s_period = np.concatenate((s,s,s))
 +
    t_period = np.arange(0,3,0.01)
 +
    py.plot(t_period,s_period)
 +
   
 +
    py.subplot(2,2,4)
 +
    S_period = rfft(s_period)
 +
    F_period = rfftfreq(len(s_period),0.01)
 +
    py.stem(F_period,np.abs(S_period)/len(S_period))
 +
</source>
  
−
py.subplot(2,2,3)
+
-->
−
s_period = np.concatenate((s,s,s))
 
−
t_period = np.arange(0,3,0.01)
 
−
py.plot(t_period,s_period)
 
  
−
py.subplot(2,2,4)
+
===Długość sygnału a rozdzielczość widma FFT ===
−
S_period = FFT.fft(s_period)
+
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.
−
F_period = FFT.fftfreq(len(s_period),0.01)
 
−
py.stem(F_period,np.abs(S_period)/len(S_period))
 
−
</source>
 
  
−
-->
+
* jaki dostęp między binami częstotliwości mamy dla 1 s sygnału próbkowanego 10Hz?
 +
* jaki dostęp między binami częstotliwości mamy dla 1 s sygnału próbkowanego 100Hz?
 +
* jaki dostęp między binami częstotliwości mamy dla 1 s sygnału próbkowanego 1000Hz?
 +
* jaki dostęp między binami częstotliwości mamy dla 10 s sygnału próbkowanego 10Hz?
 +
* jaki dostęp między binami częstotliwości mamy dla 100 s sygnału próbkowanego 10Hz?
  
−
==== Długość sygnału a rozdzielczość widma FFT ====
+
Zatem zwiększenie długości sygnału w czasie poprawia "rozdzielczość"  reprezentacji częstotliwościowej sygnału.  
−
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.
 
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ę =====
 
=====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.  
 
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.  
Linia 264: Linia 295:
 
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ć?
 
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!
 
  * Sygnał przedłużony periodycznie nie zawiera żadnej dodatkowej informacji względem pojedynczego okresu!
 +
-->
 
<!--
 
<!--
 
<source lang = python>
 
<source lang = python>
Linia 307: Linia 339:
 
-->
 
-->
  
−
=====Przedłużanie sygnału zerami =====
+
====Przedłużanie sygnału ====
−
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.  
+
=====Przedłużanie przez cykliczne powielenie=====
 +
Zobaczmy co się stanie jesli przedłużymy sygnał prze jego periodyczne przedłużenie. Efekty takiego przedłużania proszę zbadać przy użyciu poniższego kodu:
 +
<source lang = python>
 +
# -*- coding: utf-8 -*-
 +
import pylab as py
 +
import numpy as np
 +
from numpy.fft import rfft, rfftfreq
 +
 
 +
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)
 +
 
 +
Fs =100
 +
T =0.1
 +
 
 +
(s,t) = sin(f = 10.0, T=T, Fs=Fs)
 +
 
 +
 
 +
py.figure()
 +
py.subplot(2,2,1)
 +
py.plot(t,s)
 +
py.subplot(2,2,2)
 +
S = rfft(s)/len(s)
 +
F = rfftfreq(len(s),1/Fs)
 +
py.stem(F,np.abs(S))
 +
 
 +
z= np.zeros(len(s))
 +
py.subplot(2,2,3)
 +
n = 10
 +
s_period = np.hstack(n*(s,))# n razy powtarzamy s
 +
t_period = np.arange(0,T*n,1/Fs)
 +
py.plot(t_period,s_period)
 +
 
 +
py.subplot(2,2,4)
 +
S_period = rfft(s_period)/len(s)
 +
F_period = rfftfreq(len(s_period),1/Fs)
 +
py.stem(F_period,np.abs(S_period)/(len(s_period)))
 +
py.stem(F,np.abs(S),linefmt='r-', markerfmt='ro')
 +
 
 +
py.show()
 +
</source>
 +
 
 +
=====Przedłużanie zerami=====
 +
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ć:
 +
* dla sygnału sinusoidalnego o dł. 0.1s i częstości 10Hz próbkowanego 100 Hz
 +
* dla sygnału sinusoidalnego o dł. 0.1s i częstości 22Hz próbkowanego 100 Hz
 +
* dla sygnału będącego suma dwóch powyższych
 
* Jak można zinterpretować wyniki tego eksperymentu w świetle [[Twierdzenia_o_splocie_i_o_próbkowaniu_(aliasing)#Twierdzenie_o_splocie|twierdzenia o splocie]]?
 
* Jak można zinterpretować wyniki tego eksperymentu w świetle [[Twierdzenia_o_splocie_i_o_próbkowaniu_(aliasing)#Twierdzenie_o_splocie|twierdzenia o splocie]]?
 +
<!--
 
<source lang = python>
 
<source lang = python>
 
# -*- coding: utf-8 -*-
 
# -*- coding: utf-8 -*-
 +
 +
"""
 +
Created on Fri Oct 21 15:51:33 2016
 +
 +
@author: admin
 +
"""
  
 
import pylab as py
 
import pylab as py
 
import numpy as np
 
import numpy as np
−
import numpy.fft as FFT
+
from numpy.fft import rfft, rfftfreq
  
 
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
 
def sin(f = 1, T = 1, Fs = 128, phi =0 ):
Linia 327: Linia 419:
 
s = np.sin(2*np.pi*f*t + phi)
 
s = np.sin(2*np.pi*f*t + phi)
 
return (s,t)
 
return (s,t)
−
+
 
−
(s1,t) = sin(f = 15.0, T =0.1, Fs = 100, phi = 0)
+
Fs =100
−
(s2,t)= sin(f = 20.0, T =0.1, Fs = 100, phi = 0)
+
T =0.1
 +
n = 10
 +
(s1,t) = sin(f = 10.0, T=T, Fs=Fs)
 +
(s2,t)= sin(f = 22.0, T =0.1, Fs = Fs)
 
s=s1+s2
 
s=s1+s2
−
py.clf()
+
 
 +
py.figure()
 
py.subplot(2,2,1)
 
py.subplot(2,2,1)
 
py.plot(t,s)
 
py.plot(t,s)
 
py.subplot(2,2,2)
 
py.subplot(2,2,2)
−
S = FFT.fft(s)
+
S = rfft(s)/len(s)
−
F = FFT.fftfreq(len(s),0.01)
+
F = rfftfreq(len(s),1/Fs)
−
py.stem(F,np.abs(S)/len(S))
+
py.stem(F,np.abs(S))
−
py.xlim((-50,50))
+
 
−
py.ylim((0,0.7))
 
 
z= np.zeros(len(s))
 
z= np.zeros(len(s))
 
py.subplot(2,2,3)
 
py.subplot(2,2,3)
 
s_period = np.concatenate((s,z,z,z,z,z,z,z,z,z))
 
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)
+
t_period = np.arange(0,n*T,1/Fs)
 +
py.plot(t_period,s_period)
 +
 
 +
py.subplot(2,2,4)
 +
S_period = rfft(s_period)/len(s)
 +
F_period = rfftfreq(len(s_period),1/Fs)
 +
py.stem(F_period,np.abs(S_period))
 +
py.stem(F,np.abs(S),linefmt='r-', markerfmt='ro')
 +
 
 +
py.show()
 +
 
 +
 
 +
py.figure()
 +
py.subplot(2,2,1)
 +
py.plot(t,s)
 +
py.subplot(2,2,2)
 +
S = rfft(s)/len(s)
 +
F = rfftfreq(len(s),1/Fs)
 +
py.stem(F,np.abs(S))
 +
 
 +
z= np.zeros(len(s))
 +
py.subplot(2,2,3)
 +
 
 +
s_period = np.hstack((s,s,s,s,s,s,s,s,s,s))# n razy powtarzamy s
 +
t_period = np.arange(0,T*n,1/Fs)
 
py.plot(t_period,s_period)
 
py.plot(t_period,s_period)
  
 
py.subplot(2,2,4)
 
py.subplot(2,2,4)
−
S_period = FFT.fft(s_period)
+
S_period = rfft(s_period)/len(s_period)
−
F_period = FFT.fftfreq(len(s_period),0.01)
+
F_period = rfftfreq(len(s_period),1/Fs)
−
py.stem(F_period,np.abs(S_period)/len(S))
+
py.stem(F_period,np.abs(S_period))
−
py.stem(F,np.abs(S)/len(S),linefmt='r-', markerfmt='ro')
+
py.stem(F,np.abs(S),linefmt='r-', markerfmt='ro')
−
py.xlim((-50,50))
+
 
−
py.ylim((0,0.7))
 
 
py.show()
 
py.show()
 
</source>
 
</source>
 +
-->
 +
 +
=Co musimy z tego zapamiętać?=
 +
* Sygnał może być reprezentowany w dziedzine czasu lub w dziedzinie częstości
 +
* Jak wyglada widmo delty?
 +
* Jak wygląda widmo sinusa, którego całkowita ilość okresów mieści się w badanym fragmencie, a jak jeśli niecałkowita?
 +
* Jak długość sygnału wpływa na rozdzielczość widma?
 +
* Jakie częstości występują w widmie sygnału periodyzowanego cyklicznie?
 +
* Jaki efekt daje przedłużanie zerami?
 +
 +
[[Analiza_sygnałów_-_ćwiczenia]]/Fourier_2

Aktualna wersja na dzień 15:18, 10 lis 2016

Analiza_sygnałów_-_ćwiczenia/Fourier_2


Odwracalność transformaty

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.

Zadanie 1

  • 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,5 s próbkowany 100 Hz, obliczyć jego transformatę Fouriera za pomocą fft, a następnie zrekonstruować przebieg czasowy za pomocą ifft. Sygnał oryginalny i zrekonstruowany wykreślić na jednym rysunku. Uwaga: funkcja ifft zwraca wektor liczb zespolonych. Sprawdź jaka jest jegeo część urojona. Na wykresie rekonstrukcji przedstaw jego część rzeczywistą.



Badanie rozdzielczości sygnałami testowymi

  • Poniżej będziemy zajmować się sygnałami rzeczywistymi, więc stosujemy funkcje z rodziny Real FFT:

https://docs.scipy.org/doc/numpy/reference/routines.fft.html

  • W poniższych przykładach jako widmo będziemy rozumieli widmo amplitudowe, tzn wartość bezwzględną ze współczynników szeregu Fouriera.

Widmo sinusoidy i delty

Najprostsza sytuacja: Badamy współczynniki zwracane przez fft dla sinusoid o różnych częstościach.

Zadanie 2

  • Proszę kolejno wygenerować sinusoidy o długości 1s próbkowaną 32Hz i częstościach 1,10, 16 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 16Hz? Czy w tych skrajnych przypadkach faza sygnału ma wpływ na wynik transformaty?

Zadanie 3

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


Efekt nieciągłości funkcji

Zadanie 4

  • 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 wraz zprawidłową osią częstości,
    • 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 wraz zprawidłową osią częstości
  • Powtórz te same kroki dla sinusa o częstości 10.3 Hz.

Pytania:

  1. Czym różnią się przedłużenia sinusoidy 10 Hz od sinusoidy 10.3 Hz? Proszę zwrócić uwagę na miejsca sklejania sygnałów.
  2. Skąd bierze się widoczna różnica w widmie sinusoidy 10 Hz i 10.3 Hz?


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.

  • jaki dostęp między binami częstotliwości mamy dla 1 s sygnału próbkowanego 10Hz?
  • jaki dostęp między binami częstotliwości mamy dla 1 s sygnału próbkowanego 100Hz?
  • jaki dostęp między binami częstotliwości mamy dla 1 s sygnału próbkowanego 1000Hz?
  • jaki dostęp między binami częstotliwości mamy dla 10 s sygnału próbkowanego 10Hz?
  • jaki dostęp między binami częstotliwości mamy dla 100 s sygnału próbkowanego 10Hz?

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żanie sygnału

Przedłużanie przez cykliczne powielenie

Zobaczmy co się stanie jesli przedłużymy sygnał prze jego periodyczne przedłużenie. Efekty takiego przedłużania proszę zbadać przy użyciu poniższego kodu: 

# -*- coding: utf-8 -*-
import pylab as py
import numpy as np
from numpy.fft import rfft, rfftfreq

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)	

Fs =100	
T =0.1

(s,t) = sin(f = 10.0, T=T, Fs=Fs)


py.figure()
py.subplot(2,2,1)
py.plot(t,s)
py.subplot(2,2,2)
S = rfft(s)/len(s)
F = rfftfreq(len(s),1/Fs)
py.stem(F,np.abs(S))

z= np.zeros(len(s))
py.subplot(2,2,3)
n = 10
s_period = np.hstack(n*(s,))# n razy powtarzamy s
t_period = np.arange(0,T*n,1/Fs)
py.plot(t_period,s_period)

py.subplot(2,2,4)
S_period = rfft(s_period)/len(s)
F_period = rfftfreq(len(s_period),1/Fs)
py.stem(F_period,np.abs(S_period)/(len(s_period)))
py.stem(F,np.abs(S),linefmt='r-', markerfmt='ro')

py.show()
Przedłużanie zerami

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

  • dla sygnału sinusoidalnego o dł. 0.1s i częstości 10Hz próbkowanego 100 Hz
  • dla sygnału sinusoidalnego o dł. 0.1s i częstości 22Hz próbkowanego 100 Hz
  • dla sygnału będącego suma dwóch powyższych
  • Jak można zinterpretować wyniki tego eksperymentu w świetle twierdzenia o splocie?

Co musimy z tego zapamiętać?

  • Sygnał może być reprezentowany w dziedzine czasu lub w dziedzinie częstości
  • Jak wyglada widmo delty?
  • Jak wygląda widmo sinusa, którego całkowita ilość okresów mieści się w badanym fragmencie, a jak jeśli niecałkowita?
  • Jak długość sygnału wpływa na rozdzielczość widma?
  • Jakie częstości występują w widmie sygnału periodyzowanego cyklicznie?
  • Jaki efekt daje przedłużanie zerami?

Analiza_sygnałów_-_ćwiczenia/Fourier_2

Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=Ćwiczenia_2_2&oldid=5976”