%matplotlib inline
import numpy as np
import scipy.signal as sig
import scipy.interpolate as interp
import matplotlib.pyplot as plt
import math

np.set_printoptions(suppress=True)

from util import set_style
set_style()

Interpolacja ¶

Wartości sygnału cyfrowego są znane tylko w dyskretnych punktach czasu. W przypadku najczęściej spotykanych równomiernie próbkowanych sygnałów, odstęp między próbkami jest stały i równy okresowi próbkowania $T = 1/f_s$, gdzie $f_s$ oznacza częstotliwość próbkowania. W ciągu jednej sekundy mamy $f_s$ próbek sygnału. Zdarza się jednak, że potrzebujemy sygnału spróbkowanego z większą lub mniejszą częstotliwością, albo potrzebujemy wartości w dowolnej chwili, nie tylko w punktach wyznaczonych przez okres próbkowania.

Interpolacja (interpolation) sygnału cyfrowego polega na obliczeniu wartości sygnału w dowolnym punkcie pomiędzy próbkami. Trzeba tu wyraźnie zaznaczyć, że interpolacja nie tworzy nowych danych. Sygnał po interpolacji nie będzie identyczny z sygnałem spróbkowanym z większą częstotliwością. Interpolacja sygnału jedynie estymuje wartości w danych punktach na podstawie wartości istniejących próbek.

Zwiększenie częstotliwości próbkowania ¶

Mamy sygnał spróbkowany z częstotliwością $f_{s1}$. Chcemy mieć ten sygnał spróbkowany z większą częstotliwością $f_{s2} = L * f_{s1}$, gdzie $L$ jest liczbą naturalną większą od 1. Na początek, utworzymy funkcję generującą sygnał testowy – sumę sinusów pokrywających pasmo do częstotliwości Nyquista. Następnie generujemy sygnal testowy z częstotliwością próbkowania 16 kHz.

