%matplotlib inline
import numpy as np
import scipy.signal as sig
from scipy.io import wavfile
import matplotlib.pyplot as plt
from matplotlib import patches

from util import set_style
set_style()

# częstotliwość próbkowania
fs = 48000

np.set_printoptions(suppress=True)

Splot i korelacja ¶

Splot (convolution) jest jedną z najważnieszych operacji przetwarzania sygnałów. Zrozumienie zasady działania splotu sprawia jednak problemy początkującym w tej dziedzinie. Często pada również pytanie "do czego ten splot jest potrzebny". Spróbujemy odpowiedzieć na to pytanie.

Splot liniowy (linear convolution) w dziedzinie cyfrowej jest operacją wykonywaną na dwóch ciągach próbek, jego wynikiem jest nowa sekwencja próbek. Przedstawimy prosty przykład splotu dwóch sekwencji: $A = [1, 2, 3]$ oraz $B = [1, 2]$. Zasada obliczenia splotu jest następująca. Odwracamy jedną z sekwencji, np. B, w czasie ("tyłem do przodu") i wyrównujemy z drugą sekwencją (A) tak, że ostatnia próbka B po odwróceniu jest wyrównana z pierwszą próbką A. Mnożymy przez siebie te pokrywające się próbki i to jest pierwsza wartość wyniku splotu.

A      1  2  3          1 * 1 = 1
B   2  1                wynik = [1]

W kolejnym kroku przesuwamy odwrócony sygnał B w prawo o jedno miejsce. Teraz po dwie próbki każdego sygnału pokrywają się – mnożymy pokrywające się próbki przez siebie i dodajemy wyniki – to jest druga wartość wyniku.

A      1  2  3          1 * 2 + 2 * 1 = 4
B      2  1             wynik = [1, 4]

Przesuwamy ponownie odwrócony sygnał B, powtarzamy operacje. Kończymy obliczenia, jeżeli żadne próbki sygnałów nie pokrywają się.

A      1  2  3          2 * 2 + 3 * 1 = 7
B         2  1          wynik = [1, 4, 7]

A      1  2  3          3 * 2 = 6
B            2  1       wynik = [1, 4, 7, 6]

Operację splotu oznaczamy symbolem gwiazdki: $A \ast B = [1, 4, 7, 6]$. Splot jest operacją przemienną, mogliśmy odwrócić i przesuwać sygnał $A$, uzyskując taki sam wynik.

Splot możemy obliczyć za pomocą funkcji numpy.convolve.

print(np.convolve([1, 2, 3], [1, 2]))

[1 4 7 6]

Funkcja convolve może zwracać wynik na trzy sposoby, określane za pomocą trzeciego parametru, mode:

mode='full' (tryb domyślny) zwraca pełny wynik splotu,
mode='same' zwraca liczbę próbek równą długości dłuższego z sygnałów, wybranych ze środka wyniku,
mode='valid' zwraca tylko te próbki, dla których przynajmniej jeden z sygnałów jest w pełni pokryty przez drugi.

print('full: ', np.convolve([1, 2, 3], [1, 2], 'full'))
print('same: ', np.convolve([1, 2, 3], [1, 2], 'same'))
print('valid:', np.convolve([1, 2, 3], [1, 2], 'valid'))

full:  [1 4 7 6]
same:  [1 4 7]
valid: [4 7]

Moduł scipy.signal ma także własną funkcję scipy.signal.convolve, która działa identycznie do opisanej powyżej, a ponadto posiada dodatkowe możliwości, które omówimy później.

Jeżeli długości splatanych sygnałów wynoszą $L$ i $M$, wynik splotu liniowego ma długość $L+M-1$. Taka będzie też długość wyniku otrzymanego w trybie full.

W praktycznych zastosowaniach często jeden z sygnałów jest dłuższy od drugiego. Dłuższy sygnał jest często zmiennym sygnałem poddawanym przetwarzaniu, krótszy sygnał jest wtedy stały i jest on nazywany jądrem splotu (convolution kernel).

Systemy LTI, splot z odpowiedzią impulsową ¶

Bardzo ważną klasą systemów w przetwarzaniu sygnałów są systemy liniowe niezmiennicze w czasie (LTI, linear time-invariant). Liniowość jest rozumiana następująco. Jeżeli sygnał wejściowy $x_1(n)$ wywołuje na wyjściu systemu sygnał $y_1(n)$, a sygnał $x_2(n)$ wywołuje na wyjściu $y_2(n)$, to w systemie liniowym sygnał wejściowy $a_1 \cdot x_1(n) + a_2 \cdot x_2(n)$ wywołuje na wyjściu sygnał $a_1 \cdot y_1(n) + a_2 \cdot y_2(n)$ ($a_1$ i $a_2$ są liczbami rzeczywistymi). Warunek niezmienności w czasie jest następujący: jeżeli sygnał wejściowy $x(n)$ wywołuje na wyjściu systemu sygnał $y(n)$, to sygnał $x(n-M)$ wywołuje na wyjściu systemu sygnał $y(n-M)$, dla dowolnego $M$ (odpowiedź systemu na pobudzenie nie zależy od czasu, w którym to pobudzenie nastąpiło).

System LTI jest jednoznacznie opisany przez odpowiedź impulsową (impulse response). Jest to odpowiedź systemu na pobudzenie w formie impulsu jednostkowego (unit impulse), nazywanego też deltą Kroneckera, czyli sygnału mającego wartość jeden dla $n=0$ i zero dla pozostałych indeksów. Odpowiedź impulsową oznaczamy $h(n)$. Jeżeli znamy odpowiedź impulsową systemu LTI, to odpowiedź systemu na pobudzenie dowolnym sygnałem możemy obliczyć jako splot tego sygnału z odpowiedzią impulsową.

$$y(n) = h(n) \ast x(n)$$

Tyle teorii, a jak to się ma do praktyki? Weźmy pod uwagę filtr o skończonej odpowiedzi impulsowej (FIR). Jest to przykład systemu LTI. Znamy odpowiedź impulsową filtru – jest to po prostu zbiór współczynników filtru. Z powyższych rozważań wynika, że możemy obliczyć wynik filtracji wykonując splot sygnału ze zbiorem współczynników filtru. Sprawdźmy to.

