%matplotlib inline
import numpy as np
import scipy.signal as sig
from scipy.io import wavfile
import matplotlib.pyplot as plt
from util import set_style
set_style()
# częstotliwość próbkowania
fs = 48000
Analiza częstotliwościowa¶
Cyfrowe sygnały mogą być analizowane w dziedzinie czasu oraz w dziedzinie częstotliwości. Algorytmy działające w dziedzinie czasu operują bezpośrednio na próbkach cyfrowego sygnału. Analiza częstotliwościowa sygnału, stosowana w wielu operacjach cyfrowego przetwarzania sygnałów, polega na przekształceniu sygnału z dziedziny czasu do dziedziny częstotliwości. Widmo sygnału (ang. spectrum) to reprezentacja częstotliwościowa sygnału, czesto mówi się więc o analizie widmowej.
Jean-Baptiste Joseph Fourier odkrył, że dowolny sygnał okresowy można przedstawić w postaci szeregu Fouriera, czyli sumy sygnałów trygonometrycznych (kosinusów i sinusów) o różnych amplitudach i częstotliwościach. Operację przekształcenia sygnału z dziedziny czasu do dziedziny częstotliwości nazywa się przekształceniem (transformacją) Fouriera (Fourier transform). Wynik tej operacji, czyli częstotliwościową reprezentację sygnału, nazywa się transformatą Fouriera. W przypadku sygnałów dyskretnych (cyfrowych) mówimy o dyskretnym przekształceniu Fouriera (DFT - Discrete Fourier Transform). W praktyce nie oblicza się transformaty z definicji, stosuje się algorytm szybkiego przekształcenia Fouriera (FFT - Fast Fourier Transform).
Aby obliczyć FFT sygnału w języku Python, możemy skorzystać z funkcji fft zawartej w module numpy.fft. Obliczmy FFT dla sygnału sinusoidalnego o częstotliwości 1000 Hz.
n = np.arange(2048)
sinus = np.sin(2 * np.pi * n * 1000 / fs)
widmo = np.fft.fft(sinus)
W podręcznikach opisuje się często algorytm FFT typu radix-2, w którym rozmiar transformaty musi być potęgą dwójki. Algorytm wykorzystywany przez NumPy jest elastyczny. Naszybciej działa wtedy, gdy liczba próbek jest potęgą dwójki. Będzie jednak działać również dla innych rozmiarów transformaty, jedynie obliczenia mogą trwać dłużej. Najwięcej obliczeń jest potrzebnych wtedy, gdy liczba próbek jest dużą liczbą pierwszą, wtedy konieczne jest liczenie DFT z definicji.
Obliczanie FFT¶
W praktyce korzystamy z gotowych procedur obliczających FFT, ale warto pokzać w jaki sposób można samodzielnie obliczyć FFT. Funkcja napisana w języku Python będzie znacznie wolniejsza od tej z modułu NumPy, ale w podobny sposób można zaimplementować tę operację w systemie, w którym procedura obliczania FFT nie jest dostępna. Pokażemy tutaj najprostszy, klasyczny algorytm Cooleya-Tukeya typu radix-2, z podziałem czasowym (DIT - decimation in time), operujący na ciągach wartości sygnału o długości będącej potęgą liczby 2.
Mamy $N$ próbek sygnału dyskretnego ($N$ musi być potęgą dwójki). Dzielimy ten ciąg na dwa ciągi: $x_p(n)$, $x_n(n)$, składające się z próbek sygnału $x(n)$ o indeksach odpowiednio parzystych (0, 2, 4, ...) i nieparzystych (1, 3, 5, ...). Obliczamy transformaty Fouriera tych ciągów: $X_p(n)$, $X_n(n)$, wywołując rekurencyjnie procedurę obliczania FFT. Następnie składamy transformatę całego ciągu próbek, zgodnie z zależnością:
$$X(k) = X_p(k) + W_k^N X_n(k)$$ $$X(k+N/2) = X_p(k) - W_k^N X_n(k)$$ $$W_k^N = e^{-2\pi jk / N}$$
dla $0 \le k < N/2$. Tak więc obliczanie FFT polega na kolejnym dzieleniu próbek na ciągi parzyste i nieparzyste, aż do otrzymania dwupunktowych ciągów, dla których $X(0) = x(0) + x(1)$, $X(1) = x(0) - x(1)$. Następnie transformaty ciągów parzystych i nieparzystych są składane w całość.
def oblicz_fft(x):
'''Oblicza FFT dla sygnału x o długości będącej potęgą 2.'''
N = len(x)
if N == 2:
return np.array([x[0] + x[1], x[0] - x[1]], dtype=np.complex)
Xp = oblicz_fft(x[::2])
Xn = oblicz_fft(x[1::2])
w = np.exp(-2 * np.pi * 1j * np.arange(N // 2) / N)
X = np.hstack((Xp + w * Xn, Xp - w * Xn))
return X
x = np.random.random(16)
print('Własna funkcja FFT:')
print(oblicz_fft(x))
print('numpy.fft.fft:')
print(np.fft.fft(x))
Widzimy, że otrzymaliśmy wynik zgodny z tym podawanym przez funkcję z modułu NumPy.
Widmo amplitudowe¶
Na podstawie N próbek sygnału cyfrowego obliczamy N wartości zespolonego widma. Ponieważ widmo sygnału jest sygnałem zespolonym, zazwyczaj analizujemy moduł widma - jest to widmo amplitudowe (magnitude spectrum):
$$\left|H(k)\right|=\sqrt{\Re(H(k))^2+\Im(H(k))^2}$$
Moduł zespolonej wartości obliczamy funkcją np.abs. Wykonajmy wykres widma.
widmo_amp = np.abs(widmo)
plt.plot(widmo_amp)
plt.xlabel('indeks widma')
plt.ylabel('amplituda widma')
plt.title('Pełne widmo sygnału sinusoidalnego f=1 kHz')
plt.show()
Na podstawie N próbek sygnału otrzymujemy N wartości widmowych, pokrywających zakres od 0 do $f_s$, gdzie $f_s$ jest częstotliwością próbkowania, a $f_s/2$ jest częstotliwością Nyquista. Odstęp między dwoma punktami widma wynosi więc $f_s/N$, w naszym przypadku: 23,4375 Hz.
Sygnał sinusoidalny zawiera tylko jedną składową widmową o ustalonej częstotliwości. Zatem widmem sinusa, obserwowanym w zakresie od 0 Hz do częstotliwości Nyquista, jest jeden prążek widmowy. Na wykresie pełnego widma widzimy dwa prążki. Jeżeli sygnał poddawany transformacji Fouriera jest rzeczywisty, wtedy zawsze otrzymamy dwie kopie widma, druga kopia będzie lustrzanym odbiciem pierwszej (dokładnie: będzie jej zespolonym sprzężeniem). Często potrzebujemy tylko jednej kopii widma. Funkcja np.fft.rfft oblicza widmo rzeczywistego sygnału. Funkcja np.fft.fftfreq oblicza częstotliwość dla każdego punktu widmowego, argumentami funkcji są: liczba wartości widma (długość transformaty) oraz okres próbkowania. Ponadto skalujemy widmo dzieląc je przez połowę długości transformaty, aby wartości widmowe reprezentowały energię poszczególnych składowych (przez połowę długości, ponieważ energia sygnału rozkłada się równomiernie na dwie części widma, a my obserwujemy tylko jedną z nich).
widmo_amp = np.abs(np.fft.rfft(sinus)) / 1024
f = np.fft.rfftfreq(2048, 1/fs)
plt.plot(f, widmo_amp)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda widma')
plt.title('Widmo "rzeczywiste" sygnału sinusoidalnego 1 kHz')
plt.show()
Na podstawie $N$ próbek sygnału otrzymujemy $(N/2)+1$ wartości zespolonego widma jeżeli $N$ jest parzyste, $(N-1)/2$ wartości jeżeli $N$ jest nieparzyste. Pierwsza wartość transformaty jest składową stałą, czyli sumą wartości sygnału (jest ona zawsze rzeczywista). Jeżeli $N$ jest parzyste, ostatnia wartość wyniku reprezentuje składową dla częstotliwości Nyquista (również jest rzeczywista).
Obliczmy i wykreślmy widmo dla sygnału będącego sumą trzech sinusów o częstotliwościach: 1000, 3000, 5000 kHz i amplitudach: 0,5; 0,3; 0,2. Widmo jest wykreślone dla częstotliwości do 10000 Hz.
trzysin = (0.5 * np.sin(2 * np.pi * n * 1000 / fs) +
0.3 * np.sin(2 * np.pi * n * 3000 / fs) +
0.2 * np.sin(2 * np.pi * n * 5000 / fs))
widmo_trzysin = np.fft.rfft(trzysin)
plt.plot(f, np.abs(widmo_trzysin) / 1024)
plt.xlim(0, 10000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda widma')
plt.title('Widmo sumy trzech sinusów')
plt.show()
Zgodnie z oczekiwaniami, widzimy trzy prążki na częstotliwościach składowych sinusów.
Widmo amplitudowe na powyższych wykresach zostało pokazane w skali liniowej. Bardzo często wykreśla się widmo w skali logarytmicznej, obliczając 20 logarytmów dziesiętnych z amplitudy widma. Jednostką poziomu widma jest decybel (widmo w skali decybelowej).
plt.plot(f, 20 * np.log10(np.abs(widmo_trzysin) / 1024))
plt.xlim(0, 10000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda widma [dB]')
plt.title('Widmo sumy trzech sinusów - skala decybelowa')
plt.show()
Widmo fazowe¶
Widmo fazowe przedstawia fazę poszczególnych składowych widmowych. Fazę zespolonego widma oblicza się jako:
$$arg(H(k))=arctan\left(\frac{\Im(H(k))}{\Re(H(k))}\right)$$
Do obliczenia kąta fazowego używamy funkcji np.angle. Wynik jest podawany w radianach. Przykład dla sygnału sinusoidalnego 1000 Hz:
widmo_faz = np.angle(np.fft.rfft(sinus))
plt.plot(f, widmo_faz)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('faza [rad]')
plt.title('Widmo fazowe sygnału sinusoidalnego')
plt.show()
Przy analizie widma fazowego, pamiętajmy że faza jest cykliczna, zawija się wokół wartości $2\pi$. Jeżeli nie chcemy mieć przeskoków fazy, możemy rozwinąć fazę za pomocą funkcji np.unwrap. W powyższym przypadku rozwijanie fazy nie zmieni wyniku.
Widmowa gęstość mocy, periodogram¶
Badając widmo sygnału, pokazujemy tutaj widmo amplitudowe, czyli moduł widma, w skali decybelowej. W wielu dziedzinach stosuje się zamiast tego widmową gęstość mocy (PSD, power spectral density, a wynik analizy nazywa się periodogramem. Metoda ta ma zastosowanie szczególnie do sygnałów losowych oraz sygnałów z dodanym szumem. Tradycyjna metoda obliczania periodogramu operuje w dziedzinie czasu, obecnie periodogram oblicza się jako iloczyn widma sygnału i zespolonego sprzężenia widma. Dla sygnałów rzeczywistych, periodogram jest równy kwadratowi modułu widma. Przedstawia on rozklad mocy widmowej sygnału na jednostkę częstotliwości.
Do obliczenia periodogramu służy funkcja scipy.signal.periodogram. Istotne parametry to kolejno: analizowany sygnał, częstotliwość próbkowania, typ funkcji okna, rozmiar FFT. Ponadto parametr scaling decyduje o tym, czy periodogram jest skalowany do wartości mocy na jednostkę częstotliwości ('density'), czy nie ('spectrum'). Funkcja zwraca: wartości częstotliwości oraz wartości widmowej gęstości mocy.
fper, Pxx = sig.periodogram(trzysin, fs, 'hamming', 2048, scaling='density')
plt.semilogy(fper, Pxx)
plt.xlim(0, 10000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('widmowa gęstość mocy')
plt.title('Periodogram sumy trzech sinusów')
plt.show()
Rozdzielczość częstotliwościowa i czasowa analizy¶
W poprzednich przykładach rozmiar transformaty wynosił 2048. Jakie znaczenie ma ta wartość? Na podstawie $N$ próbek sygnału otrzymujemy $N$ wartości widma w zakresie od 0 do $f_s$. Zatem odstęp pomiędzy dwoma wartościami widma wynosi:
$$d_f = \frac{f_s}{N}$$
Wartość $d_f$ nazywa się rozdzielczością częstotliwościową transformaty Fouriera. Jest to najmniejsza wykrywalna różnica pomiędzy dwoma składowymi widmowymi. Dwie składowe częstotliwościowe sygnału odległe od siebie o mniej niż $d_f$ będą reprezentowane w widmie przez tę samą wartość.
Dla $f_s$ = 48 kHz, rozdzielczość częstotliwościowa zmienia się ze zwiększaniem rozmiaru transformaty następująco:
| N | 128 | 256 | 512 | 1024 | 2048 | 4096 |
|---|---|---|---|---|---|---|
| $d_f$ [Hz] | 375,00 | 187,50 | 93,75 | 46,88 | 23,44 | 11,72 |
Aby lepiej zilustrować znaczenie rozdzielczości czestotliwościowej, spróbujmy wykryć częstotliwość prążka widmowego dla pojedynczego sinusa 1000 Hz używając różnych rozmiarów transformaty. Funkcja np.argmax wyszukuje indeks maksymalnej wartości widma.
y = sinus = np.sin(2 * np.pi * np.arange(4096) * 1000 / fs)
for N in (128, 256, 512, 1024, 2048, 4096):
w = np.abs(np.fft.rfft(y[:N]))
indeks_prazka = np.argmax(w)
czestotliwosc_prazka = indeks_prazka * (fs / N)
blad = 100 * np.abs(czestotliwosc_prazka - 1000) / 1000
print('N = {:4}: częstotliwość = {:>8.3f} Hz, błąd = {:5.2f}%'.format(N, czestotliwosc_prazka, blad))
Widać że na skutek skończonej rozdzielczości FFT nie dostajemy nigdy dokładnej wartości częstotliwości prążka. Raz jest ona za mała, raz za duża. Jednak można zauważyć, że większy rozmiar transformaty to zawsze mniejszy błąd estymacji częstotliwości. Potrzebowaliśmy rozmiaru co najmniej 2048 aby błąd spadł poniżej 1%.
Czy zatem większy rozmiar transformaty jest zawsze lepszy? Jeżeli widmo sygnału jest niezmienne w czasie, tak jak w przypadku sinusa, wtedy tak. W praktycznych przypadkach (np. analiza sygnału mowy) widmo sygnału jest jednak zmienne w czasie. Większy rozmiar transformaty wymaga większej liczby próbek, a więc dłuższego odcinka czasowego. Czas, z którego pobieramy $N$ próbek wynosi:
$$d_t = \frac{N}{f_s}$$
Wartość $d_t$ jest rozdzielczością czasową FFT. Określa ona minimalny odstęp czasowy pomiędzy dwoma zdarzeniami w sygnale, który umożliwia rozróżnienie tych dwóch zdarzeń. Jeżeli np. w sygnale pojawią się dwa impulsy odległe od siebie o mniej niż $d_t$, to przy obliczaniu FFT z odcinka zawierającego oba impulsy zostaną one "uśrednione".
Poniższa tabela pokazuje zarówno rozdzielczość częstotliwościową, jak i czasową.
| N | 128 | 256 | 512 | 1024 | 2048 | 4096 |
|---|---|---|---|---|---|---|
| $d_f$ [Hz] | 375,00 | 187,50 | 93,75 | 46,88 | 23,44 | 11,72 |
| $d_t$ [ms] | 2,67 | 5,33 | 10,67 | 21,33 | 42,67 | 85,33 |
Rozdzielczość częstotliwościowa jest więc odwrotnością rozdzielczości czasowej. Zwiększamy rozmiar transformaty i poprawiamy rozdzielczość częstotliwościową. Ale jednocześnie pogarszamy rozdzielczość czasową. Dlatego duże rozmiary FFT stosujemy dla sygnałów o stabilnym widmie, gdy zależy nam na dokładnym rozróżnieniu składowych częstotliwościowych. Z kolei małe rozmiary FFT stosujemy przy analizie szybkozmiennego widma, godząc się na mniej dokładną analizę składowych. W przypadkach pośrednich szukamy kompromisu, rozmiary FFT równe 1024 i 2048 są bardzo często stosowane w praktycznych analizach sygnałów takich jak sygnały muzyczne spróbkowane z częstotliwością 48 lub 44,1 kHz.
Należy jeszcze wspomnieć o uzupełnianiu próbek zerami (zero padding). Możemy wziąć np. 1024 próbki sygnału, wstawić na koniec 1024 zera i policzyć transformatę z 2048 próbek. W ten sposób pozornie zwiększamy dwukrotnie rozdzielczość częstotliwościową względem tranformaty o rozmiarze 1024 bez uzupełniania zerami, zachowując rozdzielczość czasową transformaty 1024. Jest to pozorny zysk, ponieważ nie otrzymujemy w ten sposób więcej danych. W uzyskanym widmie, co druga wartość będzie odpowiadała wartościom widma transformaty o rozmiarze 1024, a pozostałe wartości widma są interpolowane. Niemniej jednak, metoda ta jest w pewnych zastosowaniach użyteczna.
Okna czasowe¶
Na zamieszczonych wcześniej wykresach widma sygnału złożonego z trzech sinusów można zauważyć, że prążek dla składowej 2000 Hz jest "ostry", podczas gdy prążki dla 1000 i 3000 Hz są rozmyte. Dlaczego tak się dzieje? Pamiętajmy, że FFT jest obliczane dla sygnału okresowego. Zakłada się, że część sygnału poddawana transformacji Fouriera stanowi okres tego sygnału. Jeżeli $N$ próbek sinusa obejmuje całkowitą wielokrotność jego okresu, to po obliczeniu FFT z tych próbek zobaczymy na wykresie jeden "ostry" prążek. Jeżeli ten warunek nie zostanie spełniony i "utniemy" okres sygnału gdzieś w punkcie pośrednim, to obliczone FFT będzie słuszne dla sygnału, którego okres wygląda dokładnie tak jak te $N$ próbek, a nie jak rzeczywisty okres sygnału. Przy zapętleniu tego bloku próbek wystąpi przeskok amplitudy między ostatnią a pierwszą próbką. W efekcie zamiast ostrego prążka zobaczymy prążek "rozlany" na boki. Zjawisko to nazywa się przeciekiem widma (spectral leakage).
Popatrzmy na przypadek, w którym każdy z tych sinusów analizujemy osobno.
sin1000 = np.sin(2 * np.pi * n * 1000 / fs)
sin2000 = np.sin(2 * np.pi * n * 2000 / fs)
sin3000 = np.sin(2 * np.pi * n * 3000 / fs)
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(sin1000)) / 1024), label='1000 Hz')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(sin2000)) / 1024), label='2000 Hz')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(sin3000)) / 1024), label='3000 Hz')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnałów sinusoidalnych')
plt.xlim(0, 5000)
plt.ylim(-60, 0)
plt.legend()
plt.show()
plt.figure()
i = np.arange(2000, 2048)
plt.plot(i, sin1000[-48:], label='1000 Hz')
plt.plot(i, sin2000[-48:], label='2000 Hz')
plt.plot(i, sin3000[-48:], label='3000 Hz')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Przebieg czasowy sygnałów sinusoidalnych')
plt.legend()
plt.show()
for fsin in (1000, 2000, 3000):
print('Częstotliwość = {} Hz: {:.3f} okresów na 2048 próbek'.format(fsin, 2048 / (fs / fsin)))
Dla częstotliwości 3000 Hz, blok 2048 próbek zawiera dokładnie 128 okresów sinusa. Dlatego na tej częstotliwości widzimy w widmie wyraźny prążek. Inaczej jest dla częstotliwości 1000 i 2000 Hz: tutaj blok 2048 próbek ucina okres w czasie jego trwania, a efektem w obu przypadkach jest przeciek widmowy.
Zjawisko przecieku widma można też wytłumaczyć następująco. Jeżeli bierzemy z sygnału $N$ próbek do obliczenia FFT, to można to potraktować jako mnożenie sygnału przez funkcję okna prostokątnego (boxcar), czyli funkcję, która ma wartość 1 dla próbek branych do FFT i 0 dla pozostałych. Widmo okna prostokątnego jest proporcjonalne do funkcji $sinc$, czyli $sin(x)/x$. Mnożeniu w dziedzinie czasu odpowiada splot w dziedzinie widma. Jeżeli okno prostokątne nie obejmuje wielokrotności okresu sygnału, rzeczywiste widmo sygnału i widmo okna prostokątnego "rozjeżdżają się" przy obliczaniu splotu. Wynikiem tego są zafalowania charakterystyki widmowej.
Jak duży problem stanowią przecieki widmowe? W przypadku sygnału sinusoidalnego - niewielki. Ale w przypadku sygnałów o złożonym widmie, przecieki od składowych sumują się z innymi składowymi widmowymi, zniekształacając ich amplitudy. Nie da się tego problemu w pełni uniknąć, ale można zredukować jego skutki stosując funkcje okna inne niż prostokątne. Zaprojektowano specjalne funkcje okna, których wspólną cechą jest amplituda bliska zeru na skraju okna. W ten sposób redukuje się efekt skoku amplitudy przy zapętleniu okna.
Spośród wielu funkcji okna, najczęściej stosowane są trzy: Hamminga (np.hamming), von Hanna ("Hanning", np.hanning) oraz Blackmana (np.blackman). Kształt tych okien pokzuje poniższy wykres. Okno Blackmana ma węższy kształt niż pozostałe. Okno Hamminga nie "dochodzi" do wartości 0, jest nieco uniesione.
plt.plot(np.hamming(2048), label='np.hamming')
plt.plot(np.hanning(2048), label='np.hanning')
plt.plot(np.blackman(2048), label='np.blackman')
plt.title('Typowe funkcje okien czasowych')
plt.legend()
plt.show()
Zastosowanie okna czasowego polega na przemnożeniu bloku próbek przez wybraną funkcję okna przed obliczeniem transformaty Fouriera. Poniżej pokazano wynik obliczania FFT dla sygnału sinusoidalnego 1000 Hz po zastosowaniu każdego z okien.
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(sin1000)) / 1024), label='prostokątne')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(np.hamming(2048) * sin1000)) / 1024), label='Hamminga')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(np.hanning(2048) * sin1000)) / 1024), label='von Hanna')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(np.blackman(2048) * sin1000)) / 1024), label='Blackmana')
plt.xlim(750, 1250)
plt.ylim(-90, 0)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału sin f=1 kHz dla różnych typów okna')
plt.legend()
plt.show()
Dla okna prostokątnego, poziom przecieków mierzony na częstotliwości 750 Hz wynosi ok. -30 dB. Zastosowanie okna Hamminga powoduje stłumienie przecieków na tej częstotliwości o ok. 20 dB. Okno von Hanna daje większe tłumienie, ok. 50 dB, ale przy samym prążku (ok. 920 Hz) tłumienie jest gorsze niż dla okna Hamminga. Wreszcie, okno Blackmana daje największe tłumienie dla 750 Hz: niecałe 60 dB.
Zatem czy okno Blackmana jest najlepsze? Popatrzmy co się dzieje dla sygnału sinusoidalnego o częstotliwości 3000 Hz, a więc w sytuacji, gdy dla okna o długości 2048 próbek nie ma przecieków. Wykres pokazuje częstotliwości wokół 3000 Hz.
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(sin3000)) / 1024), label='prostokątne')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(np.hamming(2048) * sin3000)) / 1024), label='Hamminga')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(np.hanning(2048) * sin3000)) / 1024), label='von Hanna')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(np.blackman(2048) * sin3000)) / 1024), label='Blackmana')
plt.xlim(2900, 3100)
plt.ylim(-90, 0)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału sin f=3 kHz dla różnych typów okna')
plt.legend()
plt.show()
Widać że każda ze specjalnych funkcji okna znacznie poszerza prążek względem okna prostokątnego. Okna Hamminga i von Hanna dają praktycznie identyczny wynik. Okno Blackmana powoduje największe poszerzenie prążka. Zatem większe stłumienie przecieków przez funkcję okna wiąże się z większym poszerzeniem głównego prążka. Wybór funkcji okna zależy więc od potrzeb, nie ma okna uniwersalnego i optymalnego. Jeżeli chodzi nam o jak najlepsze wyodrębnienie prążków z tła, dobre będzie okno Blackmana, natomiast jeżeli priorytetem jest rozróżnienie prążków położonych blisko siebie, lepszym wyborem może być okno Hamminga czy von Hanna.
Krótkookresowe przekształcenie Fouriera (STFT)¶
Do tej pory wykonywaliśmy analizę częstotliwościową dowolnie wybranego fragmentu sygnału. Ma to sens jeżeli wimo sygnału jest niezmienne. Ale praktyczne sygnały, np. mowa lub muzyka, mają widmo zmieniające się w czasie. Wybierając różny fragment sygnału do analizy otrzymamy różne wyniki. Aby to zilustrować, obliczymy widmo sygnału świergotowego (chirp). Jest to sygnał sinusoidalny, przy czym częstotliwość sygnału wzrasta lub maleje, liniowo lub logarytmicznie. Sygnał świergotowy jest używany np. do pomiarów urządzeń. Aby wygenerować taki sygnał, musimy sięgnąć po moduł scipy.signal i użyć funkcji chirp. Poniżej tworzymy sygnał trwający 1 sekundę, o częstotliwości wzrastającej liniowo od 20 Hz do 5000 Hz. Wykreślmy widmo dla dwóch fragmentów sygnału.
chirp = sig.chirp(np.arange(fs) / fs, 20, 1, 5000, 'linear')
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(chirp[10000:12048])) / 1024))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału chirp od próbki nr 10 000')
plt.xlim(0, 5000)
plt.ylim(-60, 0)
plt.show()
plt.figure()
plt.plot(f, 20 * np.log10(np.abs(np.fft.rfft(chirp[30000:32048])) / 1024))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału chirp od próbki nr 30 000')
plt.xlim(0, 5000)
plt.ylim(-60, 0)
plt.show()
Widzimy, że w oby przypadkach "złapaliśmy" inną chwilową częstotliwość sygnału. Nie wystarczy więc zbadać widma w jednej chwili, musimy przeanalizować cały sygnał.
Krótkookresowe przekształcenie Fouriera (STFT - Short Time Fourier Transform) polega na podzieleniu sygnału na odcinki czasowe za pomocą funkcji okna, obliczeniu widma dla każdego odcinka osobno oraz połączeniu wyników w całość. Można wyobrazić sobie okno przesuwające się po osi czasu. Wynik obliczenia STFT jest trójwymiarowy: czas – częstotliwość – amplituda. Jest on zwykle przedstawiany albo w formie dwuwymiarowego spektrogramu, albo na trójwymiarowym wykresie typu waterfall.
W praktyce bardzo często stosuje się zakładkowanie (overlapping) okien czasowych. Polega to na tym, że okno nie jest przesuwane o swoją długość, lecz o część długości, zwykle 25%, 50% lub 75%. Niektóre próbki są więc poddawane analizie w więcej niż jednym oknie. Robi się tak z dwóch powodów. Po pierwsze, redukuje się efekty brzegowe - zniekształcenia w skrajnych obszarach okna wynikające z mnożenia amplitud próbek przez wartości funkcji okna bliskie zeru. Po drugie, zwiększa się w ten sposób rozdzielczość czasową analizy.
Do obliczenia spektrogramu użyjemy funkcji spectrogram z modułu scipy.signal. Mamy również do dyspozycji funkcję scipy.signal.stft, dla której istnieje przekształcenie odwrotne scipy.signal.istft. Użyjemy okna Hamminga o rozmiarze 2048, z zakładkowaniem 75% (1536 próbek). Funkcja zwraca trzy macierze, kolejno: wartości częstotliwości, wartości czasu dla kolejnych okien oraz obliczone amplitudy widmowe.
f, t, Sxx = sig.spectrogram(chirp, fs=fs, window=np.hamming(2048), nperseg=2048, noverlap=1536,
scaling='spectrum', mode='magnitude')
plt.pcolormesh(t, f, 20 * np.log10(Sxx))
plt.xlabel('czas [s]')
plt.ylabel('częstotliwość [Hz]')
plt.title('Spektrogram sygnału chirp')
plt.ylim(0, 5000)
plt.colorbar()
plt.show()
Spektrogram przedstawiany jest następująco: na osi poziomej mamy czas, na osi pionowej - częstotliwość, natomiast amplituda widma jest przedstawiona za pomocą określonej mapy barw (przedstawionej na słupku powyżej). Mapę kolorów można wybrać za pomocą parametru cmap funkcji pcolormesh. Widzimy na wykresie liniową zmianę częstotliwości analizowanego sygnału.
Popatrzmy jeszcze na mocno powiększony fragment spektrogramu.
plt.pcolormesh(t, f, 20 * np.log10(Sxx))
plt.xlabel('czas [s]')
plt.ylabel('częstotliwość [Hz]')
plt.title('Fragment spektrogramu dla okna o rozmiarze 2048')
plt.xlim(0, 0.1)
plt.ylim(0, 200)
plt.show()
Widoczne prostokąty uwidaczniają rozdzielczość czasową (szerokość prostokąta) i częstotliwościową (wysokość prostokąta). Wszystkie elementy sygnału leżące wewnątrz jednego prostokąta "zlewają się" w jeden wynik. Popatrzmy jak wygląda widmo tego samego sygnału dla okna o rozmiarze 512.
f, t, Sxx2 = sig.spectrogram(chirp, fs=fs, window=np.hamming(512), nperseg=512, noverlap=384,
scaling='spectrum', mode='magnitude')
plt.pcolormesh(t, f, 20 * np.log10(Sxx2))
plt.xlabel('czas [s]')
plt.ylabel('częstotliwość [Hz]')
plt.title('Fragment spektrogramu dla okna o rozmiarze 512')
plt.xlim(0, 0.1)
plt.ylim(0, 200)
plt.show()
Widać że prostokąty stały się znacznie węższe (lepsza rozdzielczość czasowa), ale zarazem stały się dużo wyższe (gorsza rozdzielczość częstotliwościowa). Zmiana rozmiaru okna czasowego zmienia więc proporcje szerokości do wysokości każdego prostokąta, ale nie da się prostokąta "zmniejszyć". Te wykresy ilustrują w praktyce opisany wcześniej problem zależności obu rozdzielczości FFT od siebie.
Dokonajmy teraz analizy rzeczywistego dźwięku - jest to nagranie klarnetu odgrywającego dźwięk o wysokości $a^1$ (440 Hz). Do wczytania dźwięku użyjemy funkcji scipy.io.wavfile.read. Zwraca ona dwie wartości: częstotliwość próbkowania i tablicę zawierającą wartości próbek.
wav_fs, klarnet = wavfile.read('pliki/klarnet.wav')
f, t, Sxx = sig.spectrogram(klarnet, fs=wav_fs, window=np.hamming(2048), nperseg=2048, noverlap=1536,
scaling='spectrum', mode='magnitude')
plt.figure(figsize=(16,8
))
plt.pcolormesh(t, f, 20 * np.log10(Sxx))
plt.xlabel('czas [s]')
plt.ylabel('częstotliwość [Hz]')
plt.title('Spektrogram dźwięku klarnetu')
plt.ylim(0, 8000)
plt.colorbar()
plt.show()
Spektrogram pokazuje zmiany widma dźwięku klarnetu. W początkowej fazie dźwięku, tzw. fazie ataku, dźwięk zaczyna się tworzyć gdy muzyk wprawia powietrze wewnątrz instrumentu w drgania - widmo jest w tej fazie niestabilne. Następnie mamy fazę podtrzymania (muzyk cały czas dmucha w instrument). Widzimy harmoniczną strukturę dźwięku - częstotliwość podstawową 440 Hz i jej całkowite wielokrotności - harmoniczne. Około 1,2 s. muzyk przestaje pobudzać instrument i dźwięk zaczyna się wygaszać (faza wybrzmiewania lub zwolnienia). Wszystkie te informacje możemy odczytać z histogramu dzięki analizie STFT.
Odwrotne przekształcenie Fouriera (IFFT)¶
Wiemy już, że FFT dokonuje przekształcenia sygnału z dziedziny czasu do dziedziny częstotliwości. Mamy również odwrotne przekształcenie Fouriera (w wersji "szybkiej": IFFT, Inverse Fast Fourier Transform), które, jak łatwo się domyślić, wykonuje przekształcenie w odwrotnym kierunku. Bardzo często stosujemy ciąg operacji: FFT - przetwarzanie w dziedzinie częstotoliwości - IFFT i powrót do sygnału w dziedzinie czasu.
W module np.fft mamy funkcję np.fft.ifft dla ogólnego przypadku oraz np.fft.irfft jeżeli wynikiem ma być sygnał rzeczywisty. W tym drugim przypadku, aby otrzymać $N$ próbek sygnału, widmo podawane do funkcji powinno mieć długość $(N/2)+1$.
Spróbujmy wykorzystać IFFT do wygenerowania sygnału sinosuidalnego o częstotliwości 1000 Hz. Prążek widmowy odpowiadający częstotliwości $f$ ma w widmie indeks $N\cdot f/f_s$. Dla $N=2048$ i $f_s=48000$ Hz, dostajemy $n=42,666...$. Musimy tę wartość zaokrąglić. Biorąc prążek o indeksie $n=43$ otrzymujemy częstotliwość $f_s\cdot n/N$, czyli 1007.81 Hz. Aby otrzymać amplitudę sygnału równą 1, amplituda widmowa prążka musi być równa połowie długości transformaty, natomiast faza dla sygnału sinus musi wynosić $-\pi$. Tworzymy więc wektor zer o typie liczb zespolonych i wstawiamy zespoloną amplitudę prążka na właściwe miejsce.
widmo = np.zeros(1025, dtype=np.complex)
widmo[43] = -1024j
wynik = np.real(np.fft.irfft(widmo))
plt.plot(wynik[:240])
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał sinusoidalny wygenerowany metodą IFFT')
plt.show()
Okres sygnału o częstotliwości 1000 Hz wynosi (dla częstotliwości próbkowania 48000 Hz) 480 próbek, więc 240 próbek sygnału powinno objąć dokładnie 5 okresów sygnału. Ponieważ dokonaliśmy zaokrąglenia, okres uległ nieznacznemu skróceniu. Widać jednak że istotnie otrzymaliśmy w ten sposób sygnał sinusoidalny.
Estymacja częstotliwości prążków widmowych¶
Na zakończenie rozdziału pokażemy bardziej zaawansowany algorytm. Będzie on wyszukiwał częstotliwości prążków widma. Weźmy widmo fragmentu nagrania klarnetu użytego już wcześniej.
fragment_klarnet = klarnet[30000:32048]
fragment_klarnet = fragment_klarnet / np.max(np.abs(fragment_klarnet))
widmo_klarnet = 20 * np.log10(np.abs(np.fft.rfft(fragment_klarnet * np.hamming(2048))) / 1024)
f = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f, widmo_klarnet)
plt.xlim(0, 10000)
plt.ylim(-90, 0)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda widma [dB]')
plt.title('Widmo dźwięku klarnetu')
plt.show()
W pierwszym kroku musimy znaleźć początek i koniec każdego prążka. Możemy użyć progowania i szukać dla każdego prążka pierwszej i ostatniej wartości powyżej progu. Wartość progowa musi być tak dobrana, aby każda para sąsiednich prążków była rozdzielona. W powyższym przykładzie, na podstawie wykresu możemy wybrać wartość progową -45 dB.
plt.plot(f, widmo_klarnet)
plt.xlim(0, 10000)
plt.ylim(-45, 0)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda widma [dB]')
plt.title('Widmo progowane dźwięku klarnetu')
plt.show()
Algorytm działa następująco. Szukamy wartości widma powyżej progu:
ponad = (widmo >= prog)
Zmienna ponad będzie miała wartość True dla wartości powyżej lub równych progu oraz False dla pozostałych. Konwertujemy tę zmienną na liczby całkowite: True i False są przekształcane na odpowiednio 1 i 0.
ponad = (widmo >= prog).astype(np.int)
Następnie obliczamy pochodną, czyli po prostu różnicę bieżącej i poprzedniej wartości:
pochodna = np.diff(ponad)
Pochodna przyjmie wartości +1 na początku prążka (zmiana z 0 na 1) oraz -1 na końcu prążka (zmiana z 1 na 0), pozostałe wartości będą zerowe. Znajdujemy indeksy początku i końca każdego prążka za pomocą funkcji np.where. Funkcja zawsze zwraca listę, dlatego bierzemy jej pierwszy element [0]. Dodajemy 1 ponieważ pochodna pomija pierwszą wartość.
poczatki = np.where(pochodna == 1)[0] + 1
konce = np.where(pochodna == -1)[0] + 1
Teraz dla każdego prążka szukamy indeksu maksimum. Funkcja zip zwraca po jednym elemencie z każdego argumentu, u nas: pary (poczatek, koniec). W każdym zakresie [poczatek:koniec] szukamy maksimum funkcją np.argmax. Znaleziony indeks jest wyrażony w lokalnym zakresie, więc musimy dodać początek zakresu aby uzyskać indeks globalny. Na koniec, zamieniamy indeks na częstotliwość w Hz.
maksima = []
for poczatek, koniec in zip(poczatki, konce):
p = np.argmax(widmo[poczatek:koniec]) + poczatek
maksima.append(p * fs / rozmiar_okna)
Oto pierwsza wersja całej funkcji:
def maksima_widma(widmo, prog, rozmiar_okna, fs):
ponad = (widmo >= prog).astype(np.int)
pochodna = np.diff(ponad)
poczatki = np.where(pochodna == 1)[0] + 1
konce = np.where(pochodna == -1)[0] + 1
maksima = []
for poczatek, koniec in zip(poczatki, konce):
p = np.argmax(widmo[poczatek:koniec]) + poczatek
maksima.append(p * fs / rozmiar_okna)
return maksima
maksima = maksima_widma(widmo_klarnet, -45, 2048, wav_fs)
for m in maksima:
print('f = {:8.3f}, współczynnik = {:5.2f}'.format(m, m / maksima[0]))
Wykryta częstotliwość podstawowa (430,664 Hz) znacznie różni się od oczekiwanej (440 Hz), stosunki częstotliwości prążków harmonicznych do częstotliwości podstawowej również znacząco różnią się od liczb całkowitych. Dlaczego? Pamiętajmy o rozdzielczości FFT, u nas wynosi ona aż 21,5 Hz. Rzeczywiste maksima prążków mogą leżeć pomiędzy wartościami widma.
Mamy znalezione indeksy maksimów wśród wartości obliczonego widma. Możemy wziąć te indeksy, oraz wartość poprzednią i następną. Trzy punkty pozwolą nam dopasować parabolę, po czym możemy znaleźć jej maksimum. Jeżeli $k$ oznacza indeks maksimum prążka znaleziony poprzednio, wtedy:
$$a=\left|H(k-1)\right|$$ $$b=\left|H(k)\right|$$ $$c=\left|H(k+1)\right|$$
$$k'=k+\frac{a-c}{2(a-2b+c)}$$
gdzie $k'$ jest interpolowaną pozycją maksimum prążka, którą zamieniamy na częstotliwość tak jak poprzednio.
Źródło: https://www.dsprelated.com/freebooks/sasp/Quadratic_Interpolation_Spectral_Peaks.html.
def maksima_widma(widmo, prog, rozmiar_okna, fs):
ponad = (widmo >= prog).astype(np.int)
pochodna = np.diff(ponad)
poczatki = np.where(pochodna == 1)[0] + 1
konce = np.where(pochodna == -1)[0] + 1
maksima = []
for poczatek, koniec in zip(poczatki, konce):
p = np.argmax(widmo[poczatek:koniec]) + poczatek
a, b, c = widmo[p - 1:p + 2]
k = 0.5 * (a - c) / (a - 2 * b + c)
maksima.append((p + k) * fs / rozmiar_okna)
return maksima
maksima = maksima_widma(widmo_klarnet, -45, 2048, wav_fs)
for m in maksima:
print('f = {:8.3f}, współczynnik = {:5.2f}'.format(m, m / maksima[0]))
Widzimy, że otrzymaliśmy niemal idealne współczynniki harmonicznych. Zatem rzeczywista częstotliwość podstawowa dźwięku wynosi 441,335 Hz, jest więc nieznacznie odstrojona o 5,22 centów od założonej.
Podsumowanie¶
- Przekształcenie Fouriera przenosi sygnał z dziedziny czasu do dziedziny częstotliwości.
- Dla sygnałów cyfrowych stosujemy szybkie przekształcenie Fouriera FFT.
- Widmo (spektrum) to wynik obliczenia FFT dla bloku próbek sygnału. Widmo amplitudowe jest modułem widma, widmo fazowe - argumentem widma.
- Przecieki widma powstają gdy blok próbek poddawany FFT nie jest wielokrotnością okresu sygnału.
- Funkcje okna czasowego redukują przecieki widmowe. Różne okna (np. Hamminga, Blackmana) tłumią przecieki w różnym stopniu, ale większe tłumienie przecieków to zazwyczaj większa szerokość głównego prążka.
- Rozmiar okna czasowego decyduje o rozdzielczości FFT: dłuższe okno poprawia rozdzielczość częstotliwościową, ale pogarsza rozdzielczość czasową.
- Krótkookresowe przekształcenie Fouriera STFT stosujemy do analizy dłuższych sygnałów, po podzieleniu go na okna czasowe, które zwykle nakładają się na siebie (zakładkowanie).
- Spektrogram to wykres czas - częstotliwość - amplituda widma, pokazujący wynik STFT.
- Odwrotne przekształcenie Fouriera (szybka wersja: IFFT) zamienia widmo sygnału w postać czasową sygnału.
- Jeżeli chcemy znaleźć częstotliwości prążków widma, szukamy maksimów wśród wartości widma, a następnie stosujemy interpolację.
- Jeżeli nie wiemy jakie zastosować parametry analizy, to dla sygnałów dźwiękowych o częstotliwości próbkowania 48 lub 44,1 kHz, na początek warto spróbować: okno Hamminga, długość 2048 próbek, zakładkowanie STFT 75%.