Podstawy analizy i przetwarzania sygnałów - laboratorium¶
Celem zajęć jest praktyczne zapoznanie się z wybranymi, prostymi metodami przetwarzania i analizy sygnałów. Sygnałami, stanowiącymi przedmiot przedstawianych analiz są sygnały audio. Aby możliwe było zbudowanie aplikacji korzystającej z narzędzi programowych, pozwalających na realizację operacji dedykowanych analizie i przetwarzaniu sygnałów, konieczne jest 'załadowanie' odpowiednich bibliotek. Początek programu do analizy będzie więc zawierać deklaracje importu odpowiednich bibliotek, przy czym oprócz typowo stosowanych w pracy z Pythonem biblitek 'numpy' i 'matplotlib', dodatkowymi modułami oferującymi funkcje specjalizowane będzie 'pyaudio' oraz 'scipy.io.wavfile':
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
import math
import pyaudio #sudo apt-get install python-pyaudio
Ćwiczenie 1. Wczytanie i odegranie pliku audio¶
W ramach ćwiczenia pliki audio będą wczytywane z plików, przy czym używanym formatem reprezentacji sygnału w pliku będzie format 'wav' - plik zawiera nagłówek, informujący o atrybutach (takich jak liczba kanałów, częstotliwość próbkowania itd.) oraz 'surowe' dane. Wczytanie pliku to wykonanie metody 'read' zaimportowanego obiektu wavfile:
fs, dane = wavfile.read('./up1.wav')
Argumentem metody odczytu jest ścieżka do pliku, zaś wynikiem jej wykonania są dwa elementy: pierwszym jest częstotliwość próbkowania z jaką zarejestrowany został sygnał, zaś drugim, jest treść sygnału. Każda próbka jest reprezentowana dwoma bajtami, jako liczba całkowita (typu np.int16), tak więc zmienna 'dane' to wektor zmiennych dwubajtowych.
Wczytany plik może zostać odsłuchany z wykorzystaniem narzędzi programowych dostarczonych przez bibliotekę 'pyaudio'. Aby zapewnić możliwość wielokrotnego wykorzystania w programie funkcjonalności odegrania pliku, w przykładach towarzyszących modułowi została zdefiniowana funkcja 'play(.)', którą należy umieścić w kodzie programu (dla czytelności struktury programu, warto, by uczynić to w początkowej sekcji kodu programu):
def play(Dane,fs) :
PyAudio = pyaudio.PyAudio
p = PyAudio()
stream = p.open(format = p.get_format_from_width(Dane.itemsize),
channels = 1,
rate = fs,
output = True)
stream.write(Dane,num_frames=np.size(Dane)) # konieczny 2 arg, inazcej ucina
stream.stop_stream()
stream.close()
p.terminate()
Wywołanie funkcji z argumentami: danymi, zawartymi w wektorze 'dane' i częstotliwością próbkowania 'fs' spowoduje odegranie pliku
play(dane,fs)
Podanie niewłaściwej częstotliwości próbkowania powoduje zmianę szybkości odtwarzania danych, co prowadzi do zauważalnego zniekształcenia dżwięku:
play(dane,int(0.7*fs))
Zastosowana wartość współczynnika skalującego spowalnia odtwarzanie, jednocześnie w zauważalny sposób zmienia brzmienie głosów mówiących osób. Głosy stają się 'niższe', a więc dokonany zabieg jest równoważny procedurze liniowego zmniejszenia częstotliwości wszystkich komponentów sygnału mowy.
Dyskretna transformacja Fouriera (DFT)¶
Dziedziną większości procedur przetwarzania i analizy sygnałów audio jest właśnie dziedzina częstotliwości, w szczególności, gdy rozważane sygnały są tworzone jako wynik superpozycji wielu komponentów okresowych o różnych częstotliwościach (jak na przykład sygnał mowy, sygnał EKG, dźwięki instrumentów muzycznych i inne). Dyskretna Transformacja Fouriera (DFT) to metoda dekompozycji sygnału na komponenty harmoniczne. DFT jest odpowiednikiem ciągłej Transformacji Fouriera, a więc przekształcenia:
$G(\omega) = \int^{\infty}_{-\infty} f(x) e^{-j\omega x}dx$
adresowanym dla przypadku przetwarzania skończonych zbiorów danych dyskretnych (a więc jedynej możliwej do reprezentacji w komputerze formy sygnałów). Zakładając, że dany jest sygnał w postaci N-elementowego zbioru próbek $x_k, k \in [0... N-1]$, pobieranych w jednakowych odstępach czasu $T_s$ (a więc z częstotliwością próbkowania $f_s$), DFT to przekształcenie pozwalające na ocenę podobieństwa rozważanego sygnału do zbioru zespolonych funkcji harmonicznych o dyskretnych częstotliwościach z przedziału od 0 Hz ('składowa stała') do $f_s/2$.
$X_l = \sum_{k=0}^{N-1} x_k e^{j2\pi f_l k} = \sum_{k=0}^{N-1} x_k e^{j2\pi \frac{l k}{N}}$
Okres najwolniej zmiennej harmonicznej, do której porównywany jest sygnał to cały czas trwania tego sygnału (a więc $N\cdot T_s$), okres kolejnej jest dwukrotnie mniejszy (w czasie trwania sygnału następują dwa pełne okresy tej sinusoidy), okres kolejnej jest trzy razy mniejszy itd., aż do okresu równego $2\cdot T_s$. W konsekwencji, częstotliwości komponentów harmonicznych $f_l$ zmieniają się od $\frac{1}{N\cdot T_s}$, przez $\frac{2}{N\cdot T_s}$, $\frac{3}{N\cdot T_s}$ ... aż do $\frac{1}{2\cdot T_s}$. Oznacza to, że widmo przebiegu rzeczywistego złożonego z $N$-próbek będzie miało $N/2$ różnych komponentów $X_l$.
Przykład: wyznaczenie DFT danego sygnału x¶
Załóżmy, że dany jest sygnał x będący ważoną sumą trzech sinusoid: x1, x2 i x3 o różnych częstotliwościach: $x = 2x_1 + x_2 + 0.5x_3$, z któregp pobrano osiem próbek:
x1 = np.sin(np.arange(32)*2*np.pi/32) # f1 = fmin = 1/(NTs), Ts=1
x2 = np.sin(np.arange(32)*2*np.pi*2/32) # f2 = 1/(N/2 ) = 2*fmin
x3 = np.sin(np.arange(32)*2*np.pi*3/32) # f3 = 1/N/3 = 3*fmin
x = 2*x1 + x2 + 0.5*x3
plt.plot(range(32),x,linestyle='--', label='2*x1+x2')
# punkty wybrane do wektora 'dane'
s = x[0:32:4]
plt.plot(np.arange(0,32,4),s,linestyle='None', marker='o')
plt.plot(range(32),2*x1,linewidth=0.5, label='2*x1')
plt.plot(range(32),x2,linewidth=0.5, label='x2')
plt.plot(range(32),0.5*x3,linewidth=0.5, label='0.5*x3')
plt.legend()
Próbki, które będą poddane analizie (zaznaczone na wykresie znacznikami) są umieszczone w zmiennej 'dane' i mają wartości:
x[0:32:4]
DFT ciągu 'dane', umieszczone w zmiennej X (patrz kod poniżej), to wektor zawierający liczby zespolone, określające osiem współczynników Fouriera, oceniających podobieństwo tego wektora do funkcji harmonicznych o kolejnych częstotliwościach: X[0] - składowa stała, X[1] - częstotliwość $1/(N\cdot T_1)$ (co odpowiada okresowi funkcji harmonicznej o okresie $T_1 = N = 8$), X[2] - częstotliwość $2/(N\cdot T_1)$ (okres $T_2 = 4 = T_1/2$) i X[3] - częstotliwość $3/(N\cdot T_1)$ (okres $T_3 = 8/3 = T_1/3$). Wszystkich częstotliwości jest tyle, ile jest elementów ciągu wejściowego, czyli 8. Ponieważ największa częstotliwość komponentu harmonicznego, jaka może zostać reprezentowana w ciągu to N/2 (najmniejszy okres to $T=2 = T_1/ (N/2)$) komponenty widma począwszy od N/2+1 do N, nie odpowiadają częstotliwościom większym, a są równe sprzężonym wartościom zwierciadlanego odbicia wcześniejszego fragmentu widma, gdzie osią symetrii jest próbka o indeksie N/2.
X = np.fft.fft(s)
plt.plot(range(8),np.absolute(X),linestyle='None', marker='o')
Wykres modułu widma pokazuje 'skład' analizowanego sygnału - jak widać poprawnie odtworzone zostały propoecje między składowymi sygnału - największą wartość ma komponent o najniższej częstotliwości (1/8 N), dwa razy mniejszą ma kompoment o częstotliwości dwa razy większej (2/8 N) a cztery razy mniejszą ma komponent o częstotliwości trzy razy większej (3/8 N).
Jeżeli wyznaczymy DFT dla sygnału x (a więc, mającego 32-elementy), a nie dla jego ośmiu próbek (zawartych w ciągu s), uzyskamy widmo zawierające 32 komponenty, ale jego struktura pozostanie podobna - sygnał zawiera trzy komponenty częstotliwościowe, o określonych wcześniej, różnych amplitudach.
N=len(x)
X = np.fft.fft(x)
plt.plot(range(N),np.absolute(X),linestyle='None', marker='o')
Ćwiczenie 2: modyfikacja sygnału w dziedzinie widma¶
Widmo sygnału jest często wykorzystywane jako podstawa analizy sygnału lub podstawa do wykonywania różnych manipulacji na sygnale. Na przykład, jeżeli zarejestrowany sygnał jest sygnałem EKG, to komponent harmoniczny o największej amplitudzie będzie odpowiadał tętnu: automatyczna analiza tętna realizowana w 'smartwatch-ach' to detekcja sygnału (mogą to być zmiany ćiśnienia na ściankę zegarka, mogą być sygnały elektryczne rejestrowane przez zegarek), obliczenie jego DFT i określenie częstotliwości komponentu o najwyższej amplitudzie.
Celem manipulacji dokonywanych w dziedzinie widma jest najczęściej ekstrakcja lub usuwanie wybranych komponentów częstotliwościowych widma sygnału. Powodem wydzielania określonych komponentów widma jest np. separacja sygnału użytecznego (niosącego treść) od sygnałów używanych do transmisji treści (tzw. demodulacja). Powodem usuwania wybranych komponentów jest często poprawianie jakości sygnału (jeżeli usuwane komponenty stanowią elementy zakłócające).
Zmiany komponentów widma sygnału muszą uwzględniać charakter widma - jeżeli zmianie podlega komponent 'i'-ty, taka sama modyfikacja musi dotyczyć komponentu 'N-i'. Rekonstrukcja z widma sygnału oryginalnego lub zmodyfikowanego to wyznaczenie odwrotnej dyskretnej transformacji Fouriera, czyli IDFT. Przykładowe manipulacje widmem sygnału x oraz ich efekty zostały przedstawione poniżej.
# rekonstrukcja sygnału oryginalnego
xr = np.fft.ifft(X).real
plt.figure()
plt.plot(range(N),xr,label='xr')
# usunięcie (wyzerowanie) pierwszego komponentu widma i rekonstrukcja
X[1]=X[N-1]=0
x23 = np.fft.ifft(X).real
plt.plot(range(N),x23,label='x23')
# podwojenie trzeciego komponentu widma i rekonstrukcja
X[3]*=2
X[N-3]*=2
x23 = np.fft.ifft(X).real
plt.plot(range(N),x23,label='x23x2')
plt.legend()
Ćwiczenie 3: usuwanie zakłóceń z sygnału¶
Sygnały podlegające analizom bardzo często zawierają oprócz składnika będącego przedmniotem zainteresowania (sygnału 'użytecznego') również inne komponenty. Ponieważ komponenty te utrudniają percepcję sygnału użytecznego lub przeprowadzanie jego analizy, uznawane są za zakłócenie (często nazywane 'szumem'), a jednym z istotnych elementów przetwarzania sygnałów jest eliminacja lub redukcja zakłóceń. Sygnały zakłócające często mają charakter addydtywny ponieważ stanowią typowe 'tło' dla rejestracji sygnału użytecznego (np. dźwięki z otoczenia występujące podczas nagrywanej rozmowy, sygnały elektromagnetyczne istniejące w środowisku rejestracji EKG itp.). Przykładowe sygnały: mowa i szum są przedstawione poniżej:
fs, dane = wavfile.read('./m_sam1.wav')
play(dane,fs)
plt.plot(range(len(dane)),dane)
fs, szum = wavfile.read('./szum.wav') # częstotliwość próbkowania jest taka sama, więc użyta jest ta sama zmienna
play(szum,fs)
plt.plot(range(len(szum)),szum)
Widma sygnału i szumu różnią się zasadniczo - o ile widmo sygnału jest skomplikowane i wskazuje na bardzo złożone proporcje zawartości różnych komponentów częstotliwościowych, o tyle w przypadku wczytanego szumu, wszystkie komponenty mają duże amplitudy o niewielkim zróżnicowaniu.
plt.subplot(121)
plt.plot(range(len(dane)),np.absolute(np.fft.fft(dane)))
plt.subplot(122)
plt.plot(range(len(szum)),np.absolute(np.fft.fft(szum)))
Jeżeli zakłócenia są sygnałami okresowymi, które dodają się do sygnału użytecznego (są addytywne), wykorzystanie reprezentacji sygnału w dziedzinie Fouriera może pozwolić na prostą metodę ich usuwania: należy zidentyfikować w widmie komponenty zakłócenia, a następnie je wyzerować lub zredukować. To postępowanie nazywa się 'filtracją' - 'filtr' to algorytm działający jak sito, eliminujące z widma komponenty o określonych częstotliwościach. Załóżmy, że dany jest następujący, silnie zakłócony sygnał, z którego chcemy wydzielić sygnał użyteczny:
fs, dane_s = wavfile.read('./m_sam2n.wav')
play(dane_s,fs)
plt.plot(range(len(dane_s)),dane_s)
Widmo wczytanego sygnału ma postać:
widmo = np.absolute(np.fft.fft(dane_s))
plt.plot(range(len(widmo)),widmo)
Jak widać, w zakresie wysokich częstotliwości (od ok. 3 kHz) w sygnale znajduje się zbiór komponentów o bardzo dużych i niemal jednakowych amplitudach - taka struktura sugeruje, że stanowią one zakłócenie. Usnięcie tego zakłócenia może zostać dokonane przez wyzerowanie wszystkich komponentów widma odpowiadających szumowi (od ok. 3.5 kHz do 8 kHz), a wynik takiej operacji jest następujący
fs, dane_f = wavfile.read('./m_sam2nf.wav')
play(dane_f,fs)