b = sig.firwin(30, 1000, fs=fs)
szum = np.random.random(1000) * 2 - 1
# filtracja za pomocą funkcji lfilter
wynik_filtr = sig.lfilter(b, 1, szum)
# splot
wynik_splot = np.convolve(szum, b)[:1000]
print('Wyniki są równe:', np.allclose(wynik_filtr, wynik_splot))

Wyniki są równe: True

Funkcja numpy.allclose zwraca True jeżeli różnica wszystkich elementów podanych argumentów jest mniejsza od progu ($10^{-5}$). Mamy potwierdzenie: możemy użyć funkcji numpy.convolve do przeprowadzenia filtracji FIR metodą splotu z odpowiedzią impulsową filtru i uzyskujemy taki sam wynik jak z funkcji scipy.signal.lfilter.

Taka jest więc odpowiedź na pytanie "do czego jest potrzebny splot". Znając odpowiedź impulsową dowolnego systemu LTI możemy obliczyć wynik "przepuszczenia" sygnału przez ten system, po prostu obliczając splot liniowy tego sygnału ze znaną odpowiedzią impulsową systemu.

"Szybki" splot za pomocą FFT ¶

W poprzednich przykładach splataliśmy ze sobą dwa sygnały w dziedzinie czasu. W podobny sposób możemy wykonać splot dwóch transformat Fouriera sygnałów. Pokażemy teraz jak obliczyć splot w inny sposób. Należy zapamiętać następującą zasadę:

splotowi dwóch sygnałów w dziedzinie czasu odpowiada mnożenie ich transformat Fouriera,
mnożeniu dwóch sygnałów w dziedzinie czasu odpowiada splot ich transformat Fouriera.

Z powyższej zasady wynika, że możemy obliczyć wynik splotu sygnałów poprzez obliczenie ich transformat Fouriera, przemnożenie ich przez siebie i obliczenie odwrotnej transformaty. Nie możemy jednak po prostu wymnożyć transformat sygnałów, ponieważ w ten sposób obliczymy splot kołowy (circular convolution). Chcąc obliczyć splot liniowy, musimy wydłużyć oba sygnały do długości $L+M-1$, gdzie $L$ i $M$ to długości splatanych sygnałów, poprzez uzupełnienie ich zerami na końcu.

Wykonajmy tą metodą splot tych samych sygnałów co w pierwszym przykładzie.

widmo_a = np.fft.fft([1, 2, 3, 0])
widmo_b = np.fft.fft([1, 2, 0, 0])
widmo_splot = widmo_a * widmo_b
splot = np.real(np.fft.ifft(widmo_splot))
print(splot)

[ 1.  4.  7.  6.]

Widzimy, że otrzymaliśmy prawidłowy wynik.

Operację obliczania splotu w ten sposób często nazywa się szybkim splotem (fast convolution). W powyższym przykładzie oba sygnały były bardzo krótkie, więc splot w dziedzinie czasu jest liczony szybciej. W praktycznych zastosowaniach, splatane sygnały są zwykle dużo dłuższe, a wtedy liczba operacji mnożenia i dodawania potrzebna do obliczenia splotu w dziedzinie czasu gwałtownie wzrasta. W takich przypadkach obliczenie splotu za pomocą transformat FFT sygnałów jest zazwyczaj szybsze. Wzrost szybkości jest szczególnie zauważalny na procesorach sygnałowych, które są w stanie obliczać FFT bardzo szybko, zwłaszcza jeżeli długość transformaty jest potęgą dwójki. Ponadto często zdarza się, że splot jest obliczany wielokrotnie, przy czym jeden sygnał się zmienia, a drugi (jądro splotu) pozostaje stały. W takim przypadku wystarczy obliczyć transformatę jądra splotu tylko raz, co jeszcze bardziej przyspiesza obliczenia.

W systemach LTI, transformatę Fouriera odpowiedzi impulsowej systemu nazywa się transmitancją systemu (transfer function), oznaczaną $H(f)$. Transformata Fouriera sygnału na wyjściu układu jest więc równa transformacie sygnału wejściowego przemnożonej przez transmitancję systemu.

$$Y(f) = H(f) \cdot X(f)$$

Funkcja scipy.signal.fftconvolve realizuje splot metodą mnożenia transformat. Znaczenie trzeciego parametru (mode) jest takie samo jak dla zwykłego splotu (full, valid, same). Funkcja ta sama uzupełnia sygnały zerami.

print(sig.fftconvolve([1, 2, 3], [1, 2]))

[ 1.  4.  7.  6.]

Wspomnieliśmy wcześniej, że funkcja scipy.signal.convolve ma dodatkową właściwość w porównaniu z np.convolve. Otóż ma ona parametr method, który określa sposób obliczania splotu. Dla method='direct' obliczany jest splot w dziedzinie czasu, dla method='fft' – za pomocą FFT. Natomiast domyślną wartością parametru method jest 'auto', w tym przypadku funkcja sama decyduje o tym która metoda będzie szybsza. Możemy też sami to sprawdzić za pomocą funkcji scipy.signal.choose_conv_method. Na przykład, wybór metody splotu z jądrem o długości 100, przy różnej długości sygnału:

print('N=1000: method =', sig.choose_conv_method(np.arange(1000), np.arange(100)))
print('N=2000: method =', sig.choose_conv_method(np.arange(2000), np.arange(100)))

N=1000: method = direct
N=2000: method = fft

Splot blokowy, overlapp-add ¶

W praktyce czesto występuje sytuacja, w której sygnał poddawany splotowi nie jest skończony. Próbki sygnału cały czas napływają na wejście systemu i muszą być splatane ze stałym jądrem splotu. W takim przypadku potrzebne jest zastosowanie splotu blokowego: sygnał wejściowy jest dzielony na bloki, każdy z bloków jest splatany z jądrem. Nie można jednak traktować każdego bloku próbek niezależnie, należy wziąć pod uwagę próbki, które znajdowały się w poprzednim bloku (niezerowy warunek początkowy). Ponadto, jeżeli blok próbek ma długość $L$, a jądro ma długość $M$, w wyniku splotu otrzymujemy $L+M-1$ próbek, podczas gdy oczekujemy takiej samej liczby próbek, jaka była na wejściu ($L$).

