Ćwiczenia 3

Analiza_sygnałów_-_ćwiczenia/Fourier_3

Twierdzenie o splocie

Przypomnijmy sobie poznane na wykładzie Twierdzenie o splocie:

[math] g(t)=\left(s * h\right)(t)\quad \Rightarrow \quad G(f)=S(f)\cdot H(f) [/math]

To twierdzenie działa też w drugą stronę:

[math] G(f)=\left(S * H\right)(f)\quad \Rightarrow \quad g(t)=s(t)\cdot h(t) [/math]

W praktyce oznacza to tyle, że jeśli w jednej dziedzinie jakieś dwa sygnały przez siebie przemnożymy, to w drugiej dziedzinie transformaty tych sygnałów splatają się. Własność ta ma bardzo ważne konsekwencje, np. przy estymacji widma skończonego fragmentu sygnału. Dlaczego?

Wyobraźmy sobie, że mamy nieskończenie długi sygnał. Oprócz niego mamy też funkcję, która jest niezerowa tylko na skończonym odcinku. Funkcję taką będziemy nazywać oknem. Pobranie fragmentu sygnału można wyobrazić sobie jako efekt pomnożenia badanego sygnału przez okno. Ta operacja mnożenia w dziedzinie czasu, w dziedzinie częstości odpowiada splotowi widma sygnału z widmem okna. Aby uzyskać sygnał o skończonej długości odrzucamy wyzerowane odcinki. FFT widzi taki skończony odcinek jako periodyczne przedłużenie. Widać, że w praktyce estymowania widma zawsze mamy do czynienia z widmami będącymi splotem widma sygnału i widma okna, ponieważ zawsze pracujemy z sygnałami o skończonej długości.

Zastosowania

• pozwala na zamianę splotu na mnożenie
• daje wgląd w okienkowanie
• łatwiej można zrozumieć działanie filtrów

Okienka

Wstęp

• Jakie są dostępne okna w scipy.signal? Proszę odwiedzić stronę z dokumentacją.
• Jak wyglądają te okienka?

from scipy import signal
w = signal.bartlett(10)
plt.plot(w,'.')
plt.show()

Badanie własności okien

Zadanie 1: Własności różnych okien

Widmo amplitudowe w skali dB sygnału [math]s[/math], jego transformata Fouriera [math]S[/math]:

[math] P_[dB] = 20*np.log10( abs(S) ) [/math]

Proszę zwrócić uwagę na charakterystyczne elementy: szerokość piku głównego, szybkość zanikania listków bocznych, zera.

• A jakie własności ma poznane na poprzednich zajęciach okno prostokątne? Najprościej można je zrobić tak: window = np.ones(N). Proszę wykonać analogiczne rysunki jak dla okna Bartletta.

• Proszę porównać przebiegi czasowe i własności widmowe pozostałych okien dostępnych w scipy.signal, m.in.:
  • blackman(M)
  • hamming(M)
  • hanning(M)
  • kaiser(M, beta)

Przydatna może być tu funkcja do obliczania widma (proszę przeanalizować kod i w razie wątpliwości poprosić o objaśnienie):

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


from numpy.fft import rfft, rfftfreq
import numpy as np
import pylab as py
import scipy.signal as ss

def widmo_dB(s, N_fft, F_samp):
    S = rfft(s,N_fft)
    S_dB = 20 * np.log10(np.abs(S))
    F = rfftfreq(N_fft, 1.0/F_samp)
    return (S_dB,F)
    
def rysuj_wlasnosci(s, okno, N_fft, F_samp,tytul):
    okno = okno/np.linalg.norm(okno)
    s_okienkowany = s* okno
    
    py.figure()
    py.subplot(2,2,1)
    py.plot(t,okno)
    py.title(tytul)

    py.subplot(2,2,2)
    S_okienka_zera, skalaF_zera = widmo_dB(okno, 1024, Fs)
    S_okienka, skalaF = widmo_dB(okno, len(s), Fs)
    
    py.plot(skalaF_zera,S_okienka_zera)
    py.plot(skalaF,S_okienka,'r.-')
    py.title('widmo okienka')
    py.ylim((-90,25))  
    
    py.subplot(2,2,3)
    py.plot(t,s_okienkowany)
    py.title('sygnał okienkowany')
    
    py.subplot(2,2,4)
    S_okienkowany_zera, skalaF_zera = widmo_dB(s_okienkowany, 1024, Fs)
    S_okienkowany, skalaF = widmo_dB(s_okienkowany,len(s), Fs)
    py.plot(skalaF_zera,S_okienkowany_zera)
    py.plot(skalaF,S_okienkowany,'r.-')
    py.plot([f,f],[-90,25],'r')
    py.ylim((-90,25))
    py.title('widmo sygnału okienkowanego')

  