def generuj_sygnal(n, f0, fs):
    '''Generuje n próbek sygnału złożonego z sinusów, o częstotliwości pierwszego sinusa równej f0, 
    przy częstotliwości próbkowania fs.
    '''
    t = np.arange(n)
    sygnal = np.zeros(n)
    for i in range(1, fs // (2 * f0)):
        sygnal += np.sin(2 * np.pi * t * i * f0 / fs) / (i**2)
    return sygnal

x_16 = generuj_sygnal(2048, 1000, 16000)
f_16 = np.fft.rfftfreq(2048, 1 / 16000)
sp_16 = 20 * np.log10(np.abs(np.fft.rfft(x_16 * np.hamming(2048))) /1024)
plt.plot(f_16, sp_16)
plt.xlabel('Częstotliwość [Hz]')
plt.ylabel('Amplituda [dB]')
plt.title('Widmo sygnału testowego')
plt.show()

Chcemy teraz zwiększyć częstotliwość próbkowania do 48 kHz, czyli trzykrotnie. Potrzebujemy więc trzy razy więcej próbek. Na początek "doróbmy" dodatkowe próbki wstawiając po dwa zera między każdą parę próbek. Wykonując tę operację, zmieniamy energię sygnału na jednostkę czasu, musimy zatem przemnożyć amplitudy sygnału przez $L$=3. Popatrzmy na widmo powstałego sygnału.

x_48 = np.zeros(2048 * 3)
x_48[1::3] = 3 * x_16
f_48 = np.fft.rfftfreq(2048, 1 / 48000)
sp_48 = 20 * np.log10(np.abs(np.fft.rfft(x_48[:2048] * np.hamming(2048))) / 1024)
plt.plot(f_48, sp_48)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału po wstawieniu zer')
plt.show()

Co się tutaj stało? Pamiętajmy, że widmo sygnału cyfrowego obserwowane w szerokim zakresie częstotliwości jest powtarzalne. Dla sygnału rzeczywistego, mamy tutaj podstawową kopię widma w zakresie 0-8 kHz, odwróconą kopię w zakresie 8-16 kHz, normalną kopię w zakresie 16-24 kHz, itd. Przesuwając częstotliwość Nyquista z 8 na 24 kHz, objęliśmy zakresem trzy kopie widma.

Rozwiązanie problemu jest proste: należy usunąć składowe widmowe powyżej częstotliwości Nyquista oryginalnego sygnału za pomocą filtru dolnoprzepustowego. Ponieważ współczynnik $L$ zwiększania częstotliwości próbkowania jest całkowity, częstotliwość graniczna zawsze będzie ułamkiem częstotliwości próbkowania: $f_{s1}/(2L)$. Filtry tego typu nazywane są filtrami interpolacyjnymi (interpolation filter). Ponieważ część współczynników takiego typu jest zerowa, filtry interpolacyjne mogą być implementowane w prostszy sposób. Np. filtr FIR rzędu 30 dla $L=2$ ma współczynniki:

print(sig.firwin(31, 0.5))

[-0.0017004   0.          0.00293733 -0.         -0.00673009  0.
  0.01409389 -0.         -0.02678504  0.          0.04909896 -0.
 -0.09693833  0.          0.31561956  0.50080823  0.31561956  0.
 -0.09693833 -0.          0.04909896  0.         -0.02678504 -0.
  0.01409389  0.         -0.00673009 -0.          0.00293733  0.         -0.0017004 ]

Odfiltrujemy dwie niepotrzebne kopie widma za pomocą filtru FIR rzędu 300.

b_up3 = sig.firwin(301, 1 / 3)
x_48 = np.zeros(2048 * 3)
x_48[1::3] = x_16
x_48 = 3 * sig.filtfilt(b_up3, 1, x_48)
f_48 = np.fft.rfftfreq(2048, 1 / 48000)
sp_48 = 20 * np.log10(np.abs(np.fft.rfft(x_48[:2048] * np.hamming(2048))) / 1024)
plt.plot(f_48, sp_48)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału po interpolacji L=3')
plt.show()

Kopie widma zostały odfiltrowane. Zwróćmy uwagę na fakt, że powyżej częstotliwości 8 kHz brak jest składowych, które byłyby tam obecne dla sygnału spróbkowanego z większą częstotliwością. Interpolacja sygnału nie jest w stanie odtworzyć pełnego widma.

Sygnały w dziedzinie czasu wyglądają następująco.

fig, ax = plt.subplots(2)
ax[0].plot(x_16[200:300], label="fs=16 kHz")
ax[1].plot(x_48[600:700], label="fs=48 kHz")
for a in ax:
    a.set_xlabel('nr próbki')
    a.set_ylabel('amplituda')
    a.legend(loc='upper right')
ax[0].set_title('Sygnał testowy przed i po interpolacji')
plt.show()

Poniżej fragmenty sygnału (10 próbek dla pierwszego i 30 próbek dla drugiego, pokazujące wartości próbek, zaznaczonych kropkami na wykresach.

fig, ax = plt.subplots(2)
ax[0].plot(x_16[200:210], 'o-', label="fs=16 kHz")
ax[1].plot(x_48[600:630], 'o-', label="fs=48 kHz")
for a in ax:
    a.set_xlabel('nr próbki')
    a.set_ylabel('amplituda')
    a.legend()
ax[0].set_title('Sygnały przed i po interpolacji')
plt.show()

Interpolację sygnału ze stopniem $L$, która powoduje $L$-krotne zwiększenie liczby próbek, można traktować na dwa sposoby.

Jeżeli interpolowany sygnał jest odczytywany z częstotliwością próbkowania $L$ razy większą niż oryginalny, zachowana zostaje struktura widmowa sygnału.
Jeżeli częstotliwość próbkowania pozostaje taka sama, sygnał w dziedzinie czasu zostaje rozciągnięty $L$ razy, a widmo jest skalowane $L$ razy w dół, co powoduje $L$-krotne zmniejszenie częstotliwości składowych widmowych.

Operację wstawienia zer i filtracji dolnoprzepustowej możemy też przeprowadzić za pomocą funkcji scipy.signal.upfirdn:

x_48 = sig.upfirdn(b_up3, x_16, 3)

Jeżeli chcemy przeprowadzić tę operację bez konieczności projektowania filtru, możemy użyć funkcji scipy.signal.resample_poly. Jak widać poniżej, w tym przypadku pozostają szczątkowe składowe widmowe.

x_48 = sig.resample_poly(x_16, 3, 1)
f_48 = np.fft.rfftfreq(2048, 1 / 48000)
sp_48 = 20 * np.log10(np.abs(np.fft.rfft(x_48[:2048] * np.hamming(2048))) / 1024)
plt.plot(f_48, sp_48)
plt.xlabel('Częstotliwość [Hz]')
plt.ylabel('Amplituda [dB]')
plt.title('Widmo sygnału interpolowanego funkcją resample_poly')
plt.show()

Zmniejszenie częstotliwości próbkowania - decymacja ¶

W przypadku odwrotnym do przedstawionego powyżej, chcemy zmniejszyć częstotliwość próbkowania $M$ razy. Na przykład, mamy sygnał spróbkowany z częstotliwością 48 kHz, ale system wymaga podania sygnału o częstotliwości próbkowania 16 kHz. Tutaj sprawa jest nieco prostsza, zamiast dodawać nowe próbki, musimy usunąć część istniejących próbek. Operacja zmniejszania częstotliwości próbkowania nazywa się decymacją sygnału (decimation).

Zmniejszenie częstotliwości próbkowania trzy razy wymaga trzy razy mniejszej liczby próbek. Wybierzmy zatem co trzecią próbkę sygnału, odrzucając pozostałe, i popatrzmy na widmo.

x_48 = generuj_sygnal(2048 * 3, 1500, 48000)
x_16 = x_48[::3]
f_16 = np.fft.rfftfreq(2048, 1 / 16000)
sp_16 = 20 * np.log10(np.abs(np.fft.rfft(x_16 * np.hamming(2048))) / 1024)
plt.plot(f_16, sp_16)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału po wybraniu co 3. próbki')
plt.show()

Widać że składowych w widmie jest za dużo. W tym przypadku, kopie widma zostały przesunięte na nowe częstotliwości. Podstawowa część widma pozostaje w zakresie 0-24 kHz. Druga, odwrócona kopia widma przesuwa się z zakresu 24-48 kHz na -16-16 kHz. Podobnie inne kopie widma. Następuje aliasing (zakładkowanie) widma.

Rozwiązanie tego problemu wymaga zapewnienia, że nie ma w sygnale składowych powodujących aliasing. Konieczne jest zastosowanie filtru dolnoprzepustowego o częstotliwości granicznej równej $f_{s1}/(2M)$, gdzie $f_{s1}$ jest oryginalną częstotliwością próbkowania. Filtrację należy zastosować przed usunięciem próbek.

b_down = sig.firwin(201, 1 / 3)
x_16 = sig.filtfilt(b_down, 1, x_48)[::3]
f_16 = np.fft.rfftfreq(2048, 1 / 16000)
sp_16 = 20 * np.log10(np.abs(np.fft.rfft(x_16 * np.hamming(2048))) / 1024)
plt.plot(f_16, sp_16)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału po decymacji')
plt.show()

Możemy również użyć funkcji scipy.signal.decimate, jak również funkcji upfirdn i resample_poly z tego modułu.

x_16 = sig.upfirdn(b_down, x_48, 1, 3)

x_16 = sig.resample_poly(x_48, 1, 3)

Sygnały w dziedzinie czasu:

fig, ax = plt.subplots(2)
ax[0].plot(x_48[600:700], label="fs=48 kHz")
ax[1].plot(x_16[200:300], label="fs=16 kHz")
for a in ax:
    a.set_xlabel('nr próbki')
    a.set_ylabel('amplituda')
    a.legend(loc='upper right')
ax[0].set_title('Sygnał przed i po decymacji')
plt.show()

Decymację sygnału ze stopniem $M$, która powoduje $M$-krotne zmniejszenie liczby próbek, można traktować na dwa sposoby.

Jeżeli decymowany sygnał jest odczytywany z częstotliwością próbkowania $M$ razy mniejszą niż oryginalny, zachowana zostaje struktura widmowa sygnału.
Jeżeli częstotliwość próbkowania pozostaje taka sama, sygnał w dziedzinie czasu zostaje skrócony $M$ razy, a widmo jest skalowane $M$ razy w górę, co powoduje $M$-krotne zwiększenie częstotliwości składowych widmowych.

Przepróbkowanie sygnału ¶

W wielu przypadkach potrzebujemy zmienić częstotliwość próbkowania o stopień nie będący liczbą całkowitą. Dla przykładu, na płycie kompaktowej (CD) zastosowano częstotliwość próbkowania 44100 Hz, taki też standard przyjął się na komputerach osobistych. Z kolei w profesjonalnych zastosowaniach próbkuje się dźwięk z częstotliwością 48000 Hz.

Konwersję częstotliwości próbkowania, czyli przepróbkowanie sygnału (resampling) z wartości $f_{s1}$ na $f_{s2}$ można przeprowadzić stosując najpierw interpolację, a następnie decymację:

najpierw wstawiamy $L-1$ zerowych próbek pomiędzy każdą parę istniejących próbek,
następnie stosujemy filtr dolnoprzepustowy o częstotliwości granicznej równej mniejszej z wartości: $f_{s1}/(2L)$ lub $f_{s2}/(2M)$,
po czym wybieramy co $M$-tą próbkę.

Wartości $L$ i $M$ są liczbami całkowitymi spełniającymi warunek:

$$\frac{f_{s2}}{f_{s1}} = \frac{L}{M}$$

Warto zauważyć, że przy samej interpolacji stosuje się filtr dolnoprzepustowy po wstawieniu zer, a przy decymacji stosuje się filtr przed usunięciem próbek. Przy przepróbkowaniu mamy więc dwie filtracje jedna po drugiej. Nie ma sensu stosować dwóch osobnych filtrów dolnoprzepustowych, wystaczy jeden.

Aby wyznaczyć wartości $L, M$, umożliwiające konwersję częstotliwości próbkowania z 44100 Hz na 48000 Hz, szukamy największego wspólnego podzielnika. Użyjemy do tegu funkcji gcd z modułu math Pythona.

d = math.gcd(44100, 48000)
L = 48000 // d
M = 44100 // d
print('L = {}, M = {}'.format(L, M))

L = 160, M = 147

Jak widać, obie częstotliwości zostały wybrane dość niefortunnie. Wartości $L$ i $M$ są dużymi liczbami. Konwersja z 44,1 kHz na 48 kHz wymaga: wstawienia 160 zer między każdą parę próbek, zastosowania filtru dolnoprzepustowego o częstotliwości granicznej 22050 Hz, a następnie wybrania po jednej próbce co 147 próbek. Analogicznie odbywa się konwersja w drugą stronę. Wykonanie przepróbkowania sygnału wymaga dużej liczby operacji i dużego obszaru pamięci.

Do przepróbkowania sygnału możemy użyć wspomnianej już wyżej funkcji scipy.signal.resample_poly.

x_44 = generuj_sygnal(2048, 1000, 44100)
x_48 = sig.resample_poly(x_44, 160, 147)
f_48 = np.fft.rfftfreq(2048, 1 / 48000)
sp_48 = 20 * np.log10(np.abs(np.fft.rfft(x_48[:2048] * np.hamming(2048))) / 1024)
plt.plot(f_48, sp_48)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału po funkcji resample_poly')
plt.show()

Widać, że ta funkcja nie tłumi dostatecznie prążków powyżej 22050 Hz.

Przepróbkowanie sygnału za pomocą FFT ¶

Obserwując pokazane wcześniej wykresy widmowe można zauważyć, że operacja przepróbkowania sygnału poszerza widmo o "puste" wartości (interpolacja), albo je "ucina" (decymacja). Można więc przeprowadzić te operacje bezpośrednio w dziedzinie widma.

W przypadku interpolacji (zwiększania częstotliwości próbkowania) ze stopniem $L$:

obliczamy transformatę sygnału,
poszerzamy transformatę $L$-krotnie, dopisując zera na końcu,
wartości widma, za wyjątkiem "zerowej" (składowej stałej), mnożymy przez $L$ aby zachować energię sygnału,
obliczamy odwrotną transformatę.

Poniżej przykład dla trzykrotnego zwiększenia częstotliwości próbkowania z 16 kHz na 48 kHz. Stosujemy tu funkcje obliczające widmo sygnału rzeczywistego. Dla bloku 2048 próbek, funkcja rfft zwraca 1025 wartości, ostatnia z nich jest wartością dla częstotliwości Nyquista. Poszerzamy widmo o 2 * 1024 zerowych próbek i mnożymy pozostałe wartości widma przez 3, pozostawiając jednak niezmienioną składową stałą.

x_16 = generuj_sygnal(2048, 1000, 16000)
sp_16 = np.fft.rfft(x_16)
sp_48 = np.zeros(3 * (len(sp_16) - 1) + 1, dtype=np.complex)
sp_48[:len(sp_16)-1] = sp_16[:-1]
sp_48[1:len(sp_16)-1] *= 3
f_48 = np.fft.rfftfreq(3 * 2048, 1 / 48000)
plt.semilogy(f_48, np.abs(sp_48) / 1024)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału po interpolacji')
plt.show()

x_48 = np.real(np.fft.irfft(sp_48))
fig, ax = plt.subplots(2)
ax[0].plot(x_16[200:300], label="fs=16 kHz")
ax[1].plot(x_48[600:700], label="fs=48 kHz")
for a in ax:
    a.set_xlabel('nr próbki')
    a.set_ylabel('amplituda')
    a.legend(loc='upper right')
ax[0].set_title('Sygnał przed i po interpolacji FFT')
plt.show()

Uzyskany wynik jest poprawny, taki jak w przypadku poprzedniej metody interpolacji.

Decymację ze stopniem $M$ przeprowadza się następująco:

obliczamy transformatę sygnału,
wartości dla częstotliwości większych od 0 i mniejszych od częstotliwości Nyquista dzielimy przez $M$ aby zachować energię sygnału,
usuwamy część widma powyżej nowej częstotliwości Nyquista,
obliczamy odwrotną transformatę sygnału.

Poniżej przykład dla trzykrotnego zmniejszenia częstotliwości próbkowania z 48 kH na 16 kHz. Dla bloku 6144 (3*2048) próbek, funkcja rfft zwraca 3073 wartości. Po decymacji, długość transformaty dla funkcji rfft powinna wynosić 1025. Bierzemy pierwszych 1025 wartości do obliczenia odwrotnej transformaty, zmniejszając wartości widma trzykrotnie.

x_48 = generuj_sygnal(2048 * 3, 1000, 48000)
sp_48 = np.fft.rfft(x_48)
sp_16 = sp_48[:(len(sp_48) - 1) // 3 + 1]
#sp_16[-1] = sp_48[-1]
sp_16[1:-1] /= 3
f_16 = np.fft.rfftfreq(2048, 1 / 48000)
plt.semilogy(f_16, np.abs(sp_16) / 1024)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda')
plt.title('Widmo sygnału po decymacji')
plt.show()

x_16 = np.real(np.fft.irfft(sp_16))
fig, ax = plt.subplots(2)
ax[0].plot(x_48[600:700], label="fs=48 kHz")
ax[1].plot(x_16[200:300], label="fs=16 kHz")
for a in ax:
    a.set_xlabel('nr próbki')
    a.set_ylabel('amplituda')
    a.legend(loc='upper right')
ax[0].set_title('Sygnał przed i po decymacji FFT')
plt.show()

Przepróbkowanie ze współczynnikiem $L/M$ możemy przeprowadzić następująco:

bierzemy blok próbek o długości $k \cdot M$,
rozszerzamy albo obcinamy widmo tak jak wcześniej, aby uzyskać długość $k \cdot L$,
skalujemy wartości widmowe ze współczynnikiem $L/M$,
obliczamy odwrotną transformatę.

W powyższych przypadkach uprościliśmy nieco procedurę, pomijając składową widmową na częstotliwości Nyquista. Przy parzystej długości transformaty, składowa ta "nie ma pary" (jest rzeczywista i występuje raz w widmie), zatem przenosząc ją na niższą częstotliwość, powinniśmy ją rozdzielić na dwie składowe, dzieląc amplitudę przez 2.

Operację przepróbkowania sygnału za pomocą FFT możemy przeprowadzić za pomocą funkcji scipy.signal.resample. Pierwszym parametrem funkcji jest sygnał poddawany przepróbkowaniu, drugim – docelowa liczba próbek przepróbkowanego sygnału. Poniżej pokazano porównanie działania funkcji resample_poly i resample. Widać że w paśmie do pierwotnej częstotliwości Nyquista widmo jest identyczne, ale funkcja resample nie pozostawia składowych powyżej tej częstotliwości, podczas gdy funkcja resample nie stłumiła w pełni tych składowych.

x_44 = generuj_sygnal(2048, 1000, 44100)
x_48_poly = sig.resample_poly(x_44, 160, 147)
x_48 = sig.resample(x_44, int(len(x_44) * 160 / 147))
f_48 = np.fft.rfftfreq(2048, 1 / 48000)
sp_48_poly = 20 * np.log10(np.abs(np.fft.rfft(x_48[:2048] * np.hamming(2048))) / 1024)
plt.plot(f_48, 20 * np.log10(np.abs(np.fft.rfft(x_48_poly[:2048] * np.hamming(2048))) / 1024), label='resample_poly')
plt.plot(f_48, 20 * np.log10(np.abs(np.fft.rfft(x_48[:2048] * np.hamming(2048))) / 1024), label='resample')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Porównanie funkcji przepróbkowujących sygnał')
plt.legend()
plt.show()

Nadpróbkowanie ¶

Nadpróbkowanie (oversampling) oznacza spróbkowanie sygnału z częstotliwością wyższą niż docelowa. Na przykład, sygnały dźwiękowe mają pasmo częstotliwości ograniczone do ok. 20 kHz, ponieważ tyle wynosi zakres słyszenia osób o zdrowym słuchu. Dlatego praktycznie stosowane są częstotliwości próbkowania 44,1 kHz lub 48 kHz. Ale w profesjonalnym przetwarzaniu dźwięku próbkuje się sygnał z większą częstotliwością, np. 96 kHz lub 192 kHz (nadpróbkowanie razy 2 lub razy 4). Najczęściej nadpróbkowanie polega na zastosowaniu częstotliwości próbkowania będącej całkowitą wielokrotnością docelowej wartości. Sygnał jest przetwarzany w nadpróbkowanej formie, a następnie decymowany do docelowej częstotliwości próbkowania.

Nadpróbkowanie sygnału zwiększa liczbę operacji, które trzeba wykonać aby przetworzyć sygnał (mamy więcej próbek). W jakim celu stosuje się nadpróbkowanie? Są dwie główne korzyści:

zwiększenie rozdzielczości częstotliwościowej,
zmniejszenie ryzyka aliasingu.

Jak omówiono w rozdziale dotyczącym analizy widmowej, rozdzielczość częstotliwościowa analizy jest równa $f_s/N$, gdzie $N$ jest rozmiarem transformaty. Widać więc, że N-krotne nadpróbkowanie sygnału poprawia N-krotnie rozdzielczość częstotliwościową, a przez to zwiększa dokładność operacji przetwarzania sygnału wykonywanych w dziedzinie częstotliwości.

Nadpróbkowanie zwiększa również "margines" widma, redukując zjawisko aliasingu. W sygnale spróbkowanym z normalną częstotliwością może wystąpić aliasing w trakcie jego przetwarzania, jeżeli przeprowadzane są operacje prowadzące do przesunięcia składowych widmowych powyżej częstotliwości Nyquista. W sygnale nadpróbkowanym, składowe te znajdą się poniżej częstotliwości Nyquista nadpróbkowanego sygnału i mogą zostać odfiltrowane przed końcową decymacją.

Sygnał piłokształtny jest sygnałem harmonicznym o nieoggraniczonym widmie. Generowanie takiego sygnału "z definicji" niemal zawsze spowoduje powstanie aliasingu. Poniżej przykład sygnału piłokształtnego o częstotliwości 1250 Hz, wygenerowanego przy częstotliwości próbkowania 48 kHz. Widoczny jest bardzo silny aliasing.

pila1 = sig.sawtooth(2 * np.pi * np.arange(2048) * 1250 / 48000)
w_pila1 = 20 * np.log10(np.abs(np.fft.rfft(pila1 * np.hamming(2048))) / 1024)
f = np.fft.rfftfreq(2048, 1 / 48000)
plt.plot(f, w_pila1)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału piłokształtnego z aliasingiem')
plt.ylim(bottom=-50)
plt.show()

Poniżej ten sam przykład dla sygnału piłokształtnego wygenerowanego z częstotliwością próbkowania 192 kHz (czterokrotne nadpróbkowanie) i zdecymowanego. W poprzednim przypadku, wszystkie składowe powyżej 24 kHz powodują aliasing. Tutaj aliasing powstanie dla składowych powyżej 96 + (96 - 24) = 168 kHz, zatem margines widma jest znacznie większy. Aliasing nadal występuje, jednak amplitudy odbitych prążków są znacznie mniejsze niż przy braku nadpróbkowania (poniżej -50 dB).

pila2 = sig.sawtooth(2 * np.pi * np.arange(8192) * 1250 / 192000)
pila2 = sig.decimate(pila2, 4)
w_pila2 = 20 * np.log10(np.abs(np.fft.rfft(pila2 * np.hamming(2048))) / 1024)
plt.plot(f, w_pila2)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału piłokształtnego po nadpróbkowaniu i decymacji')
plt.ylim(bottom=-70)
plt.show()

Interpolacja wartości sygnału ¶

Poprzednie przykłady dotyczyły zwiększania lub zmniejszania częstotliwości próbkowania całego sygnału. W niektórych przypadkach chcemy mieć jednak możliwość odczytu dowolnej wartości sygnału "pomiędzy próbkami". Ponieważ sygnał cyfrowy jest zdefiniowany tylko w dyskretnych punktach czasu, musimy użyć interpolacji sygnału aby oszacować jego wartość dla ułamkowego indeksu próbki.

Załóżmy, że mamy zbiór próbek opisujący jeden okres sygnału okresowego. Niech tym sygnałem będzie suma pięciu składowych sinusoidalnych, a okres będzie opisany za pomocą 256 próbek.

N = 256
n = np.arange(N)
tablica = (np.sin(2 * np.pi * n / N) + np.sin(2 * np.pi * 2 * n / N) + np.sin(2 * np.pi * 3 * n / N) +
         np.sin(2 * np.pi * 4 * n / N) + np.sin(2 * np.pi * 5 * n / N)) / 5
plt.plot(n, tablica)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Próbki okresu sygnału')
plt.show()

Jeżeli będziemy odczytywać tę tablicę próbek w sposób cykliczny, otrzymamy sygnał okresowy. Indeksy odczytywanych próbek są kolejnymi liczbami naturalnymi. Po dojściu do końca tablicy, "zawijamy" indeks, co odpowiada operacji modulo długość tablicy. W języku Python, operator modulo oznacza się symbolem procentu (%).

indeks = np.arange(1500) % N
sygnal = tablica[indeks]
plt.plot(sygnal)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał wygenerowany przez powtarzanie okresu')
plt.show()

Jeżeli tablica próbek o długości $N$ jest odczytywana z częstotliwością $f_s$, to jeden okres sygnału trwa $N \cdot f_s$ sekund. Zatem wynikowa czestotliwość odczytywanego sygnału jest równa:

$$f = \frac{f_s}{N}$$

Dla $N$=256 i $f_s$=1 kHz, otrzymujemy $f$=3,91 Hz.

Co się stanie, gdy odczytamy z tablicy co drugą próbkę? Odpowiada to przypadkowi, gdy tablica zawiera $N/2$ próbek okresu. Odczytany okres sygnału ulega dwukrotnemu skróceniu, a zatem częstotliwość sygnału wzrasta dwukrotnie. Jeżeli krok przesuwania indeksu próbki oznaczymy jako $k$, to wynikowa częstotliwość sygnału wynosi:

$$f = k \cdot \frac{f_s}{N}$$

k = 2
indeks = (k * np.arange(1500)) % N
sygnal = tablica[indeks]
plt.plot(sygnal)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał przy odczycie tablicy z krokiem 2')
plt.show()

Należy pamiętać, że jeżeli sygnał zapisany tablicy ma jakiekolwiek składowe widmowe powyżej częstotliwości $f_s/4$, odczytanie go z dwukrotnie większą częstotliwością spowoduje aliasing widma.

W praktycznych zastosowaniach potrzebujemy uzyskać dowolną częstotliwość sygnału odczytywanego z tablicy. Jeżeli potrzebujey częstotliwości $f$=1 Hz, to z podanej wyżej zależności można obliczyć krok przesuwania indeksu:

fs = 1000
f = 1
k = N * f / fs
print('Krok indeksu dla f={} Hz: k={}'.format(f, k))

Krok indeksu dla f=1 Hz: k=0.256

Widzimy więc, że krok odczytu próbki jest ułamkowy. Aby uzyskać wartości sygnału "pomiędzy próbkami", potrzebna jest interpolacja.

Najprostszym podejściem jest zaokrąglenie indeksu do najbliższej wartości całkowitej. Operację tę nazywa się interpolacją najbliższego sąsiada (nearest neighbour) lub interpolacją zerowego rzędu.

n = (k * np.arange(fs // f)) % N
indeks = (np.round(n).astype(np.int)) % N
sygnal = tablica[indeks]
n2 = np.arange(fs // f)
wzor = (np.sin(2 * np.pi * n2 * f / fs) + np.sin(2 * np.pi * 2 * n2 * f / fs) + np.sin(2 * np.pi * 3 * n2 * f / fs) +
         np.sin(2 * np.pi * 4 * n2 * f / fs) + np.sin(2 * np.pi * 5 * n2 * f / fs)) / 5
e0 = np.sqrt(np.mean((wzor - sygnal)**2))
print('Interpolacja n. sąsiada: e =', e0)

plt.plot(wzor, label='wzorcowy')
plt.plot(sygnal, '+', label='interpolowany')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.xlim(50, 80)
plt.ylim(0.7, 0.8)
plt.title('Interpolacja najbliższego sąsiada')
plt.legend()
plt.show()

Interpolacja n. sąsiada: e = 0.00743366188857

Widać, że metoda najbliższego sąsiada powoduje występowanie próbek o tej samej amplitudzie. Błąd interpolacji został obliczony jako pierwiastek ze średniej kwadratów róznic sygnałów: interpolowanego oraz wygenerowanego bezpośrednio. Dokładność interpolacji tą metodą zależy od liczby próbek w tablicy (im większy odstęp między próbkami, tym większy błąd) i od kroku odczytu. Metoda ta może być wystarczająca, gdy okres sygnału jest reprezentowany przez dużą liczbę próbek (np. 2048) i gdy częstotliwość jest zwiększana (krok większy od 1).

Interpolacja liniowa ¶

Interpolacja liniowa (linear interpolation), czyli interpolacja pierwszego rzędu, oblicza wartość sygnału zakładając, że kształt sygnału pomiędzy dwoma znanymi punktami jest liniowy. Jest to prosta metoda interpolacji, ale jest ona dokładniejsza niż metoda najbliższego sąsiada. Załóżmy, że znamy wartości sygnału liniowego: (0, 0) i (1, 1), chcemy obliczyć wartość w punkcie 0,5. Metoda najbliższego sąsiada zaokrągli indeks i otrzymamy wartość (0,5; 0) – błąd będzie równy 0,5. Metoda interpolacji liniowej obliczy średnią obu wartości i otrzymamy tu poprawny wynik (0,5; 0,5). Oczywiście, błąd będzie tym większy, im bardziej nieliniowy jest interpolowany przebieg.

Znając wartości sygnału $(x_1, y_1)$ i $(x_2, y_2)$, a chcąc obliczyć wartość sygnału w punkcie $x$, $x_1 \le x \lt x_2$, możemy użyć prostej zależności (zakładamy tu, że $x$ jest indeksem próbki, czyli $x_2-x_1=1$):

$$x = y_1 + (x-x_1)(y_2-y_1)$$

Do przeprowadzenia interpolacji liniowej użyjemy funkcji numpy.interp. Pierwszym argumentem są punkty interpolacji (ułamkowe indeksy), drugim – punkty $x$, w których sygnał jest określone (indeksy próbek w tablicy), trzecim – wartości $y$ sygnału w tych punktach (wartości próbek w tablicy). Podanie parametru period umożliwia poprawną interpolację w buforze kołowym.

n = np.arange(N)
indeks = (k * np.arange(fs // f)) % N
sygnal = np.interp(indeks, n, tablica, period=N)
e1 = np.sqrt(np.mean((wzor - sygnal)**2))
print('Interpolacja liniowa:    e =', e1)
print('Interpolacja n. sąsiada: e =', e0)

plt.plot(wzor, label='wzorcowy')
plt.plot(sygnal, '+', label='interpolowany')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.xlim(50, 80)
plt.ylim(0.7, 0.8)
plt.title('Interpolacja liniowa')
plt.legend()
plt.show()

Interpolacja liniowa:    e = 0.000243262879532
Interpolacja n. sąsiada: e = 0.00743366188857

Widoczne jest, że próbki leżą znacznie bliżej idealnej krzywej, nie ma "schodków", które występowały dla metody najbliższego sąsiada, jednak nadal interpolowane wartości nie pokrywają się w pełni z przebiegiem sygnału. Błąd interpolacji jest znacznie mniejszy.

Metoda interpolacji liniowej jest prosta obliczeniowo, dobrze sprawdza się w przypadkach, gdy okres sygnału jest reprezentowany przez dużą liczbę próbek i gdy w sygnale nie występują gwałtowne zmiany kształtu przebiegu czasowego.

Interpolacja krzywymi sklejanymi ¶

W przypadku interpolacji liniowej, pary punktów są łączone prostymi, czyli wielomianami pierwszego rzędu. Możliwe jest też zastosowanie wielomianów wyższego rzędu. Do trójek punktów można dopasować wielomian drugiego rzędu (parabolę), do czterech punktów – wielomian trzeciego rzędu, itd. Jest jednak warunek: poszczególne krzywe dopasowane do różnych punktów muszą się łączyć ze sobą w sposób ciągły. Po "sklejeniu" odcinków krzywych muszą on tworzyć ciągły przebieg, bez załamań. W praktyce, dla odcinków krzywych stopnia $n$, same odcinki, jak też ich pochodne do stopnia $n-1$ włącznie, muszą zachowywać ciągłość w punkcie połączenia. Powstałe w ten sposób przebiegi noszą nazwę krzywych sklejanych n-tego stopnia (splines, nieformalnie po polsku używa się terminu "splajny").

Do przeprowadzenia interpolacji sygnału za pomocą krzywych sklejanych musimy sięgnąć do pakietu scipy.interp (tutaj dołączony pod nazwą interp) i funkcji interp1d. Pierwsze dwa parametry funkcji to znane wartości $x$ i $y$ sygnału. Trzeci parametr, kind oznacza sposób przeprowadzenia interpolacji. Podając tu nazwę 'nearest' lub 'linear', uzyskamy odpowiednio interpolację metodą najbliższego sąsiada lub liniową. Zastosowanie krzywych sklejanych wymaga podania nazwy: 'zero' (interpolacja stopnia zerowego), 'slinear' (pierwszego, liniowa), 'quadratic' (drugiego, kwadratowa) lub 'cubic' (trzeciego, kubiczna). Najczęściej w praktyce używane są krzywe sklejane stopnia trzeciego. Zerowy stopień krzywej oznacza, że nie ma wymagań co do ciągłości krzywej (może ona być łamana).

Funkcja interp1d działa nieco inaczej od funkcji numpy.interp, ponieważ zwraca funkcję, której możemy użyć do obliczania wartości interpolowanego sygnału. Nie ma ona również parametru period, więc musimy rozszerzyć tablicę próbek dodając kopię pierwszego elementu na koniec tablicy, aby umożliwić interpolację między ostatnią a pierwszą próbką okresu. Poniżej przykład dla interpolacji krzywymi sklejanymi trzeciego rzędu.

indeks = (k * np.arange(fs // f)) % N
tablica2 = np.append(tablica, tablica[0])
interpolator = interp.interp1d(np.arange(N + 1), tablica2, 'cubic', assume_sorted=True)
sygnal = interpolator(indeks)
e3 = np.sqrt(np.mean((wzor - sygnal)**2))
print('Interpolacja kubiczna:   e =', e3)
print('Interpolacja liniowa:    e =', e1)
print('Interpolacja n. sąsiada: e =', e0)

plt.plot(wzor, label='wzorcowy')
plt.plot(sygnal, '+', label='interpolowany')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Interpolacja krzywymi sklejanymi 3. stopnia')
plt.xlim(50, 80)
plt.ylim(0.7, 0.8)
plt.legend()
plt.show()

Interpolacja kubiczna:   e = 6.01172069485e-08
Interpolacja liniowa:    e = 0.000243262879532
Interpolacja n. sąsiada: e = 0.00743366188857

Interpolowane punkty są dobrze dopasowane do przebiegu funkcji, błąd interpolacji jest bardzo mały. Porównajmy jeszcze dokładność wszystkich dostępnych metod interpolacji.

i = np.arange(len(tablica2))
for typ in ('nearest', 'linear', 'zero', 'slinear', 'quadratic', 'cubic'):
    interpolator = interp.interp1d(i, tablica2, typ, assume_sorted=True)
    sygnal = interpolator(indeks)
    e = np.sqrt(np.mean((wzor - sygnal)**2))
    print('{:10} e = {}'.format('{}:'.format(typ), e))

nearest:   e = 0.007433661888570736
linear:    e = 0.00024326287953216756
zero:      e = 0.014690824167433392
slinear:   e = 0.00024326287953216615
quadratic: e = 2.5983045357530926e-06
cubic:     e = 6.011720694847838e-08

Najmniej dokładna jest metoda krzywych sklejanych nieciągłych, zerowego rzędu, w następnej kolejności metoda najbliższego sąsiada. Interpolacja liniowa, zwykła i krzywymi sklejanymi, daje identyczny wynik. Krzywe sklejane rzędu drugiego i trzeciego dają wielokrotnie mniejszy błąd interpolacji niż poprzednie metody, nawet w przypadku badanego sygnału, który ma bardzo łagodnie zmieniający się przebieg.

Interpolacja za pomocą funkcji sinc ¶

Jedną z najdokładniejszych metod interpolacji sygnału jest ta wykorzystująca funkcje sinc. Funkcja ta, której nazwa pochodzi od łacińskiego terminu sinus cardinalis, jest równa wartości sinusa argumentu, podzielonej przez ten argument, przy czym dla zerowego argumentu, wartość funkcji sinc jest równa jedności. Transformata Fouriera funkcji okna prostokątnego jest proporcjonalna do funkcji sinc. Zastosowanie interpolacji z użyciem funkcji sinc jest optymalnym podejściem do rekonstrukcji sygnału cyfrowego.

Funkcja numpy.sinc oblicza wartość funkcji według zależności: $sinc(x) = sin(\pi x)/(\pi x)$. Poniżej na wykresie pokazano fragment sygnału zapisanego w tablicy (10 próbek) oraz 10 funkcji sinc ustawionych w taki sposób, że maksimum każdej funkcji przypada dla indeksu jednej z próbek, a amplituda funkcji sinc w tym miejscu jest równa wartości tej próbki. Dla każdej z próbek sygnału, oznaczonej symbolem x na wykresie, jedna z funkcji sinc ma maksimum, pozostałe funkcje sinc mają wartość zero. Czyli: sumując wartości wszystkich funkcji sinc dla każdego indeksu reprezentowanego przez próbkę, otrzymamy dokładnie wartości sygnału.

Interpolacja za pomocą funkcji sinc (potocznie nazywana "interpolacją sinkową") polega na zsumowaniu wszystkich funkcji sinc dla dowolnie wybranego indeksu. Na wykresie pokazano pionową linią indeks równy 1,7. Sumując wartości wszystkich funkcji sinc w tym punkcie otrzymamy interpolowaną wartość sygnału dla tego indeksu.

syg = tablica[60:70]
plt.plot(np.arange(0, 10), syg, '#c0c0c0')
plt.plot(np.arange(0, 10), syg, 'x', c='k')
t = np.arange(0, 10, 0.1)
for i, x in enumerate(syg):
    plt.plot(t, x * np.sinc(np.arange(-i, 10 - i, 0.1)))
plt.axvline(1.7, color='#808080', lw=1)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Ilustracja metody interpolacji funkcjami sinc')
plt.show()

Metoda interpolacji funkcjami sinc jest najbardziej złożoną obliczeniowo z przedstawionych tutaj, ze względu na dużą liczbę obliczeń i duże wymagania pamięciowe, ponieważ trzeba uwzględnić osobną funkcję sinc dla każdej próbki w sygnale.

Niestety, nie ma w pakietach NumPy i SciPy funkcji dokonującej interpolacji opisywaną metodą. Poniższa funkcja pochodzi ze strony https://gist.github.com/endolith/1297227. Funkcja numpy.tile tworzy macierz powielając argument zadaną liczbę razy. Pierwsze wywołanie np.tile bierze wektor punktów interpolacji i powiela go rzędami tyle razy, ile jest próbek w tablicy. Drugie wywołanie funkcji bierze wektor kolumnowy indeksów próbek w tablicy i powiela go kolumnami tyle razy, ile jest interpolowanych punktów. Obie macierze są odejmowane od siebie, co odpowiada przesuwaniu funkcji sinc. Następnie wykonywane jest mnożenie macierzowe (operator @) wektora wartości próbek w tablicy przez wartości sinc dla obliczonej uprzednio macierzy, operacja ta odpowiada skalowaniu funkcji sinc przez wartości próbek i sumowanie wszystkich funkcji sinc. Ten algorytm można też zapisać za pomocą pętli for, jednak wersja macierzowa ma znacznie krótszy zapis i jest szybsza (ale zajmuje znacznie więcej pamięci).

def sinc_interp(x, xp, fp):
    '''Interpoluje sygnał zdefiniowany w punktach xp, z wartościami fp, obliczając wartości w punktach x.'''
    sincM = np.tile(x, (len(xp), 1)) - np.tile(xp[:, np.newaxis], (1, len(x)))
    y = fp @ np.sinc(sincM)
    return y

sygnal = sinc_interp(indeks, np.arange(N), tablica)
esinc = np.sqrt(np.mean((wzor - sygnal)**2))
print('Interpolacja sinc:       e =', esinc)

plt.plot(wzor, label='wzorcowy')
plt.plot(sygnal, '+', label='interpolowany')
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Interpolacja sinc')
plt.xlim(50, 80)
plt.ylim(0.7, 0.8)
plt.legend()
plt.show()

Interpolacja sinc:       e = 0.000549854792911

Otrzymaliśmy gorszą dokładność niż dla krzywych sklejanych. Wynika to ze skończonej precyzji obliczeń i kumulowania błędów przy mnożeniu dużych macierzy. Jednak dla sygnałów znacznie bardziej złożonych niż ten wykorzystany w przykładach, interpolacja funkcjami sinc może dać lepszą dokładność.

Podsumowanie ¶

$L$-krotne zwiększenie częstotliwości próbkowania (interpolację) uzyskamy wstawiając $L-1$ zer między każdą parę próbek i stosując filtr dolnoprzepustowy o częstotliwości granicznej równej oryginalnej częstotliwości Nyquista.
$M$-krotne zmniejszenie częstotliwości próbkowania (decymację) uzyskamy stosując filtr dolnoprzepustowy o częstotliwości granicznej równej docelowej częstotliwości Nyquista i wybierając z sygnału co $M$-tą próbkę.
Przepróbkowanie sygnału w stosunku $L/M$ możemy wykonać stosując interpolację stopnia $L$, a następnie decymację stopnia $M$, przy czym wystarczy jeden raz zastosować filtr dolnoprzepustowy.
Zmianę częstotliwości próbkowania sygnału można przeprowadzić również przy pomocy transformacji Fouriera (FFT).
Nadpróbkowanie sygnału jest przydatne w jego przetwarzaniu, zwiększa rozdzielczość częstotliwościową i zmniejsza efekty ewentualnego aliasingu.
Jeżeli sygnał okresowy po interpolacji/decymacji jest odczytywany z oryginalną częstotliwością próbkowania, uzyskujemy zmniejszenie/zwiększenie częstotliwości sygnału.
Oczytanie wartości próbek sygnału dla niecałkowitych indeksów jest możliwe dzięki interpolacji.
Interpolacja metodą najbliższego sąsiada jest najmniej dokładna, można ją stosować gdy dostatecznie dużo próbek opisuje kształt sygnału.
Interpolacja liniowa jest mało dokładna, ale w niektórych przypadkach wystarczająca, gdy przebieg sygnału nie zmienia się gwałtownie.
Interpolacja krzywymi sklejanymi trzeciego stopnia wymaga wykonania większych obliczeń, ale potrafi dokładnie interpolować złożone kształty przebiegu sygnału.
Interpolacja funkcjami sinc jest optymalną metodą rekonstrukcji sygnału, ale jest złożona obliczeniowo.