Metoda overlap-add (OLA) rozwiązuje te problemy poprzez zakładkowanie wyników cząstkowych splotów: $M-1$ ostatnich próbek wyniku splotu dla poprzedniego bloku jest dodawanych do pierwszych próbek wyniku splotu dla bieżącego bloku. Algorytm można opisać następująco. Stałe jądro splotu ma długość $M$, uzupełniamy je zerami do długości $L+M-1$. Tworzymy bufor o długości $M-1$ i wypełniamy go zerami. Następnie dla każdego bloku próbek sygnału wejściowego o długości $L$:

uzupełniamy blok próbek $M-1$ zerami,
obliczamy splot bloku próbek z jądrem, dowolną metodą (z definicji lub przez FFT),
do pierwszych $M-1$ próbek dodajemy zawartość bufora,
początkowe $L$ próbek wyniku jest branych jako wynik splotu bieżącego bloku,
pozostałe $M-1$ próbek jest zapisywanych w buforze.

Poniżej symulujemy blokowy splot sygnału przypadkowego z odpowiedzią impulsową filtru o długości 31. Funkcja scipy.signal.convolve samodzielnie uzupełnia sygnały zerami.

L = 100
M = 31
b = sig.firwin(M, 1000, fs=fs)
szum = np.random.random(10 * L) * 2 - 1
bufor = np.zeros(M - 1)
splot_bloki = []

for i in range(0, len(szum) - L + 1, L):
    blok = szum[i:i+L]
    splot = sig.convolve(blok, b)
    splot[:M-1] += bufor
    wynik_bloku = splot[:L]
    # tutaj mamy wynik splotu dla bloku, który możemy dalej przetwarzać
    splot_bloki.append(wynik_bloku)
    bufor[:] = splot[L:]

wynik = np.hstack(splot_bloki)
print('Jednakowe wyniki:', np.allclose(wynik, sig.convolve(szum, b)[:1000]))

Jednakowe wyniki: True

Wynik splotu metodą overlapp-add jest identyczny jak przy obliczeniu splotu całego sygnału na raz. Metoda overlapp-add jest często wykorzystywana w wielu algorytmach opartych na blokowym splataniu sygnału.

Istnieje również metoda overlapp-save (OLS), która polega na zakładkowaniu sygnału przed wykonaniem splotu i na odrzuceniu nadmiarowych próbek wyniku. Złożoność obliczeniowa obu metod jest prawie jednakowa.

Korelacja wzajemna ¶

Korelacja w statystyce określa współzależność zmiennych losowych. W przetwarzaniu sygnałów, korelacja wzajemna (cross-correlation), nazywana też korelacją skrośną lub krzyżową, jest miarą podobieństwa dwóch sygnałów. Korelację oblicza się jako sumę wartości obu sygnałów przemnożonych przez siebie (dot product). Obliczenia korelacji są wykonywane dla różnych wartości przesunięcia czasowego (opóźnienia, lag). Jeden z sygnałów jest przesuwany względem drugiego, dla każdego przesunięcia mnożymy przez siebie nakładające się próbki i sumujemy wyniki. Widać więc, że korelacja wzajemna jest bardzo podobna do splotu. Różnica jest taka, że przy obliczaniu splotu odwracamy jeden z sygnałów w czasie, a przy obliczaniu korelacji nie robimy tego.

W praktyce wykorzystujemy korelację wzajemną np. do wykrywania opóźnień sygnałów. Jeżeli obliczymy wartości korelacji dla różnych przesunięć, maksymalna wartość korelacji wskaże nam opóźnienie sygnałów względem siebie.

Pokażemy obliczenia korelacji wzajemnej na prostym przykładzie sygnałów: $A=[0, 0, 1, 2, 3, 2, 1]$, $B=[1, 2, 3, 2, 1]$. Widzimy, że sygnał $A$ jest opóźniony względem $B$ o dwie próbki. Do obliczenia korelacji wzajemnej możemy użyć funkcji numpy.correlate lub scipy.signal.correlate, ich działanie jest analogiczne do funkcji convolve. Jest jedna różnica: dla numpy.correlate domyślnym trybem jest 'valid', podczas gdy dla scipy.signal.correlate – 'full', co jest bardziej intuicyjne.

A = [0, 0, 1, 2, 3, 2, 1]
B = [1, 2, 3, 2, 1]
C = sig.correlate(A, B)
print(C)

[ 0  0  1  4 10 16 19 16 10  4  1]

Jak zinterpretować ten wynik? Pierwsza wartość wyniku (zerowy indeks) jest obliczona dla przypadku, gdy ostatnia próbka $B$ pokrywa się z pierwszą próbką $A$, zatem przesunięcie dla tego przypadku wynosi $M-1$, gdzie $M$ jest długością sygnału $B$. Wartość przesunięcia dające największą wartość korelacji możemy obliczyć następująco:

print('Przesunięcie =', np.argmax(C) - (len(B) - 1))

Przesunięcie = 2

plt.plot(np.arange(len(C)) - (len(B) - 1), C)
plt.axvline(0, lw=1, c='k')
plt.xlabel('opóźnienie')
plt.ylabel('korelacja')
plt.title('Korelacja wzajemna sygnałów A i B')
plt.show()

Otrzymaliśmy więc spodziewany wynik.

Zauważmy, że możemy otrzymać dokładnie ten sam wynik obliczając splot $A$ z sygnałem $B$ odwróconym w czasie.

print(sig.convolve(A, B[::-1]))

[ 0  0  1  4 10 16 19 16 10  4  1]

Funkcje correlate i convolve robią więc niemal to samo.

Pokażemy teraz bardziej praktyczny przykład. Sygnał jest przetwarzany przez filtr FIR. Przetworzony sygnał jest dodatkowo opóźniony względem oryginalnego. Widząc projekt filtru, wiemy ile to opóźnienie wynosi. Załóżmy jednak, że nie znamy parametrów filtru. Mając do dyspozycji jedynie oba sygnały, mamy je wyrównać w czasie.

h = sig.firwin(31, 1000, fs=fs)
x = sig.chirp(np.linspace(0, 1, 2001), 20, 1, 1000)
y = sig.lfilter(h, 1, x)
plt.plot(x[:200], label='oryginalny')
plt.plot(y[:200], label='po filtracji')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnały przed i po filtracji')
plt.legend()
plt.show()