#zad2 Okienko i szerokość prążka w widmie
    
Fs = 100
dt = 1/Fs
N_fft =1024

f = 10.2
T = 1

t = np.arange(0,T,dt)
s = np.sin(2*np.pi*f*t)
# blackman(M), hamming(M) hanning(M) kaiser(M, beta)

M = len(s)
okno = np.ones(M)
rysuj_wlasnosci(s, okno, N_fft, Fs,'prostokąt')

okno = ss.hanning(M)#
rysuj_wlasnosci(s, okno, N_fft, Fs,'hanning')

okno = ss.hamming(M)#
rysuj_wlasnosci(s, okno, N_fft, Fs,'hamming')

okno = ss.blackman(M)#
rysuj_wlasnosci(s, okno, N_fft, Fs,'blackman')

Zadanie 3: Wpływ okienkowania na wykrywalność składowych o różnej amplitudzie

• Wygeneruj sygnał będący sumą dwóch sinusoid o fazie 0, czasie trwania [math]T=1[/math] s, i częstości próbkowania [math]Fs=100[/math] Hz. Jedna niech ma częstość [math]f=10.2[/math] Hz. Częstość drugiej, w kolejnych zapuszczeniach skryptu proszę zmieniać od 11,4 do 16,4 Hz z krokiem co 0,5 Hz.
• Jaki jest wpływ okienek na możliwości rozróżnienia dwóch częstości?
• Proszę powtórzyć iteracje dla przypadku gdy druga z sinusoid ma 10-krotnie niższą amplitudę.

Okienkować możemy też w częstości

Zadanie 4: Konstrukcja prostego filtru

• Proszę wykreślić widmo, a następie obliczyć i wykreślić odwrotną transformatę Fouriera dla sygnału rzeczywistego zadanego w dziedzinie częstości przez (poniższy kod definiuje sygnał w dodatniej części widma):

S = np.zeros(128)
S[21] = 1

• Proszę wykreślić okienko prostokątne w częstości, pozycja 1 w binach częstości od 21 do 30, a następnie wykreślić odpowiadający mu sygnał w dziedzinie czasu.

• Koncepcyjnie najprostszym filtrem jest operacja przemnożenia widma sygnału przez okienko w dziedzinie częstości a następnie zastosowanie odwrotnej transformaty Fouriera. Co stanie się z sygnałem w dziedzinie czasu, a co w dziedzinie częstości, jeśli chcielibyśmy przefiltrować go pasmowo przy użyciu filtra pasmowo przepustowego z poprzedniego punktu?

Zadanie 5: Porównanie działania najprostszych filtrów fft na biały szum

• Proszę wytworzyć sygnał będący białym szumem o długości 256 próbek.
• Policz transformatę rzeczywistą tego sygnału.
• Przygotuj dwa okna w dziedzinie częstości (w części dodatniej) o długości 129:
  • okno prostokątne z jedynkami w punktach 20:30
  • okno Gaussowskie np.exp(-((F-25.0)/5)**2)
• Przemnóż widmo sygnału przez każde z okien i oblicz odwrotną transformatę Fouriera.
• Zilustruj wyniki rysunkiem opisanym przez poniższą tabelkę:

okno prostokątne w częstości	okno Gaussowskie w częstości
przebieg czasowy odwrotnej tr. Fouriera okna Prostokątnego	przebieg czasowy odwrotnej tr. Fouriera okna Gaussowskiego
przebieg czasowy szumu	przebieg czasowy szumu
Przebieg czasowy odwrotnej tr. Fouriera iloczynu widma szumu i okna prostokątnego	Przebieg czasowy odwrotnej tr. Fouriera iloczynu widma szumu i okna Gaussowskiego

Analiza_sygnałów_-_ćwiczenia/Fourier_3