c = sig.correlate(y[:200], x[:200])
plt.figure()
plt.plot(np.arange(len(c)) - 199, c)
plt.xlabel('przesunięcie')
plt.ylabel('korelacja wzajemna')
plt.title('Korelacja sygnałów przed i po filtracji')
plt.show()
print('Przesunięcie =', np.argmax(c) - 199)

Przesunięcie = 15

Używamy do obliczeń pierwszych 200 próbek, zakładając że opóźnienie mieści się w tym zakresie. Największą korelację znaleziono dla przesunięcia równego 15 próbek. Jest to prawidłowy wynik, filtr o 31 współczynnikach ma opóźnienie grupowe równe dokładnie 15 próbkom. Korelacja wzajemna pozwoliła więc nam na znalezienie szukanej wartości opóźnienia. Możemy teraz wyrównać sygnały w czasie:

plt.plot(x[:200], label='oryginalny')
plt.plot(y[15:215], label='po filtracji')
plt.xlabel('nr próbki')
plt.ylabel('amplituda [dB]')
plt.title('Sygnały z kompensacją opóźnienia')
plt.legend()
plt.show()

Obliczanie korelacji wzajemnej za pomocą FFT ¶

Pokazaliśmy wcześniej w jaki sposób można obliczyć splot za pomocą mnożenia transformat Fouriera sygnałów. Skoro korelacja wzajemna jest podobna do splotu, można również wykorzystać FFT do obliczenia korelacji. Najprostszą metodą jest odwrócenie w czasie jedego z sygnałów przed obliczeniem FFT. Można jednak postąpić inaczej. Dla sygnałów rzeczywistych, odwróceniu sygnału w czasie odpowiada obliczenie zespolonego sprzężenia (complex conjugate) jego transformaty. Sprzężenie liczby zespolonej obliczamy przez odwrócenie znaku części urojonej. Możemy je obliczyć za pomocą funkcji numpy.conj. Tak jak w przypadku splotu, musimy pamiętać o uzupełnieniu sygnałów zerami. Możemy użyć do tego funkcji numpy.pad. Obliczmy ten sam przykład co wcześniej dla korelacji wzajemnej z definicji.

widmo_A = np.fft.fft(np.pad(A, (0, len(B) - 1), 'constant'))
widmo_B = np.fft.fft(np.pad(B, (0, len(A) - 1), 'constant'))
widmo_kor = widmo_A * np.conj(widmo_B)
kor = np.real(np.fft.ifft(widmo_kor))
print(kor)

[ 10.  16.  19.  16.  10.   4.   1.  -0.   0.   1.   4.]

Wynik jest nieco inny niż dla obliczenia korelacji wzajemnej z definicji: jest on przesunięty cyklicznie. Pierwsza wartość wyniku jest obliczona dla zerowego przesunięcia, ostatnie cztery wartości są obliczone dla ujemnych przesunięć. Możemy uzyskać ten sam wynik co poprzednio, dokonując cyklicznego przesunięcia wyniku za pomocą funkcji numpy.roll:

print(np.roll(kor, 4))

[ -0.   0.   1.   4.  10.  16.  19.  16.  10.   4.   1.]

Aby obliczyć korelację wzajemną metodą FFT, możemy skorzystać z funkcji scipy.signal.correlate z parameterem method='fft'.

Autokorelacja ¶

Autokorelacja (autocorrelation) jest szczególnym przypadkiem korelacji wzajemnej sygnału "samego ze sobą". Autokorelacja dla danej wartości przesunięcia opisuje podobieństwo sygnału do jego przesuniętej kopii. Największe podobieństwo, co oczywiste, uzyskamy dla zerowego przesunięcia, gdy oba korelowane sygnały są identyczne. Autokorelacja sygnału dla zerowego przesunięcia jest równa energii sygnału, czyli sumie kwadratów próbek. Unormowana funkcja autokorelacji powstaje przez podzielenie wartości funkcji autokorelacji przez energię sygnału, tak że dla zerowego przesunięcia, wartość autokorelacji przyjmuje wartość maksymalną, równą 1.

Podstawowym zastosowaniem autokorelacji jest wyszukiwanie okresowości w sygnale. Jeżeli sygnał jest idealnie okresowy (np. sygnał sinusoidalny), maksima funkcji autokorelacji będą znajdowały się w punktach równych wielokrotności okresu. Dla sygnałów "pseudookresowych", w których kolejne okresy mogą różnić się od siebie (np. dźwięków instrumentów muzycznych), autokorelacja pozwoli znaleźć długość powtarzającego się fragmentu sygnału.

Obliczmy funkcję autokorelacji dla sygnału trójkątnego, wygenerowanego za pomocą funkcji scipy.signal.sawtooth. Sekwencja ta powtarza się z okresem o długości 20 próbek.

trojkat = sig.sawtooth(2 * np.pi * 5 / 100 * np.arange(100), 0.5)
plt.plot(trojkat)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał trójkątny')
plt.show()

Moduły numpy i scipy nie mają specjalnych funkcji do obliczenia autokorelacji. Musimy skorzystać z funkcji do obliczania korelacji wzajemnej. W przypadku stosowania funkcji np.correlate musimy pamiętać o ustawieniu trybu na full.

autokor = sig.correlate(trojkat, trojkat)
plt.plot(autokor)
plt.xlabel('indeks wartości')
plt.ylabel('autokorelacja')
plt.title('Funkcja autokorelacji sygnału trójkątnego')
plt.show()

Funkcja autokorelacji obliczona jako korelacja wzajemna jest zawsze symetryczna (nie ma znaczenia w którą stronę przesuwamy kopię sygnału). Potrzebujemy części, która zaczyna się od zerowego przesunięcia. Musimy więc odrzucić pierwszych $L-1$ wartości, gdzie $L$ oznacza długość sygnału. Ponadto możemy unormować funkcję autokorelacji w taki sposób, aby dla zerowego przesunięcia wartość autokorelacji wynosiła jeden.

autokor = sig.correlate(trojkat, trojkat)
autokor = autokor[len(autokor)//2:]
autokor = autokor / autokor[0]
plt.plot(autokor)
plt.xlabel('przesunięcie')
plt.ylabel('autokorelacja')
plt.title('Unormowana funkcja autokorelacji sygnału trójkątnego')
plt.show()

Powyższy wykres można zinterpretować w następujący sposób. Dla zerowego przesunięcia mamy maksymalną wartość autokorelacji. W miarę zwiększania się przesunięcia, kopie sygnału stają się coraz mniej podobne do siebie. W chwili gdy minimum jednego sygnału pokrywa się z maksimum drugiego (przesunięcie równe połowie okresu), autokorelacja przyjmuje wartość minimalną. Gdy dalej zwiększamy przesunięcie, okresy zaczynają nakładać się na siebie i wartość autokorelacji wzrasta. Gdy przesunięcie jest równe długości okresu, wartość autokorelacji osiąga lokalne maksimum. Dalej cykl się powtarza.

Na wykresie widać, że dla przesunięć równych kolejnym wielokrotnościom okresu, wartości autokorelacji są coraz mniejsze, chociaż przesunięcie sygnału o dowolną liczbę okresów daje taki sam sygnał. Wynika to z liniowego sposobu liczenia autokorelacji – dla większych przesunięć, mniej probek bierze udział w obliczeniach, stąd mniejsze wartości. Gdybyśmy policzyli funkcję autokorelacji w sposób kołowy, wszystkie maksima miałyby taką samą amplitudę. W analizie sygnałów nie ma to zwykle znaczenia, ponieważ w większości przypadków zależy nam na znalezieniu pierwszego maksimum.

Obliczanie autokorelacji za pomocą FFT ¶

Wiemy już, że korelację wzajemną możemy obliczyć za pomocą FFT. W ten sam sposób można obliczyć też autokorelację. W przypaku sygnałów rzeczywistych, mnożeniu transformaty sygnału przez jej zespolone sprzężenie odpowiada kwadrat modułu transformaty. W taki sposób jest zazwyczaj liczona autokorelacja metodą FFT.

Poniżej obliczamy autokorelację sygnału trójkątnego podaną wyżej metodą. Wynik normujemy dzieląc przez wynik otrzymany dla zerowego przesunięcia. Z obliczonego wyniku bierzemy pierwszą połowę wartości autokorelacji, druga połowa jest symetryczna i odpowiada ujemnym przesunięciom.

autokor_fft = np.real(np.fft.ifft(np.abs(np.fft.fft(trojkat))**2))
autokor_fft /= autokor_fft[0]
plt.plot(autokor_fft[:len(autokor_fft) // 2])
plt.xlabel('opóźnienie')
plt.ylabel('autokorelacja')
plt.title('F. autokorelacji s. trójkątnego obliczona przez FFT')
plt.show()

W ten sposób otrzymujemy autokorelację kołową. Amplituda maksimów dla przesunięć równych wielokrotnościom okresu jest jednakowa. Aby otrzymać autokorelację liniową, musimy uzupełnić sygnał $L-1$ zerami, gdzie $L$ jest długością analizowanego sygnału.

trojkat2 = np.pad(trojkat, (0, len(trojkat) - 1), 'constant')
autokor_fft2 = np.real(np.fft.ifft(np.abs(np.fft.fft(trojkat2))**2))
autokor_fft2 /= autokor_fft2[0]
plt.plot(autokor_fft2[:len(autokor_fft2) // 2])
plt.xlabel('opóźnienie')
plt.ylabel('autokorelacja')
plt.title('F. autokorelacji s. trójkątnego obliczona przez FFT')
plt.show()

Otrzymjemy identyczny wynik jak dla autokorelacji obliczonej wcześniej z definicji.

Aby obliczyć autokorelację metodą FFT, możemy skorzystać z funkcji scipy.signal.correlate z parameterem method='fft'.

Widma funkcji korelacji i autokorelacji ¶

Postacie widmowe funkcji korelacji wzajemnej i autokorelacji mają również interpretację fizyczną. Transformatę Foriera funkcji autokorelacji sygnału $x(n)$ oblicza się jako:

$$S_{xx}(k) = X(k) \cdot X(k)^* = \left|X(k)\right|^2$$

Funkcja $S_{xx}$ reprezentuje widmową gęstość mocy (PSD, Power spectra density), nazywaną też periodogramem, chociaż perodogram zwykle oblicza się inna metodą. Widmowa gęstość mocy opisuje rozkład częstotliwościowy mocy sygnału.

Tranformata Fouriera funkcji korelacji wzajemnej sygnałów $x(n)$ i $u(n)$ jest obliczana jako:

$$S_{xy}(n) = X(k)^* \cdot Y(k)$$

Funkcję $S_{xy}$ nazywa się widmem wzajemnym (cross-spectrum) lub wzajemną gęstością widmową (CSD, cross-spectral density). Funkcja ta opisuje wzajemną zależność wartości dwóch sygnałów w dziedzinie częstotliwości.

Koherencję (coherence) reprezentuje wzajemną zależność mocy dwóch sygnałów w dziedzinie częstotliwości. Obliczamy ją następująco:

$$C_{xy} = \frac{\left|S_{xy}\right|^2} {S_{xx} S_{yy}}$$

Wymienione wyżej funkcje są stosowane przede wszystkim do analizy sygnałów o charakterze losowym (stochastycznych), również zaszumionych sygnałów. Ze względu na naturę takich sygnałów, nie wystarczy obliczyć tych funkcji dla jednego bloku próbek. W praktyce stosuje się podział sygnału na odcinki za pomocą okien czasowych (z zakładkowaniem), obliczenie pożądanych funkcji oraz uśrednienie i skalowanie wyników z poszczególnych okien.

W module scipy.signal mamy do dyspozycji następujące funkcje:

welch – oblicza widmową gęstość mocy (PSD),
periodogram – oblicza periodogram, wywołuje funkcję welch,
csd – oblicza wzajemną gęstość widmową (CSD)
coherence – oblicza koherencję.

Znajdowanie częstotliwości podstawowej za pomocą autokorelacji ¶

W rozdziale poświęconym analizie widmowej sygnałów znajdowaliśmy częstotliwości prążków widmowych w nagraniu instrumentu muzycznego. Częstotliwość pierwszego prążka harmonicznego wyznacza częstotliwość podstawową sygnału, odpowiadającą za wysokość dźwięku na skai muzycznej. Metoda wyszukiwania częstotliwości prążka w FFT jest dla wielu osób "oczywistym" podejściem do wyznaczania częstotliwości podstawowej. Jest to jednak metoda niedokładna z powodu małej rozdzielczości FFT, co wymaga stosowania interpolacji wartości widma. Obliczymy teraz częstotliwość podstawową sygnału z użyciem autokorelacji.

Wybieramy fragment sygnału ze środka nagrania, tam gdzie dźwięk jest najbardziej stabilny (stan ustalony).

wav_fs, klarnet_wav = wavfile.read('pliki/klarnet.wav')
klarnet = klarnet_wav[30000:31024]
klarnet = klarnet / np.max(np.abs(klarnet))
plt.plot(np.arange(30000, 31024), klarnet)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Fragment nagrania dźwięku klarnetu')
plt.xlim(30000, 30500)
plt.show()

Widać, że pewien fragment sygnału powtarza się cyklicznie. Musimy znaleźć długość tego okresu. Obliczamy funkcję autokorelacji. Przed jej obliczeniem powinniśmy usunąć z sygnału składową stałą, poprzez odjęcie od sygnału jego wartości średniej. Operację tę nazywa się usuwaniem trendu (detrending). Jeżeli jej nie wykonamy, obliczona autokorelacja będzie znieształcona, a ponadto w omawianym przypadku dostaniemy błąd dzielenia przez zero przy normalizacji.

klarnet = klarnet - np.mean(klarnet)
klarnet_autokor = sig.correlate(klarnet, klarnet)
klarnet_autokor = klarnet_autokor[len(klarnet_autokor)//2:]
klarnet_autokor = klarnet_autokor / klarnet_autokor[0]
plt.plot(klarnet_autokor)
plt.xlabel('opóźnienie')
plt.ylabel('autokorelacja')
plt.title('F. autokorelacji dźwięku klarnetu')
plt.show()

Musimy teraz wyznaczyć indeks pierwszego "pełnego" maksimum funkcji. Ponieważ obliczyliśmy autokorelację liniową, jest to najwyższe spośród pełnych maksimów, ale nie możemy użyć funkcji max z powodu maksimum w zerze. Aby znaleźć szukany indeks, musimy wykonać kilka sztuczek z funkcją autokorelacji.

Najpierw kilka ogólnych uwag o zakresach w języku Python. Zakresy definiujemy w formie [start:stop:krok], gdzie start to indeks pierwszej wartości, stop to indeks za ostatnią wartością zakresu (przedział półotwarty), krok określa o ile zwiększany jest indeks. Pominięcie którejkolwiek wartości oznacza użycie domyślnych wartości: zero dla start, długość sekwencji dla stop i jeden dla krok.

Algorytm działa następująco.

Zerujemy wartości mniejsze od zera.
Bierzemy pierwszą połowę wartości i "rozciagamy" ją na pełny zakres.
Odejmujemy rozciągniętą kopię od oryginalnej funkcji.
Ponownie zerujemy ujemne wartości.

Najpierw tworzymy bufor o rozmiarze równym długości obliczonej funkcji autokorelacji. Niech $n$ oznacza połowę tej długości (po zaokrągleniu w dół).

Wybieramy pierwszą połowę wartości za pomocą zakresu [:n].
Wpisujemy ją do bufora pod co drugi indeks (zaczynając od zera: 0, 2, 4, ...): [::2].
Pozostałe miejsca bufora wypełnimy wartościami średnimi:
- bierzemy znów pierwszą połowę wartości (0, 1, 2, ...): [:n],
- bierzemy wartości kolejne po wybranych wyżej (1, 2, 3, ...): [1:n+1],
- obliczamy ich średnią,
- wstawiamy w puste miejsca bufora, co drugi indeks zaczynając od 1 (1, 3, 5, ...): [1::2].

autokor2 = np.maximum(klarnet_autokor, 0)
bufor = np.zeros(autokor2.shape)
n = len(autokor2) // 2
bufor[::2] = autokor2[:n]
bufor[1::2] = (autokor2[:n] + autokor2[1:n+1]) / 2
autokor2 = np.maximum(autokor2 - bufor, 0)
plt.plot(autokor2)
plt.xlabel('opóźnienie')
plt.ylabel('autokorelacja')
plt.title('Zmodyfikowana f. autokorelacji dźwięku klarnetu')
plt.show()

Powyższe operacje spowodowały odjęcie od funkcji autokorelacji jej "zerowego" maksimum i wszystkich parzystych maksimów.

Znajdujemy indeks $m$, który wyznacza położenie maksimum tej funkcji. Okres sygnału jest równy $m \cdot T$, gdzie $T$ jest okresem próbkowania, $T = 1 / f_s$. Szukana częstotliwość jest równa $f = 1 / (mT) = f_s / m$.

m = np.argmax(autokor2)
print('Częstotliwość podstawowa = {:.2f} Hz'.format(wav_fs / m))

Częstotliwość podstawowa = 441.00 Hz

Otrzymany wynik wskazuje na wysokość dźwięku $a^1$, której odpowiada dokładna częstotliwość 440 Hz. Przy użyciu metody FFT uzyskaliśmy wynik 430,664 Hz, który znacznie różni się od spodziewanego z powodu małej rozdzielczości FFT. Stosując interpolację widma otrzymaliśmy wynik 441.335 Hz. Stosując autokorelację uzyskujemy wynik dokładniejszy niż dla FFT, bez konieczności używania interpolacji, która daje tylko przybliżenie rzeczywistych wartości. W przypadku autokorelacji możemy również zastosować interpolację kwadratową, w identyczny sposób jak dla FFT.

a, b, c = autokor2[m-1:m+2]
k = 0.5 * (a - c) / (a - 2 * b + c)
print('Częstotliwość podstawowa = {:.3f} Hz'.format(wav_fs / (m + k)))

Częstotliwość podstawowa = 441.682 Hz

Rozplot, filtracja odwrotna ¶

Mamy następującą sytuację. Użyteczny sygnał $x(n)$ jest zniekształcany przez system LTI o znanej nam odpowiedzi impulsowej $h(n)$. Na wyjściu systemu obserwujemy sygnał $y(n)$, który jest splotem sygnału $x(n)$ z odpowiedzią impulsową $h(n)$. Chcemy odzyskać oryginalny sygnał $x(n)$. W tym celu musimy wykonać operację odwrotną do splotu. Formalnie taka operacja nie jest zdefiniowana. Jednak wiemy, że splot może być obliczony przez przemnożenie transformat: $X(f) \cdot H(f)$. Odwrotnością mnożenia jest dzielenie, zatem możemy obliczyć transformatę sygnału wejściowego $X(f)$ mnożąc $Y(f)$ przez odwrotność transmitancji $H(f)$. W dziedzinie czasu operacji tej odpowiada splot $y(n)$ z odpowiedzią impulsową układu o odwróconej transmitancji. Taka operacja (przeprowadzana w dziedzinie częstotliwości lub czasu) nazywana jest rozplotem (deconvolution). Odpowiada jej operacja filtracji za pomocą filtru, którego transmitancja jest odwrotnością transmitancji układu zniekształcającego – operację tą nazywa się więc filtracją odwrotną (inverse filtering).

Aby była możliwa filtracja odwrotna, transmitancja systemu musi być odwracalna. Jeżeli odwracamy transmitancję, jej zera stają się biegunami i na odwrót. Wiemy, że warunkiem stabilności układu jest to, że moduły biegunów muszą być mniejsze od jedności. Aby filtr odwrotny był stabilny, również zera oryginalnej transmitancji muszą mieć moduł mniejszy od jedności. Taki system nazywa się minimalnofazowym. Liniowofazowe filtry FIR nie są, niestety, minimalnofazowe, filtr odwrotny jest albo na granicy stabilności, albo niestabilny.

Trzeba tu wyraźnie zaznaczyć, że w większości praktycznych przypadków nie jest możliwe odzyskanie oryginalnego sygnału. Jeżeli zastosujemy filtr tłumiący wybrane zakresy częstotliwości, filtracja odwrotna nie odzyska wytłumionych składowych, jedynie spowoduje wzmocnienie szumu.

Echo jest problemem w wielu systemach telekomunikacyjnych. W poniższym przykładzie pokażemy trywialny problem usuwania pojedynczego echa o stałym opóźnieniu i stałej amplitudzie. Sygnał złożony z trzech składowych sinusoidalnych przechodzi przez system, który dodaje echo o stałym opóźnieniu równym 100 próbek i amplitudzie równej połowie amplitudy oryginalnego sygnału. Ten proces możemy łatwo opisać za pomocą filtru FIR.

n = np.arange(1000)
x = np.sin(2 * np.pi * n * 200 / fs) + np.sin(2 * np.pi * n * 400 / fs) + np.sin(2 * np.pi * n * 600 / fs)
h = np.zeros(101)
h[0] = 1
h[100] = 0.5
y = sig.lfilter(h, 1, x)
plt.plot(x, label='oryginalny')
plt.plot(y, label='z echem')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał oryginalny i z dodanym echem')
plt.legend(loc='upper right')
plt.show()

Widzimy zniekształcenie sygnału przez echo. Spróbujemy dokonać rozplotu. Zamiast obliczać rozplot za pomocą FFT, możemy użyć ponownie funkcji lfilter dla odwróconej transmitancji (zamieniamy b z a).

x_est = sig.lfilter([1], h, y)
plt.plot(x, label='oryginalny')
plt.plot(x_est, label='odzyskany')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Wynik rozplotu sygnału')
plt.legend(loc='upper right')
plt.show()

W rozpatrywanym prostym przypadku otrzymaliśmy oryginalny sygnał. Praktyczne przypadki usuwania echa są dużo bardziej złożone i wymagają skomplikowanych podejść, np. opartych na cepstrum sygnału.

W module scipy jest dostępna funkcja scipy.signal.deconvolve. Działa ona inaczej niż opisana wyżej procedura. Tworzy ona filtr IIR, w którym współczynnikami licznika są próbki sygnału $y(n)$, natomiast współczynniki mianownika są tworzone przez odpowiedź impulsową systemu, który przetworzył sygnał. Odpowiedź impulsowa takiego filtru jest odtworzonym sygnałem. Funkcja zwraca dwa wyniki: tablicę próbek odtworzonego sygnału oraz "resztę" powstałą po obliczeniu różnicy między $y(n)$ a wynikiem splotu odtworzonego sygnału z odpowiedzią impulsową. Funkcja nie zwraca pełnego sygnału, liczba próbek wyniku jest o $M-1$ mniejsza od długości $y(n)$, gdzie $M$ jest długością odpowiedzi impulsowej. Widać, że również w tym przypadku uzyskujemy poprawny wynik rozplotu.

x_est2, res = sig.deconvolve(y, h)
plt.plot(x, label='oryginalny')
plt.plot(x_est2, label='odzyskany')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Wynik rozplotu funkcją deconvolve')
plt.legend(loc='upper right')
plt.show()
print('Średnia energia reszty:', np.mean(res**2))

Średnia energia reszty: 0.323411124889

Estymacja odpowiedzi impulsowej. Sygnał MLS ¶

W innej, również często spotykanej w praktyce sytuacji, mamy do czynienia z systemem LTI, który jest dla nas "czarną skrzynką" – nie znamy jego parametrów. Możemy podać na jego wejście dowolny sygnał i zarejestrować sygnał wyjściowy. Na tej podstawie chcemy otrzymać odpowiedź impulsową systemu. Pozornie sprawa jest prosta, wystarczy podać na wejście impuls jednostkowy. Problem polega na tym, że w praktyce bardzo trudno jest wytworzyć sygnał odpowiadający idealnemu impulsowi jednostkowemu. Na przykład, chcemy uzyskać odpowiedź impulsową pomieszczenia, np. sali koncertowej, traktując ją jako system LTI (co jest uproszczeniem). Tradycyjna metoda polega na strzelaniu z pistoletu hukowego i zarejestrowaniu dźwięku. Wystrzał z pistoletu nie jest jednak dokładnie impulsem jednostkowym, a ponadto zarejestrowany sygnał jest bardzo krótki, a przez to wynik jest obarczony dużym szumem.

Jeżeli transformaty sygnału wejściowego i wyjściowego oznaczymy przez $X(f)$ i $Y(f)$, a szukaną transmitancję przez $H(f)$, wtedy znana już zależność:

$$Y(f) = H(f) \cdot X(f)$$

może być przekształcona przez przemnożenie obu stron przez zespolone sprzężenie $X(f)$:

$$Y(f) \cdot X^*(f) = H(f) \cdot X(f) \cdot X^*(f)$$

Wykorzystując wcześniej podane zależności na obliczanie widma korelacji wzajemnej $S_{xy}(f)$ i autokorelacji $S_{xx}(f)$, możemy zapisać tę zależność w formie:

$$S_{xy}(f) = H(f) \cdot S_{xx}(f)$$

albo jako splot w dziedzinie czasu:

$$R_{xy}(n) = h(n) \ast R_{xx}(n)$$

Jeżeli na wejście badanego systemu podamy sygnał $x(n)$, którego autokorelacja jest równa impulsowi jednostkowemu, to splot odpowiedzi impulsowej z impulsem jednostkowym będzie równy tej samej odpowiedzi impulsowej. A więc obliczając korelację wzajemną między sygnałem wyjściowym a wejściowym, otrzymamy szukaną odpowiedź impulsową.

Sekwencja o maksymalnej długości (MLS, maximum length sequence) to sygnał pseudolosowy (o charakterze szumowym), składający się z sekwencji zer i jedynek, która jest być zapętlana. Jest to sygnał często stosowany w pomiarach akustycznych i spełnia on podany wyżej wymóg: jego autokorelacja liczona w sposób kołowy jest równa impulsowi jednostkowemu. Do wygenerowania sygnału MLS możemy użyć funkcji scipy.signal.max_len_seq. Jako parametr podaje się liczbę bitów $N$, wygenerowana sekwencja ma długość $2^N-1$. Długość sekwencji powinna być co najmniej dwukrotnie większa niż spodziewana długość odpowiedzi impulsowej.

Jako badany system wykorzystamy filtr FIR 30. rzędu. Generujemy sekwencję MLS na 8 bitach, o długości 255 próbek. Skalujemy ten sygnał tak, aby miał zerową wartość średnią i przyjmował wartości z zakresu od -1 do 1. Sygnał ten podajemy na wejście badanego układu. Następnie obliczamy korelację wzajemną sygnału wyjściowego i skalujemy ją, dzieląc przez długość sekwencji. Uzyskaną funkcję korelacji pokazuje poniższy wykres.

h = sig.firwin(31, (1000, 5000), pass_zero=False, fs=fs)
mls, _ = sig.max_len_seq(8)
mls = 2 * mls - 1
y = sig.lfilter(h, 1, mls)
h_est = sig.correlate(y, mls) / 255
plt.plot(h_est)
plt.xlabel('indeks próbki')
plt.ylabel('korelacja')
plt.title('F. korelacji odpowiedzi filtru z MLS')
plt.show()

Aby uzyskać pożądany wynik, musimy wziąć część funkcji autokorelacji od wartości środkowej i "przyciąć" do żądanej długości. Poniższy wykres pokazuje oryginalną odpowiedź impulsową filtru i odpowiedź estymowaną za pomocą MLS.

ind = len(h_est) // 2
plt.plot(h_est[ind:ind+31], label='estymowana')
plt.plot(h, label='oryginalna')
plt.xlabel('nr próbki')
plt.ylabel('korelacja')
plt.title('F. korelacji odpowiedzi filtru z MLS')
plt.legend()
plt.show()

Jak widać, w części środkowej uzyskaliśmy idealny wynik, drobne odstępstwa są widoczne na krańcach odpowiedzi impulsowej.

Do czego może nam się przydać odpowiedź impulsowa sali koncertowej? Mając tę odpowiedź, możemy wirtualnie "przepuścić" dowolny dźwięk przez ten system. Możemy wziąć nagranie dokonane w studio, obliczyć splot z pomierzoną odpowiedzią sali koncertowej i w ten sposób uzyskać brzmienie zbliżone do tego, jakie mielibyśmy odtwarzając dźwięk w tej sali koncertowej. Oczywiście w praktyce problem jest bardziej złożony. Mierzona odpowiedź impulsowa zawiera w sobie także właściwości mikrofonu i głośnika, zawiera szum tła, a sama sala nie jest dokładnie systemem LTI. Niemniej jednak, opisane rozwiązanie jest stosowane w procesorach efektów brzmieniowych i sprawdza się całkiem dobrze.

Podsumowanie ¶

Splot liniowy sygnałów jest liczony przez przesuwanie odwróconego w czasie sygnału względem drugiego, mnożenie pokrywających się próbek i sumowanie wyników.
Długość splotu liniowego jest równa sumie długości sygnałów pomniejszonej o 1: $L+M-1$.
Splot można obliczyć w dziedzinie częstotliwości, mnożąc transformaty sygnałów i obliczając transformatę odwrotną.
W systemach LTI (liniowych, o niezmienniczej odpowiedzi impulsowej), sygnał wyjściowy jest równy splotowi sygnału wejściowego z odpowiedzią impulsową systemu; transformata odpowiedzi impulsowej jest transmitancją systemu.
Korelacja wzajemna sygnałów jest miarą ich podobieństwa. Oblicza się ją podobnie jak splot, ale nie odwraca się sygnału w czasie. Położenie maksimum funkcji korelacji wskazuje na przesunięcie dające największe podobieństwo sygnałów.
Korelację można obliczyć w dziedzinie częstotliwości, mnożąc transformatę jednego sygnału przez zespolone sprzężenie drugiej transformaty i obliczając transformatę odwrotną. Uzyskujemy w ten sposób korelację kołową, aby uzyskać korelację liniową, sygnały należy uzupełnić do długości $L+M-1$.
Autokorelacja jest korelacją wzajemną sygnału z jego kopiami przesuniętymi w czasie. Autokorelację stosuje się do analizy sygnałów okresowych i pseudookresowych. Położenie maksimum funkcji autokorelacji wskazuje na okres sygnału.
Autokorelację można obliczyć w dziedzinie częstotliwości, obliczając kwadrat modułu trasformaty sygnału (periodogram), a następnie transformatę odwrotną.
Rozplot pozwala w niektórych przypadkach uzyskać oryginalny sygnał, przed jego zniekształceniem przez system o znanej odpowiedzi impulsowej. Rozplot odopowiada filtracji odwrotnej za pomocą filtru o transmitancji będącej odwrotnością transmitancji systemu zniekształcającego.
Estymację odpowiedzi impulsowej systemu LTI można przeprowadzić podając na wejście systemu sygnał MLS i obliczając jego korelację wzajemną z sygnałem wyjściowym.