In [1]:
%matplotlib inline
import numpy as np
import scipy.signal as sig
from scipy.io import wavfile
import matplotlib.pyplot as plt
import math
import time
In [2]:
from util import set_style
set_style()
Out[2]:
In [3]:
# częstotliwość próbkowania
fs = 48000

Generowanie cyfrowych sygnałów

W poprzednich rozdziałach zajmowaliśmy się przetwarzaniem istniejących sygnałów. W tym rozdziale omówione zostaną metody wytwarzania sygnałów, które są najczęściej używane w technice cyfrowej.

Generowanie sygnału sinusoidalnego

Sygnał sinusoidalny jest jednym z najczęściej stosowanych w przetwarzaniu sygnałów, generowaliśmy go wielokrotnie w przykładach zamieszczonych w poprzednich rozdziałach. Cyfrowe wytwarzanie sygnału sinusoidalnego nie jest jednak tak prostym zadaniem, jak mogłoby się wydawać.

Chcemy wytworzyć sygnał sinusoidalny o częstotliwości $f$, np. 200 Hz, przy założonej częstotliwości próbkowania $f_s$, np. 48 kHz. Okres sygnału wynosi $T=1/f_s$. W ciągu sekundy mamy $f_s$ próbek sygnału. Stąd wynika, że liczba próbek $N$ przypadająca na jeden okres sygnału wynosi:

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

Dla $f$=200 Hz i $f_s$=48 kHz, $N$=240 próbek. Jeden okres sygnału odpowiada zmianie pulsacji (kąta fazowego) sygnału o $2\pi$. Zatem zmiana pulsacji w ciągu jednego okresu próbkowania (z próbki na próbkę) wynosi:

$$\Delta\varphi = \frac{2\pi}{N} = \frac{2\pi f}{f_s}$$

Pulsacja $n$-tej próbki sygnału będzie równy $n\cdot\Delta\varphi$. Faza sygnału "zawija się" wokół wartości $2\pi$, więc możemy obliczyć wartości modulo $2\pi$, chociaż nie jest to konieczne. Poniższy przykład generuje 1000 próbek sygnału sinusoidalnego o częstotliwości 200 Hz, przy częstotliwości próbkowania 48 kHz.

In [4]:
f = 200
n = np.arange(1000)
faza = (n * 2 * np.pi * f / fs) % (2 * np.pi)
plt.plot(n, faza)
plt.xlabel('nr próbki')
plt.ylabel('kąt fazowy [rad]')
plt.title('Sygnał fazowy')
plt.show()

Powyższy przebieg nazywa się sygnałem fazowym (phasor), reprezentuje on zmiany kąta fazowego w czasie.

Wygenerowanie sygnału sinosuidalnego polega po prostu na obliczeniu funkcji sin z kąta fazowego dla każdej próbki.

In [5]:
sinus = np.sin(faza)
plt.plot(n, sinus)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał sinusoidalny')
plt.show()

Generowanie sygnału sinusoidalnego za pomocą szeregu Taylora

Powyższa metoda wytwarzania sygnału sinusoidalnego jest prosta, o ile mamy do dyspozycji funkcję obliczającą sinus argumentu. Nie zawsze tak jest. Implementując algorytmy na procesorach sygnałowych, mamy właściwie do dyspozycji tylko trzy operacje matematyczne: dodawanie, odejmowanie i mnożenie. W jaki sposób obliczyć wartość funkcji trygonometrycznych za pomocą tych operacji?

Najczęściej stosowana metoda stosuje rozwinięcie funkcji trygonometrycznych w szereg Taylora. Dla funkcji sinus, wzór ma postać:

$$sin(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + ...$$

Szereg jest nieskończony, więc w praktyce musimy go ograniczyć do pewnej liczby elementów. Mamy zatem metodę obliczania sinusa za pomocą dodawania, odejmowania i mnożenia. Współczynniki $1/(2k+1)!$ są stałe, można je obliczyć i zapisać w kodzie programu, więc operacja dzielenia nie jest potrzebna. Poniższy kod realizuje obliczanie funkcji sinus biorąc pod uwagę k_max składników (maksymalna potęga równa 2*k_max – 1).

In [6]:
def sinus_taylor(x, k_max):
    '''Oblicza wartości funkcji sin dla wartości kąta fazowego x, z szeregu Taylora zawierającego k_max składników.'''
    y = np.copy(x)
    xx = np.copy(x)
    x2 = x * x
    silnia = 1.0
    for k in range(1, k_max + 1):
        xx *= x2
        silnia = silnia * (2 * k) * (2 * k + 1)
        y += (-1)**k * xx / silnia
    return y

Poniżej pokazano sygnał wygenerowany dla k_max = 9.

In [7]:
sinust = sinus_taylor(faza, 9)
plt.plot(sinust)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Szereg Taylora dla sin(x) - 9 składowych')
plt.show()
print('Błąd:', np.sqrt(np.sum((sinust[:240] - sinus[:240])**2)))

Błąd: 0.00237475024966

Poniżej wykreślono wartości błędu obliczania okresu sygnału sinusoidalnego metodą szeregu Taylora, w porównaniu z tradycyjną metodą. Liczba składowych większa niż 7 nie powoduje już znaczącego spadku błędu, co wynika z coraz większych mianowników kolejnych składników.

In [8]:
k = range(3, 21)
err = [np.sqrt(np.sum((sinus_taylor(faza[:240], kk) - sinus[:240])**2)) for kk in k]
plt.plot(k, err)
plt.xlabel('liczba składowych w szeregu Taylora')
plt.ylabel('błąd obliczeń')
plt.title('Błąd obliczeń w zależności od liczby składowych')
plt.xticks(np.arange(2, 21, 2))
plt.show()

Filtr jako oscylator

Pokażemy teraz ciekawostkę: filtr IIR zaprojektowany na granicy stabilności, który działa jako oscylator, generujący sygnał sinusoidalny. Jest to filtr drugiego rzędu, o transmitancji:

$$H(z)=\frac{1}{1+a_1 z^{-1} + z^{-2}}$$

Wartość współczynnika $a_1$ zależy od pożądanej częstotliwości $f$:

$$a_1 = -2 \cos(2\pi f/f_s)$$

Filtr musi posiadać niezerowe warunki początkowe, zgodnie z zależnością:

$$y(-2) = -\sin(2\pi f/f_s)$$ $$y(-1) = 0$$

Następnie wystarczy "uruchomić" filtr, nie podając żadnego sygnału na jego wejście. Filtr posiada parę biegunów położonych na okręgu jednostkowym. Przy podanych wyżej założeniach, filtr "wpada" w oscylacje o stałej amplitudzie, które mają charakter sygnału sinusoidalnego o pożądanej częstotliwości.

W poniższej implementacji, początkowy stan filtru został obliczony za pomocą funkcji lfiltic. Następnie użyliśmy funkcji lfilterdo obliczenia przebiegu. Należało tu sztucznie podać na wejście filtru sygnał złożony z samych zer.

In [9]:
f = 200
b = [1.0]
a = [1.0, -2 * np.cos(2 * np.pi * f / fs), 1.0]
zi = sig.lfiltic(b, a, [0, -np.sin(2 * np.pi * f / fs)])
print('Moduły biegunów:', [np.abs(k) for k in sig.tf2zpk(b, a)[1]])
sinusf = sig.lfilter(b, a, np.zeros(1000), zi=zi)[0]
plt.plot(sinusf)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Odpowiedź oscylatora')
plt.show()
print('Błąd:', np.sqrt(np.sum((sinusf[:240] - sinus[1:241])**2)))

Moduły biegunów: [1.0, 1.0]

Błąd: 1.17544378039e-12

W wygenerowanym sygnale "brakuje" pierwszej, zerowej próbki. Poza tym, otrzymaliśmy bardzo dokładny przebieg, różniący się od wygenerowanego "z definicji" o $10^{-12}$. Metoda jest prosta w implementacji i nie wymaga obliczania funkcji trygonometrycznych, za wyjątkiem stanu początkowego i współczynnika $a_1$. Wadą algorytmu jest to, że filtr jest wrażliwy na błędy precyzji obliczeń w niektórych implementacjach, np. stałoprzecinkowych, filtr może wtedy stać się niestabilny lub zmienić częstotliwość oscylacji.

Poniżej przedstawiono generator (zarówno w znaczeniu generatora sygnału, jak i struktury języka Python), który oblicza próbki sygnału sinusoidalnego bez użycia modułów NumPy i SciPy. Używa on funkcji sin i cos z modułu math z biblioteki standardowej języka Python, ale zamiast nich możemy użyć funkcji obliczających te wartości za pomocą szeregu Taylora. Generator jest uruchamiany instrukcją next. Po napotkaniu instrukcji yield, generator zwraca wartość i zawiesza działanie, zapamiętując swój stan. Po kolejnym wywołaniu next, generator wznawia działanie od ostatnio zapamiętanego miejsca, aż do kolejnej instrukcji yield. Pierwsza zwracana wartość jest zerowa, kolejne są obliczane podaną wyżej metodą.

In [10]:
def singen(f, fs):
    '''Generuje próbki sinusa o częstotliwości f, z częstotliwością próbkowania fs.'''
    yield 0
    y1 = 0
    y2 = -math.sin(2 * np.pi * f / fs)
    a1 = 2 * math.cos(2 * np.pi * f / fs)
    while True:
        y = a1 * y1 - y2
        yield y
        y2 = y1
        y1 = y
In [11]:
# Przykład użycia
gen = singen(200, fs)
sinusf2 = [next(gen) for i in range(1000)]

Opis algorytmu pochodzi z dokumentu: Yao-Ting Cheng, TMS320C62x Algorithm: Sine Wave Generation, Texas Instruments Application Report SPRA708.

Generowanie sygnałów harmonicznych

Sygnałami harmonicznymi nazywamy takie, które można potraktować jako sumę sygnałów sinusoidalnych o częstotliwościach będących całkowitą wielokrotnością częstotliwości podstawowej (fundamental frequency). Te składowe sygnały nazywamy harmonicznymi sygnału, a ich suma tworzy szereg harmoniczny. Poniżej przedstawiono trzy najczęściej stosowane, "regularne" sygnały harmoniczne. Sygnał sinusoidalny również zalicza się do sygnałów harmonicznych, zawiera on tylko składową podstawową. Warto dodać, że dźwięki większości instrumentów muzycznych w stanie ustalonym są w przybliżeniu sygnałami harmonicznymi.

Sygnał piłokształtny (sawtooth) ma przebieg narastający liniowo od wartości $-A$ do $A$, gdzie $A$ jest amplitudą sygnału, po czym następuje przeskok z powrotem do wartości $-A$. Jeżeli popatrzymy na zamieszczony wcześniej wykres sygnału fazowego, zobaczymy, że jest on właśnie piłokształtny. Możemy więc wygenerować sygnał piłokształtny w podobny sposób, ustawiając wartość maksymalną (modulo) na 2 i odejmując 1, uzyskując w ten sposób sygnał o wartościach od -1 do 1. Dodajemy 1 przed obliczeniem modulo, aby przebieg zaczynał się w zerze od zerowej amplitudy (przesuwamy fazę o połowę okresu).

In [12]:
f = 200
n = np.arange(1000)
pila = ((n * 2 * f / fs + 1) % 2) - 1
plt.plot(n, pila)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał piłokształtny')
plt.show()

Sygnał prostokątny (square) przyjmuje tylko dwie wartości: $-A$ i $A$. Klasyczny sygnał prostokątny, który właściwie należałoby nazywać kwadratowym, ma równe półokresy. Sygnał taki można wygenerować za pomocą sygnału piłokształtnego, dokonując jego progowania. Amplitudy dodatnie są zamieniane na wartości $+A$, ujemne na $-A$. Do progowania możemy użyć warunku logicznego $y > 0$, zwracającego wartości logiczne True i False. Wartości logiczne można rzutować na wartości liczbowe (numpy.astype), odpowiednio 1 i 0. Aby zamienić je na wartości +1 i -1, trzeba je przemnożyć przez 2 i odjąć 1.

In [13]:
prostokat = (pila >= 0).astype(np.int) * 2 - 1
plt.plot(n, prostokat)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał prostokątny')
plt.show()

Powyższy sygnał prostokątny ma wspołczynnik kształtu, nazywany też współczynnikiem wypełnienia (fill rate) równy 0,5, czyli połowę okresu zajmuje część dodatnia. Fala prostokątna może też mieć inny współczynnik kształtu. Można go uwzględnić korygując przesunięcie sygnału piłokształtnego. Poniżej przedstawiono przykład dla współczynnika kształtu 0,8.

In [14]:
b = 0.8
faza = ((n * 2 * f / fs + 1) % 2) + (b - 1)
prostokat08 = (faza >= 0).astype(np.int) * 2 - 1
plt.plot(n, prostokat08)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał prostokątny, wsp. kształtu 0,8')
plt.show()

Sygnałem trójkątnym (triangle) nazywamy taki, w którym amplituda narasta liniowo do połowy okresu, po czym maleje liniowo do końca okresu. Jedną z metod wygenerowania takiego sygnału jest użycie sygnału piłokształtnego, obliczenie jego modułu, przeskalowanie razy 2 i odjęcie 1. Sygnał piłokształtny jest przesuwany o 1/4 okresu aby wynikowy sygnał trójkątny zaczynał się w zerze.

In [15]:
faza = ((n * 2 * f / fs - 0.5) % 2) - 1
trojkat = 2 * np.abs(faza) - 1
plt.plot(n, trojkat)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał trójkątny')
plt.show()

Sygnał trójkątny może też być obliczony z sygnału prostokątnego poprzez jego całkowanie, czyli obliczenie skumulowanej sumy. Sygnał prostokątny może być obliczony z sygnału trójkątnego poprzez jego różniczkowanie, czyli obliczenie różnicy między wartościami. Do całkowania możemy użyć funkcji np.cumsum, a do różniczkowania – funkcji np.diff. Konieczne jest zastosowanie współczynnika skalującego amplitudę sygnału, tak aby zachować zmienność amplitudy w zakresie od -1 do 1.

In [16]:
b = 0.5
s = b * fs / f
trojkat2 = 2 * np.cumsum(prostokat) / s - 1
prostokat2 = np.diff((s / 2) * (trojkat2 + 1))
plt.plot(trojkat2)
plt.plot(prostokat2)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał trójkątny i prostokątny')
plt.show()

Moduł [scipy.signal] posiada funkcje do generowania opisanych wyżej sygnałów. Funkcja sawtooth generuje sygnał piłokształtny, a także trójkątny, jeżeli wartość drugiego parametru jest równa 0,5. Funkcja square generuje sygnał prostokątny, drugi parametr pozwala ustawić współczynnik kształtu. Pierwszym argumentem obu funkcji są wartości fazowe, które należy wygenerować identycznie jak dla funkcji sin (liniowo od 0 do $2\pi$).

In [17]:
f = 200
n = np.arange(1000)
faza = n * 2 * np.pi * f / fs
fig, ax = plt.subplots(4, sharex=True)
ax[0].plot(sig.sawtooth(faza))
ax[1].plot(sig.sawtooth(faza, 0.5))
ax[2].plot(sig.square(faza))
ax[3].plot(sig.square(faza, 0.8))
ax[3].set_xlabel('nr próbki')
for a in ax:
    a.set_ylabel('ampl.')
ax[0].set_title('Sygnały harmoniczne - scipy.signal')
plt.show()

Idealne sygnały harmoniczne mają nieskończone widmo. Konsekwencją tego jest powstawanie aliasingu widma przy cyfrowym generowaniu tych sygnałów w opisany wyżej sposób, o ile częstotliwość sygnału nie jest bardzo mała w porównaniu z częstotliwością Nyquista. Powstanie aliasingu powoduje, że odbite prążki pojawiają się w podstawowej części widma i albo sumują się ze składowymi w tej części (jeżeli stosunek częstotliwości Nyquista do częstotliwości sygnału jest liczbą wymierną), albo powodują powstanie nieharmonicznego sygnału. Poniżej pokazano widmo sygnału prostokątnego o częstotliwości podstawowej 1250 Hz. Aliasing jest bardzo wyraźnie widoczny.

In [60]:
prostokat2 = sig.square(2 * np.pi * np.arange(2048) * 1250 / fs)
w_prostokat = 20 * np.log10(np.abs(np.fft.rfft(prostokat2 * np.hamming(2048))) / 1024)
f_prostokat = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f_prostokat, w_prostokat)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.ylim(bottom=-50)
plt.title('Widmo sygnału prostokątnego (nieograniczone pasmo)')
plt.show()

W ogólnym przypadku nie można generować cyfrowych sygnałów harmonicznych opisanymi wyżej metodami, ponieważ z uwagi na występowanie aliasingu, sygnały w tej formie są nieprzydatne. Praktycznie stosowany może być jedynie sygnał trójkątny dla bardzo małych częstotliwości (np. 1 Hz przy częstotliwości próbkowania 48 kHz), ponieważ w tym przypadku aliasing jest pomijalny. Sygnał taki ma zastosowanie do modulacji innych sygnałów, np. w syntezatorze muzycznym (w module LFO).

Generowanie sygnałów o ograniczonym paśmie

Aby sygnały wygenerowane cyfrowo nie były zniekształcone przez aliasing, należy zadbać o to, by żadne ze składowych sygnału nie przekroczyły częstotliwości Nyquista. Sygnały spełniające ten warunek nazywa się sygnałami ograniczonymi pasmowo (bandlimited). W przypadku sygnałów harmonicznych, należy "przyciąć" widmo sygnału do częstotliwości Nyquista. Ponieważ sygnały harmoniczne są reprezentowane przez szereg harmoniczny, można je wygenerować ograniczając granice sumowania do żądanej częstotliwości. Metoda addytywna generowania sygnałów harmonicznych polega zatem na sumowaniu składowych sinusoidalnych o określonych amplitudach, w ograniczonym zakresie częstotliwości.

Na początek rozpatrzmy sygnał trójkątny. Szereg Fouriera dla tego sygnału ma następującą postać (stosujemy uproszczony zapis: $\sin f$ oznacza sygnał sinusoidalny o częstotliwości $f$), $A$ jest amplitudą sygnału.

$$y = \frac{8A}{\pi^2} \left( \sin f - \frac{1}{9} \sin 3f + \frac{1}{25} \sin 5f - \frac{1}{49} \sin 7f + ... \right)$$

W sygnale pojawiają się więc nieparzyste wielokrotności częstotliwości podstawowej, amplituda składowej o częstotliwości $kf$ jest proporcjonalna do $1/k^2$, ze znakiem na zmianę dodatnim i ujemnym.

Poniższa funkcja generuje sygnał trójkątny opisaną wyżej metodą. Wykresy pokazują wygenerowany przebieg czasowy dla częstotliwości 200 Hz oraz widmo dla częstotliwości 1250 Hz.

In [19]:
def gen_trojkat(n, f, fs):
    '''Generuje n próbek sygnału trójkątnego o częstotliwości f,
       ograniczonego pasmowo do fs/2.'''
    faza = np.arange(n) * 2 * np.pi * f / fs
    sygnal = np.zeros(n)
    k = 1
    znak = 1
    while k * f < fs / 2:
        sygnal += znak * (1 / k**2) * np.sin(k * faza)
        k += 2
        znak *= -1
    sygnal *= 8 / np.pi**2
    return sygnal
In [20]:
trojkat_op = gen_trojkat(1000, 200, fs)
plt.plot(trojkat_op)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał trójkątny ograniczony pasmowo, f=200 Hz')
plt.show()
In [62]:
w_trojkat_op = 20 * np.log10(np.abs(np.fft.rfft(gen_trojkat(2048, 1250, fs) * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), w_trojkat_op)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.ylim(bottom=-60)
plt.title('Widmo sygnału trójkątnego ograniczonego pasmowo, f=1250 Hz')
plt.show()

W przypadku sygnału trójkątnego o tej częstotliwości, postać czasowa sygnału nie odbiega od idealnego. Trzeba jednak pamiętać, że dla większych częstotliwości, liczba składowych maleje, a dla częstotliwości większej niż połowa częstotliwości Nyquista, wygenerowany sygnał jest sinusoidalny (wyższe składowe "nie mieszczą" się w widmie. Dla sygnału o czestotliwości 1250 Hz widać 10 składowych widmowych, aliasing oczywiście nie występuje.

W następnej kolejności rozpatrzmy sygnał piłokształtny. Szereg Fouriera ma w tym przypadku postać:

$$y = -\frac{2A}{\pi} \left( \sin f + \frac{1}{2} \sin 2f + \frac{1}{3} \sin 3f + \frac{1}{4} \sin 4f + ... \right)$$

Sygnał piłokształtny zawiera składowe o parzystych i nieparzystych wielokrotnościach częstotliwości podstawowej, amplituda $k$-tej składowej jest proporcjonalna do $1/k$.

In [22]:
def gen_pila(n, f, fs):
    '''Generuje n próbek sygnału piłokształtnego o częstotliwości f,
       ograniczonego pasmowo do fs/2.'''
    faza = np.arange(n) * 2 * np.pi * f / fs - np.pi
    sygnal = np.zeros(n)
    k = 1
    while k * f < fs / 2:
        sygnal += (1 / k) * np.sin(k * faza)
        k += 1
    sygnal *= -2 / np.pi
    return sygnal
In [23]:
pila_op = gen_pila(1000, 200, fs)
plt.plot(pila_op)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał piłokształtny ograniczony pasmowo, f=200 Hz')
plt.show()

W tym przypadku postać czasowa sygnału znacznie różni się od idealnej. Pamiętajmy, że w sygnale ograniczonym pasmowo brakuje składowych o wysokich częstotliwościach, a składowe te są konieczne do prawidłowego oddania fragmentów sygnału, w których następuje gwałtowna zmiana amplitudy (tutaj: "zęby" piły). Brak tych składowych objawia się w postaci widocznego na wykresie efektu "brzęczenia" (ringing), nazywanego efektem Gibbsa. Postać widma jest zgodna z oczekiwaniami.

In [63]:
w_pila_op = 20 * np.log10(np.abs(np.fft.rfft(gen_pila(2048, 1250, fs) * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), w_pila_op)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału piłokształtnego ograniczonego pasmowo, f=1250 Hz')
plt.ylim(bottom=-60)
plt.show()

Sygnał prostokątny o współczynniku kształtu 0,5 ("kwadratowy") jest przypadkiem pośrednim pomiędzy sygnałem trójkątnym a piłokształtnym. Posiada tylko składowe o nieparzystych wielokrotnościach częstotliwości podstawowej (jak trójkątny), ale ich amplituda jest proporcjonalna do $1/k$ (jak w piłokształtnym).

$$y = \frac{4}{\pi} \left( \sin f + \frac{1}{3} \sin 3f + \frac{1}{5} \sin 5f + \frac{1}{7} \sin 7f + ... \right)$$

Tym razem wygenerujemy przebieg w inny sposób: zamiast dodawać składowe sinusoidalne w pętli while, użyjemy rachunku macierzowego. Funkcja jest znacznie mniej czytelna, ale za to działa szybciej (każde obliczenia z użyciem pętli będą wolniejsze). Wektor k zawiera indeksy składowych, od 1 do $(f_s/2)/f$ włącznie. Wektor faza zawiera wartości kąta fazowego dla składowej podstawowej (o częstotliwości f). Za pomocą funkcji numpy.outer tworzymy macierz w wyniku mnożenia wektorów k i faza (iloczyn diadyczny), tak że poszczególne wiersze macierzy zawierają wartości kąta fazowego dla kolejnych składowych. Macierz tę mnożymy przez wektor amplitud (1/k), przy czym musimy go zamienić na macierz jednokolumnową za pomocą funkcji numpy.reshape (wartość -1 oznacza, że ten rozmiar zostanie obliczony automatycznie). Następnie sumujemy macierz wierszami za pomocą funkcji np.sum, parametr axis=0 oznacza sumowanie wzdłuż pierwszego wymiaru (wzdłuż kolumn).

In [25]:
def gen_kwadrat(n, f, fs):
    '''Generuje n próbek sygnału prostokątnego o współczynniku kształtu 0,5, o częstotliwości f,
       ograniczonego pasmowo do fs/2.'''
    k = np.arange(0, fs // (2 * f), 2) + 1
    faza = np.arange(n) * 2 * np.pi * f / fs
    amplitudy = np.reshape(1 / k, (-1, 1))
    sygnal = amplitudy * np.sin(np.outer(k, faza))
    sygnal = (4 / np.pi) * np.sum(sygnal, axis=0)
    return sygnal
In [26]:
kwadrat_op = gen_kwadrat(1000, 200, fs)
plt.plot(kwadrat_op)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał prostokątny 0,5 ograniczony pasmowo, f=200 Hz')
plt.show()

Tak jak w przypadku sygnału piłokształtnego, tak i tutaj brak składowych o wysokich częstotliwościach objawia się w postaci oscylacji w punktach nieciągłości sygnału.

In [64]:
w_kwadrat_op = 20 * np.log10(np.abs(np.fft.rfft(gen_kwadrat(2048, 1250, fs) * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), w_kwadrat_op)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału prostokątnego 0,5 ograniczonego pasmowo, f=1250 Hz')
plt.ylim(bottom=-60)
plt.show()

Sygnał prostokątny o współczynniku kształtu $b$ innym niż 0,5 jest najbardziej skomplikowany w generowaniu metodą addytywną. Szereg Fouriera jest w tym przypadku opisany sumą sygnałów kosinusoidalnych:

$$y = A(2b-1) + \frac{4}{\pi} \left(\sin(b\pi) \cos f + \frac{\sin(2b\pi)}{2} \cos 2f + \frac{\sin(3b\pi)}{3} \cos 3f + ... \right)$$

Sygnał zawiera zarówno parzyste, jak i nieparzyste wielokrotności częstotliwości podstawowej.

In [28]:
def gen_prostokat(n, f, fs, b):
    '''Generuje n próbek sygnału prostokątnego o współczynniku kształtu b, częstotliwości f,
       ograniczonego pasmowo do fs/2.'''
    faza = np.arange(n) * 2 * np.pi * f / fs
    sygnal = np.zeros(n)
    k = 1
    while k * f < fs / 2:
        sygnal += (np.sin(b * k * np.pi) / k) * np.cos(k * faza)
        k += 1
    sygnal = (4 / np.pi) * sygnal + 2 * b - 1
    return sygnal
In [29]:
prostokat_op = gen_prostokat(1000, 200, fs, 0.75)
plt.plot(prostokat_op)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał prostokątny 0,75 ograniczony pasmowo, f=200 Hz')
plt.show()
In [65]:
w_prostokat_op = 20 * np.log10(np.abs(np.fft.rfft(gen_prostokat(2048, 1250, fs, 0.8) * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), w_prostokat_op)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału prostokątnego 0,75 ograniczonego pasmowo, f=1250 Hz')
plt.ylim(bottom=-60)
plt.show()

Z powyższych przykładów wypływa następujący wniosek. Cyfrowe sygnały harmoniczne o postaci czasowej zgodnej z idealnym przebiegiem są zniekształcone przez aliasing widma. Ograniczenie pasma częstotliwości usuwa aliasing, ale zniekształca postać czasową sygnału. Nie jest możliwe wygenerowanie idealnego cyfrowego sygnału harmonicznego bez zniekształceń widmowych.

Generowanie sygnałów okresowych za pomocą FFT

Sygnały harmoniczne oraz inne sygnały okresowe można również generować poprzez zbudowanie żądanego widma i obliczenie odwrotnej transformaty Fouriera. Zasadnicza wada tego podejścia to brak możliwości uzyskania dowolnej częstotliwości sygnału, dostępne częstotliwości wynikają z rozmiaru transformaty i częstotliwości próbkowania. Uzyskanie dowolnej częstotliwości sygnału wymaga zastosowania interpolacji. Zaletą metody jest szybsze działanie niż w przypadku generowania addytywnego (zwłaszcza gdy trzeba obliczyć bardzo wiele składowych) i możliwość generowania sygnałów o dowolnym widmie.

Na początek rozpatrzmy sygnał sinusoidalny. Chcemy uzyskać $N$ próbek sygnału, zatem potrzebujemy $N$ wartości widma. Ponieważ tworzymy sygnał rzeczywisty, wystarczy określić jedną połowę widma i użyć funkcji numpy.fft.irfft. Dokładnie potrzebujemy $N/2+1$ wartości widma: pierwsza wartość określa składową stałą, ostatnia - częstotliwość Nyquista (dla uproszczenia zakładamy, że $N$ jest parzyste), pozostałe wartości określają amplitudę (i fazę) składowych sygnału. Wartość widma o indeksie $k$ reprezentuje składową sygnału o częstoliwości $kf_s/N$.

Tworzymy więc bufor o długości $N/2+1$ wypełniony zerami, pamietając że posługujemy się wartościami zespolonymi (np.complex). Amplituda sygnału ma wynosić 1. Ponieważ składowa od sygnału sinusoidalnego musi pojawić się w widmie dwa razy, dzielimy amplitudę na dwa. Składowa sinusoidalna stanowi część urojoną zespolonej pulsacji widmowej, dlatego wartość składowej widmowej musi być równa $-0,5j$. Wstawiając tę wartość pod indeks $k=1$ uzyskamy jeden okres sygnału, indeks $nk$ pozwoli uzyskać częstotliwość $n$ razy większą. Po obliczeniu IFFT i pominięciu szczątkowej części urojonej, musimy jeszcze przeskalować amplitudę sygnału ze współczynnikiem $N$ (funkcja irfft dzieli wynik przez rozmiar transformaty).

In [31]:
N = 2048
widmo_sinus = np.zeros(N // 2 + 1, dtype=np.complex)
widmo_sinus[1] = -0.5j
tab_sinus = N * np.real(np.fft.irfft(widmo_sinus))
plt.plot(tab_sinus)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Okres sygnału sinusoidalnego otrzymany z IFFT')
plt.show()

Otrzymaliśmy więc oczekiwany wynik.

W ten sam sposób możemy wygenerować jeden okres sygnału harmonicznego, posługując się zapisem tych sygnałów w formie szeregu Fouriera. Dla sygnału trójkątnego musimy obliczyć amplitudy nieparzystych składowych widmowych, zgodnie ze wzorem na szereg Fouriera (ponownie dzieląc amplitudy na dwa), pozstawiając parzyste składowe równe zeru. Należy pamiętać o nie nadpisywaniu składowych: stałej i dla częstotliwości Nyquista.

In [32]:
N = 2048
widmo_trojkat = np.zeros(N // 2 + 1, dtype=np.complex)
k = np.arange(1, len(widmo_trojkat), 2)
widmo_trojkat[k] = (-4 / np.pi**2) / k**2
tab_trojkat = N * np.real(np.fft.irfft(widmo_trojkat))
plt.plot(tab_trojkat)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Okres sygnału trójkątnego otrzymany z IFFT')
plt.show()

Dla sygnału piłokształtnego należy obliczyć amplitudy wszystkich składowych widma. Ustawienie dodatnich amplitud daje "piłę" narastającą, ujemnych – opadającą.

In [33]:
N = 2048
widmo_pila = np.zeros(N // 2 + 1, dtype=np.complex)
k = np.arange(1, len(widmo_pila))
widmo_pila[k] = (1j / np.pi) / k
tab_pila = N * np.real(np.fft.irfft(widmo_pila))
plt.plot(tab_pila)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Okres sygnału piłokształtnego otrzymany z IFFT')
plt.show()

Podobnie postępujemy dla sygnalu prostokątnego, regularnego (o współczynniku kształtu równym 0,5).

In [34]:
N = 2048
widmo_kwadrat = np.zeros(N // 2 + 1, dtype=np.complex)
k = np.arange(1, len(widmo_kwadrat), 2)
widmo_kwadrat[k] = (-2j / np.pi) / k
tab_kwadrat = N * np.real(np.fft.irfft(widmo_kwadrat))
plt.plot(tab_kwadrat)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Okres sygnału prostokątnego 0,5 otrzymany z IFFT')
plt.show()

Dla sygnału prostokątnego o współczynniku kształtu $b$ należy pamiętać o tym, że składniki szeregu Fouriera są reprezentowane przez kosinusy, więc musimy ustawić rzeczywiste wartości amplitud składowych widma. Ponadto konieczne jest przesunięcie wyniku na skali amplitudy, co uzyskujemy ustawiając składową stałą na wartość $2b-1$.

In [35]:
N = 2048
b = 0.8
widmo_prostokat = np.zeros(N // 2 + 1, dtype=np.complex)
k = np.arange(1, len(widmo_prostokat))
widmo_prostokat[k] = (2 / np.pi) * np.sin(k * b * np.pi) / k
widmo_prostokat[0] = 2 * b - 1
tab_prostokat = N * np.real(np.fft.irfft(widmo_prostokat))
plt.plot(tab_prostokat)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Okres sygnału prostokątnego 0,8 otrzymany z IFFT')
plt.show()

Generatory tablicowe

Załóżmy, że posiadając wartości próbek jednego okresu sygnału okresowego (wygenerowanego lub odczytanego np. z pliku), chcemy wygenerować sygnał ciągły o dowolnej częstotliwości. Układ lub program realizujący to podejście nazywa się generatorem tablicowym (wavetable). Mając okres sygnału zapisany w tablicy, możemy uzyskać sygnał ciągły odczytując ten okres w sposób cykliczny, a zmieniając krok odczytu, tak jak to opisano w rozdziale dotyczącym interpolacji, można uzyskać inne częstotliwości sygnału.

Jeżeli tablica zawiera okres sygnału sinusoidalnego, możemy swobodnie przeskalowywać sygnał metodą interpolacji bez ryzyka powstania aliasingu. Aby uzyskać dobrą dokładność, należy zadbać o dużą liczbę próbek opisujących okres. W takim przypadku wystarczające jest użycie prostej metody interpolacji, np. liniowej.

W przypadku gdy składowe widmowe sygnału wypełniają całe pasmo częstotliwości, jak np. dla sygnału piłokształtnego wygenerowanego opisywaną wcześniej metodą, zmiana częstotliwości sygnału przez interpolację spowoduje aliasing. Dlatego należy ograniczyć jeszcze bardziej pasmo sygnału, zostawiając "rezerwę" dla interpolacji. Jeżeli np. założymy, że chcemy zwiększyć częstotliwość sygnału maksymalnie o oktawę, co oznacza podwojenie częstotliwości, musimy ograniczyć pasmo sygnału do $f_s/4$, tak że po dwukrotnym zwiększeniu częstotliwości, sygnał wypełni pasmo do $f_s/2$. Jeżeli generujemy sygnał samodzielnie, musimy zadbać o brak składowych powyżej $f_s/4$. Jeżeli mamy istniejący sygnał, musimy go przetworzyć filtrem dolnoprzepustowym bądź "obciąć" pasmo w widmie sygnału.

Dla przykladu, wygenerujmy tablicę dla sygnału piłokształtnego, zawierającą składowe do $f_s/4$, metodą FFT.

In [36]:
widmo_pila_2048 = np.zeros(1025, dtype=np.complex)
k = np.arange(1, 2048 // 4)
widmo_pila_2048[k] = (1j / np.pi) / k
tab_pila_2048 = 2048 * np.real(np.fft.irfft(widmo_pila_2048))

Odczytując cyklicznie próbki z tablicy z częsttotliwością 48 kHz, otrzymamy częstotliwość równą 48000 / 2048 = 23,44 Hz. Za pomocą interpolacji możemy uzyskać częstotliwości do wartości dwukrotnie większej (46,87 Hz), bez ryzyka aliasingu. Wygenerujmy tą metodą sygnał o częstotliwości 44 Hz.

In [67]:
N = 2048
k = 44 / (fs / N)
indeksy = (np.arange(N) * k) % N
sygnal = np.interp(indeksy, np.arange(2048), tab_pila_2048, period=N)
widmo = 20 * np.log10(np.abs(np.fft.rfft(sygnal * np.hamming(2048))) / 1024)
plt.plot(sygnal)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał piłokształtny interpolowany, f=44 Hz')
plt.figure()
plt.plot(np.fft.rfftfreq(2048, 1 / fs), widmo)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału piłokształtnego interpolowanego, f=44 Hz')
plt.ylim(bottom=-70)
plt.show()

Co zrobić jeżeli chcemy uzyskać wyższą częstotliwość, np. 90 Hz? Nie możemy dalej interpolować tej samej tablicy powyżej częstotliwości 46,88 Hz, ponieważ powstanie aliasing. Potrzebujemy nowej tablicy dla kolejnej okawy. Długość tablicy powinna być dwukrotnie mniejsza niż poprzednia, czyli 1024. Możemy wygenerować nową tablicę tak jak w poprzednim przypadku, możemy też obliczyć widmo tablicy, "przyciąć" je do rozmiaru 513, wyzerować składowe powyżej $f_s/4$ i obliczyć transformatę odwrotną.

In [68]:
N = 1024
widmo_pila_1024 = widmo_pila_2048[:513]
widmo_pila_1024[256:] = 0
tab_pila_1024 = 1024 * np.real(np.fft.irfft(widmo_pila_1024)) 
k = 90 / (fs / N)
indeksy = (np.arange(2048) * k) % N
sygnal = np.interp(indeksy, np.arange(N), tab_pila_1024, period=N)
widmo = 20 * np.log10(np.abs(np.fft.rfft(sygnal * np.hamming(2048))) / 1024)
plt.plot(sygnal)
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał piłokształtny interpolowany, f=90 Hz')
plt.figure()
plt.plot(np.fft.rfftfreq(2048, 1 / fs), widmo)
plt.xlabel('częstotliwość (Hz)')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału piłokształtnego interpolowanego, f=90 Hz')
plt.ylim(bottom=-70)
plt.show()

W ten sam sposób możemy generować kolejne poziomy tablic o zmniejszającej się długości (512, 256, ...) dla coraz wyższych częstotliwości. Generator tablicowy w opisanej tu formie składa się więc z kilku tablic (poziomów) okresu sygnału o różnych rozmiarch, generujących częstotliwości różniące się o oktawę. Generowanie sygnału polega na wyborze największej tablicy z tych, które generują częstotliwości mniejsze niż pożądana, oraz na odczycie próbek z tej tablicy w sposób cykliczny, z interpolacją.

Wadą opisywanego rozwiązania jest to, że w przypadku płynnych zmian częstotliwości odczytywanego sygnału dochodzi do "przeskoków" pomiędzy różnymi tablicami. Ponieważ sąsiednie tablice różnią się energią sygnału – jedna (po interpolacji) ma prawie pełne pasmo, druga (bez interpolacji) tylko połowę pasma – może to powodować słyszalne wahania głośności sygnału. W praktyce rozwiązuje się ten problem mieszając dwa sygnały, interpolowane "w górę" i "w dół" z sąsiednich tablic, co powoduje wyrównywanie różnic głośności. Inną metodą jest dopuszczenie do powstawania aliasingu w pewnym zakresie, np. powyżej 16 kHz, zakładając, że nie będzie on słyszalny dla większości osób.

Odtwarzanie próbek dźwięku – sampling

Samplery (sampler) to elektroniczne instrumenty muzyczne, których działanie polega na odtwarzaniu wgranych do niego brzmień – próbek dźwiękowych (samples), ze zmianą wysokości dźwięku i zapętleniem sygnału. Umożliwiają one stworzenie instrumentu muzycznego z dowolnego brzmienia, które można nagrać. Działanie samplera jest podobne do przedstawionych wcześniej generatorów tablicowych, z tą różnicą, że sygnały nie są generowane w instrumencie, a wprowadzane z zewnątrz. Próbki dźwiękowe bardzo często składają się z dwóch części.

  • Transjent (transient) – jest to początkowy fragment dźwięku, który jest odtwarzany jeden raz, od początku do końca. Może to być np. początkowy fragment dźwięku instrumentu muzycznego (tzw. faza ataku), w którym dźwięk jest budowany i która ma bardzo duże znaczenie dla charakteru dźwięku.
  • Pętla (loop) jest odtwarzania "w kółko", np. aż do zwolnienia klawisza na klawiaturze muzycznej. Stosuje się ją do powtarzalnego fragmentu dźwięku. Np. w dźwięku instrumentu muzycznego można znaleźć "pseudookres" sygnału i go zapętlić.

Niektóre próbki dźwięku mogą nie posiadać pętli, są wtedy odtwarzane jeden raz, od początku do końca. Możliwa jest też próbka bez transjentu.

Odtworzenie dźwięku w samplerze polega na "wygenerowaniu" próbek transjentu (o ile występuje), a następnie na cyklicznym odczytywaniu próbek pętli. Odbywa się to w taki sam sposób, jak w przypadku generatorów tablicowych. Po wyłączeniu dźwięku, pętla jest nadal odtwarzana, ze stopniowym wygaszaniem aplitudy. Zmiana wysokości dźwięku (transpozycja) odbywa się za pomocą interpolacji. Należy tu zwrócić szczególną uwagę na ryzyko powstania aliasingu i zapobiec mu np. przez filtrację sygnału.

Ważną rzeczą, o której należy wspomnieć, są zniekształcenia czasowe transjentu przy zmianie wysokości. Jak omówiono w rozdziale dotyczącym interpolacji, zwiększanie kroku odczytu próbek z pamięci powoduje zwiększanie częstotliwości, ale i skracanie dźwięku. Podobnie, zmniejszanie kroku to zmniejszenie częstotliwości i wydłużenie dźwięku. Duże zmiany czasu trwania transjentu mogą spowodować, że dźwięk straci swój charakter. Tak jak w generatorze tablicowym stosuje się kilka poziomów tablic dla różnych częstotliwości, tak i w samplerze stosuje się multisampling - nagrywa się wiele próbek tego samego brzmienia o różnych wysokościach, tak aby skrócić zakresy przeskalowywania próbek.

Poniżej pokazano początkowy fragment nagrania dźwięku klarnetu. W pierwszej fazie dźwięk jest tworzony, jego amplituda narasta Po pewnym czasie kształt dźwięku się stabilizuje, zmiany amplitudy wynikają ze sposobu grania na instrumencie. Na tej podstawie wyodrębniono transjent i pętlę, zaznaczone pionowymi liniami. Kolejny wykres pokazuje zawartość pętli – jeden pseudookres sygnału.

In [39]:
fs_klarnet, wav_klarnet = wavfile.read('pliki/klarnet.wav')
wav_klarnet = wav_klarnet / np.max(np.abs(wav_klarnet))
petla_poczatek = 2103
petla_koniec = 2202
plt.figure(figsize=(12, 4))
plt.plot(np.arange(4410) / fs_klarnet, wav_klarnet[:4410])
plt.axvline(petla_poczatek / fs_klarnet, c='red', lw=1)
plt.axvline(petla_koniec / fs_klarnet, c='olive', lw=1)
plt.xlabel('czas [s]')
plt.ylabel('amplituda')
plt.title('Próbka dźwięku klarnetu')
plt.show()
In [40]:
transjent = wav_klarnet[:petla_poczatek]
petla = wav_klarnet[petla_poczatek:petla_koniec]
plt.plot(petla)
plt.xlabel('nr probki')
plt.ylabel('amplituda')
plt.title('Próbka dźwięku klarnetu - pętla')
plt.show()

Częstotliwość podstawowa dźwięku, jak zmierzono w jednym z poprzednich rozdziałów, jest równa ok. 441 Hz. Wygenerujemy teraz sygnał dla częstotliwości 500 Hz, z zapętleniem sygnału i wygaszeniem na końcu.

In [41]:
k = 500 / 441
# Transjent
indeks_transjent = np.arange(0, len(transjent), k)
sygnal_transjent = np.interp(indeks_transjent, np.arange(len(transjent)), transjent)
# Pętla
indeks_petla = np.arange(k * 2000) % len(petla)
sygnal_petla = np.interp(indeks_petla, np.arange(len(petla)), petla, period=len(petla))
# Cały sygnał
sygnal = np.hstack((sygnal_transjent, sygnal_petla))
# Wygaszenie dźwięku
wygaszanie = np.linspace(1.0, 0.0, 500)
sygnal[-500:] *= wygaszanie

plt.figure(figsize=(12, 4))
plt.plot(np.arange(len(sygnal)) / fs_klarnet, sygnal)
plt.xlabel('czas [s]')
plt.ylabel('amplituda')
plt.title('Odtworzona próbka dźwięku klarnetu, f=500 Hz')
plt.show()

Sygnał został odtworzony, ze zmianą czestotliwości podstawowej. Nie zostały zachowane zmiany obwiedni amplitudowej, obecne w oryginalnym dźwięku. Można je wprowadzić do sygnału albo za pomocą sterownika, np. pokrętła na klawiaturze muzycznej, sterującego głośnością dźwięku, albo stosując sygnał modulujący.

Generowanie sygnału świergotowego

Sygnał świergotowy (chirp), nazywany również sygnałem przemiatającym (sweep), jest to sygnał o zmiennej częstotliwości chwilowej. Najprostszym przypadkiem sygnału świergotowego jest sygnał sinusoidalny, którego częstotliwość chwilowa jest przestrajana w trakcie generowania. Najczęściej stosuje się liniowe lub logarytmiczne przestrajanie częstotliwości. Sygnały świergotowe stosuje się np. do badania charakterystyki częstotliwościowej układów, podając na wejście sygnał świergotowy, rejestruje się odpowiedź systemu na pobudznia o różnych czestotliwościach.

Sygnał świergotowy przestrajany liniowo (linear chirp) możemy wygenerować w taki sam sposób, jak sygnał sinusoidalny, ale stosując inną wartość częstotliwości dla każdej próbki. Poniższy przykład generuje sygnał o czasie trwania 100 ms, przestrajany liniowo w zakresie od 1 Hz do 200 Hz. Postać widmową sygnału najwygodniej jest przedstawić za pomocą spektrogramu.

In [42]:
f = np.linspace(1, 200, fs // 10)
chirp = np.sin(2 * np.pi * np.arange(fs // 10) * f / fs)
plt.plot(np.arange(fs // 10) / fs, chirp)
plt.xlabel('czas [s]')
plt.ylabel('amplituda')
plt.title('Sygnał świergotowy, liniowy 1-200 Hz')
plt.show()
In [43]:
ff, tt, Sxx = sig.spectrogram(chirp, fs=fs, noverlap=384, nperseg=512, nfft=2048)
plt.pcolormesh(tt, ff, Sxx, shading='gouraud')
plt.ylim(0, 500)
plt.xlabel('czas [s]')
plt.ylabel('częstotliwość [Hz]')
plt.title('Spektrogram sygnału świergotowego, liniowego')
plt.show()

W module scipy.signal mamy funkcję chirp, która generuje sygnał świergotowy. Jej argumenty to kolejno: punkty obliczania sygnału, częstotliwość f0 dla pierwszej próbki, punkt t1, częstotliwość f1 dla punktu t1 oraz metoda przestrajania częstotliwości: linear (liniowa), quadratic (kwadratowa), logarithmic (logarytmiczna) lub hyperbolic (hiperboliczna). Częstotliwość chwilowa jest obliczana przez dopasowanie danej funkcji do punktów (0, f0) i (t1, f1) Poniżej przykład dla przestrajania logarytmicznego od 20 Hz do 500 Hz w ciągu 1 sekundy.

In [44]:
t = np.linspace(0, 1, fs)
chirpl = sig.chirp(t, 20, t[-1], 500, 'logarithmic')
plt.plot(t, chirpl)
plt.xlabel('czas [s]')
plt.ylabel('amplituda')
plt.title('Sygnał świergotowy, logarytmiczny 20-500 Hz')
plt.figure()
ff, tt, Sxx = sig.spectrogram(chirpl, fs=fs, noverlap=256, nperseg=512, nfft=2048)
plt.pcolormesh(tt, ff, Sxx, shading='gouraud')
plt.ylim(0, 500)
plt.xlabel('czas [s]')
plt.ylabel('częstotliwość [Hz]')
plt.title('Spektrogram sygnału świergotowego, logarytmicznego')
plt.show()

Generowanie sygnałów losowych

Sygnałami losowymi nazywamy takie, w których przyjmowanie różnych wartości jest opisane funkcją rozkładu prawdopodobieństwa. W wielu dziedzinach, sygnał losowy nazywa się szumem (noise). Generowanie sygnału losowego jest przydatne np. w testowaniu układów w szerokim zakresie częstotliwości, lub w syntezie dźwięków instrumentów perkusyjnych.

Jednym z najczęściej stosowanych sygnałów losowych jest sygnał o równomiernym rozkładzie prawdpodobieństwa (uniform noise). Sygnał taki przyjmuje wartości z ustalonego zakresu (np. od 0 do 1) z jednakowym prawdopodobieństwem. Do wygenerowania takiego sygnału można użyć funkcji numpy.random.random. Zwracane są wartości z przedziału półotwartego [0, 1). Aby uzyskać sygnał losowy, musimy te wartości przemnożyć przez 2 i odjąć 1, uzyskując w ten sposób zakres wartości [-1, 1).

In [45]:
szumb = 2 * np.random.random(2048) - 1
plt.plot(szumb[:100])
plt.xlabel('nr próbki')
plt.ylabel('amplituda')
plt.title('Sygnał szumowy o rozkładzie równomiernym')
plt.show()

Wartości sygnału zmieniają się w sposób losowy. Widmo amplitudowe szumu o rozkładzie równomiernym ma następującą postać.

In [69]:
widmo_szumb = 20 * np.log10(np.abs(np.fft.rfft(szumb)) / (len(szumb) / 2))
plt.plot(np.fft.rfftfreq(len(szumb), 1 / fs), widmo_szumb)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo amplitudowe szumu o rozkładzie równomiernym')
plt.show()

Widmo sygnału losowego o rozkładzie równomiernym jest w przybliżeniu płaskie. Tutaj obliczyliśmy widmo dla 2048 próbek – jest to zbyt krótki odcinek czasowy, aby poprawnie uchwycić losowość sygnału. Powtarzając obliczenia widma dla wielu odcinków i uśredniając wyniki, otrzymalibyśmy w przybliżeniu płaskie widmo. Charakterystyka fazowa tego sygnału ma charakter przypadkowy. Opisywany sygnał losowy często nazywa się szumem białym (white noise), przez analogię do światła białego, które również ma płaską charakterystykę widmową.

Ważną cechą sygnału losowego jest kształt jego funkcji autokorelacji.

In [47]:
autokor = sig.correlate(szumb, szumb) / np.sum(szumb**2)
plt.plot(np.arange(-len(szumb) + 1, len(szumb)), autokor)
plt.xlabel('przesunięcie [próbki]')
plt.ylabel('autokorelacja')
plt.title('Autokorelacja szumu białego')
plt.show()

Dla zerowego przesunięcia, funkcja autokorelacji zawsze przyjmuje wartość 1. Przy dowolnym przesunięciu, kopie sygnału szumowego powinny być całkowicie nieskorelowane ze sobą, z uwagi na losowy charakter sygnału, dlatego funkcja autokorelacji powinna być bliska zeru. Obliczenie autokorelacji jest więc swego rodzaju testem jakości sygnału losowego. Dla idealnego sygnału losowego, funkcja autokorelacji powinna być równa impulsowi jednostkowemu (delcie Kroneckera).

Sygnały losowe generowane przez algorytmy komputerowe są w rzeczywistości sygnałami pseudolosowymi. Są one obliczane za pomocą deterministycznego algorytmu, jednak wytwarzane wartości powinny tworzyć oczekiwany rozkład prawdopodobieństwa. Najczęściej stosowanym algorytmem generowania sygnału losowego o rozkładzie równomiernym jest liniowy generator kongruentny (LCG, linear congruential generator). Oblicza on wartości sygnału za pomocą prostej zależności:

$$x(n) = \left(ax(n-1)+c\right)\bmod m$$

Parametry generatora to: mnożnik $a$, inkrement $c$, wartość modulo $m$ oraz wartość początkowa $x(0)$, nazywana ziarnem (seed). Dla ustalonych wartości parametrów, algorytm wygeneruje zawsze taką samą sekwencję wartości. Aby uzyskiwać różne wyniki dla kolejnych wywołań generatora, zazwyczaj ustawia się ziarno na wartość ulegającą zmianie, np. z zegara systemowego. W języku Python można do tego użyć funkcji time z modułu time.

Poniżej pokazano generator wytwarzający wartości sygnału losowego o rozkładzie równomiernym metodą LCG, stosując wartości parametrów zgodnie ze standardem ANSI C. Funkcja ta działa iteracyjnie, więc w języku Python będzie bardzo wolna (jest tu pokazana w celach demonstracyjnych), ale ten sam algorytm może być zaimplementowany w kompilowanym kodzie lub w Asemblerze).

In [48]:
def generuj_szum(n, ziarno=None):
    m = 2**31
    a = 1103515245
    c = 12345
    if ziarno is None:
        ziarno = time.time() * 1e7
    ziarno = np.uint32(ziarno % m)
    y = np.empty(n, dtype=np.uint32)
    x = ziarno
    for i in range(n):
        x = (a * x + c) % m
        y[i] = x
    return y / 2**31

Przetestujemy generator tworząc po sześć liczb losowych. Jeżeli podamy tę samą wartość ziarna, generator wytworzy tę samą sekwencję liczb – widać, że jest to algorytm deterministyczny. Jeżeli nie podamy wartości ziarna, zostanie ono obliczone na podstawie aktualnego czasu. Za pomocą funkcji time.sleep wstrzymujemy działanie programu na 1 ms aby otrzymać inną wartość ziarna.

In [49]:
print('ziarno = 121212')
print(generuj_szum(6, 121212))
print(generuj_szum(6, 121212))
print(generuj_szum(6, 121212))
print('ziarno time()')
print(generuj_szum(6))
time.sleep(0.001)
print(generuj_szum(6))
time.sleep(0.001)
print(generuj_szum(6))

ziarno = 121212
[ 0.52311921  0.0774806   0.41863972  0.268099    0.36139792  0.99526336]
[ 0.52311921  0.0774806   0.41863972  0.268099    0.36139792  0.99526336]
[ 0.52311921  0.0774806   0.41863972  0.268099    0.36139792  0.99526336]
ziarno time()
[ 0.75705795  0.43901506  0.28192975  0.70836727  0.0756818   0.58025571]
[ 0.41342594  0.05272984  0.79266878  0.79592538  0.49850054  0.98540065]
[ 0.79250734  0.25469953  0.22666395  0.69201843  0.10413045  0.2873085 ]

Obliczenie funkcji autokorelacji udowadnia, że wygenerowana sekwencja próbek ma charakter sygnału losowego.

In [50]:
szumb2 = 2 * generuj_szum(2048) - 1
autokor2 = sig.correlate(szumb2, szumb2) / np.sum(szumb2**2)
plt.plot(np.arange(-len(szumb2) + 1, len(szumb2)), autokor2)
plt.xlabel('przesunięcie [próbki]')
plt.ylabel('autokorelacja')
plt.title('Autokorelacja obliczonego sygnału losowego')
plt.show()

Oprócz szumu białego, stosuje się inne sygnały losowe, np. sygnały o odpowiednio ukształtowanym widmie – szum różowy, szum brązowy – a także szumy wąskopasmowe, które można uzyskać filtrując szum biały za pomocą filtrów środkowoprzepustowych.

Modulacja sygnałów

Modulacją sygnału (modulation) nazywamy zmianę parametrów sygnału za pomocą innego sygnału. W typowych rozwiązaniach, sygnałem poddawanym modulacji, nazywanym sygnałem nośnym lub nośną (carrier), jest sygnał sinusoidalny. W takim przypadku mamy trzy parametry, które możemy modulować: amplitudę, częstotliwość i fazę.

Modulacja sygnałów analogowych jest szeroko stosowana w telekomunikacji, np. do transmisji fal radiowych. W technice cyfrowej modulacja ma mniejsze znaczenie, z uwagi na dyskretny charakter sygnałów i powstawanie aliasingu. Cyfrowa modulacja jest spotykana np. w cyfrowej syntezie dźwięku. Teoria modulacji sygnałów stanowi bardzo obszerny temat, dlatego tutaj będą przedstawione tylko bardzo podstawowe zagadnienia. Pominięty zostanie m.in. aspekt demodulacji sygnałów.

Modulacja amplitudy (AM, amplitude modulation) polega na modulowaniu amplitudy sygnału nośnego przez sygnał modulujący:

$$y(n) = (A_c + x_m(n)) \sin(2\pi n f_c / f_s) = (A_c + A_m \sin(2\pi n f_m/f_s)) \sin(2\pi n f_c / f_s)$$

Przez $f_m$ oznaczamy częstotliwość modulującą, $f_c$ – częstotliwość nośną, $A_m$ – amplitudę sygnału modulującego, $A_c$ – amplitudę sygnału nośnego. Stosunek $A_m/A_c$ nazywa się głębokością modulacji.

W poniższym przykładzie, sygnał nośny o częstotliwości 5 kHz i amplitudzie 1 jest modulowany sygnałem sinusoidalnym o częstotliwości 500 Hz i amplitudzie 0,25.

In [51]:
nosna = np.sin(2 * np.pi * np.arange(2048) * 5000 / fs)
modulator = np.sin(2 * np.pi * np.arange(2048) * 500 / fs)

mod_am = (1 + 0.25 * modulator) * nosna
plt.plot(np.arange(500) / fs, mod_am[:500])
plt.xlabel('czas [s]')
plt.ylabel('amplituda')
plt.title('Sygnał z modulacją amplitudy')
plt.show()

Modulacja amplitudy wpływa zatem na obwiednię sygnału. Widmo sygnału zmodulowanego amplitudowo wygląda następująco. Występuje składowa na częstotliwości nośnej. Składowe widmowe sygnału modulującego znajdują się po obu stronach składowej nośnej – są one nazywane wstęgą dolną i wstęgą górną. W tym przykładzie, wstęga dolna to składowa na częstotliwości 5000 - 500 = 4500 Hz, wstęga górna – na częstotliwości 5000 + 500 = 5500 Hz. W taki sam sposób powstaje widmo w przypadku gdy sygnał modulujący zawiera wiele składowych widmowych, musi to jednak być sygnał o ograniczonym paśmie.

In [77]:
widmo_am = 20 * np.log10(np.abs(np.fft.rfft(mod_am * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), widmo_am)
plt.xlim(0, 8000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału z modulacją amplitudy')
plt.ylim(bottom=-60)
plt.show()

Analogowa modulacja AM jest nadal używana do transmisji sygnału radiowego na falach średnich i długich, chociaż większość rozgłośni transmituje sygnał na falach ultrakrótkich. W przypadku fal długich, sygnał nośny o częstotliwości z zakresu 30-300 kHz jest modulowany amplitudowo sygnałem radiowym o paśmie czestotliwości ograniczonym do zakresu słyszalnego.

Modulacja częstotliwości (FM, frquency modulation) polega na modulowaniu chwilowej częstotliwości sygnału.

$$y(n) = A_c \sin\left(2\pi n (f_c+x_m(n))/f_s\right) = A_c \sin\left(2\pi n (f_c+ A_m \sin(2\pi n f_m/f_s))/f_s\right)$$

Możemy ją uzyskać zmieniając wartość częstotliwości używaną do obliczenia kąta fazowego. Poniższym przykład pokazuje modulację częstotliwości dla sygnału nośnego o częstotliwości 500 Hz i amplitudzie 1, za pomocą sygnału sinusoidalnego o częstotliwości 100 Hz i amplitudzie 30.

In [53]:
modulator = np.sin(2 * np.pi * np.arange(2048) * 100 / fs)
mod_fm = np.sin(2 * np.pi * np.arange(2048) * (500 + 30 * modulator) / fs)
plt.plot(np.arange(2048) / fs, mod_fm)
plt.xlabel('czas [s]')
plt.ylabel('amplituda')
plt.title('Sygnał z modulacją częstotliwości')
plt.show()

Widmo sygnału zmodulowanego częstotliwościowo pokazano dla sygnału nośnego o częstotliwości 5000 Hz i amplitudzie 1, zmodulowanego za pomocą sygnału sinusoidalnego o częstotliwości 500 Hz i amplitudzie 10. W tym przypadku, wstęga dolna i górna zawierają kilka składowych. Szerokość pasma zależy od indeksu modulacji, czyli stosunku amplitudy sygnału modulującego do częstotliwości sygnału nośnego.

In [76]:
modulator = np.sin(2 * np.pi * np.arange(2048) * 500 / fs)
mod_fm = np.sin(2 * np.pi * np.arange(2048) * (5000 + 10 * modulator) / fs)
widmo_fm = 20 * np.log10(np.abs(np.fft.rfft(mod_fm * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), widmo_fm)
plt.xlim(0, 8000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału z modulacją częstotliwości')
plt.ylim(bottom=-60)
plt.show()

Analogowa modulacja częstotliwości jest stosowana m.in do transmisji sygnału radiowego za pomocą fal ultrakrótkich (UKF). W Polsce, częstotliwości nośne, przydzielane rozgłośniom radiowym mieszczą się w zakresie 87,5 – 108 MHz.

Modulacja fazy (PM, phase modulation) jest podobna do modulacji częstotliwości. Różnica polega na tym, że wartości sygnału modulującego są dodawane do kąta fazowego sygnału nośnego.

$$y(n) = A_c \sin\left(2\pi n f_c/f_s + x_m(n)\right) = A_c \sin\left(2\pi n f_c/f_s + A_m \sin(2\pi n f_m/f_s))\right)$$

In [55]:
modulator = np.sin(2 * np.pi * np.arange(2048) * 100 / fs)
mod_pm = np.sin(2 * np.pi * np.arange(2048) * 500 / fs + 10 * modulator)
plt.plot(np.arange(2048) / fs, mod_pm)
plt.xlabel('czas [s]')
plt.ylabel('amplituda')
plt.title('Sygnał z modulacją fazy')
plt.show()
In [75]:
modulator = np.sin(2 * np.pi * np.arange(2048) * 500 / fs)
mod_pm = np.sin(2 * np.pi * np.arange(2048) * 5000 / fs + 1 * modulator)
widmo_pm = 20 * np.log10(np.abs(np.fft.rfft(mod_pm * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), widmo_pm)
plt.xlim(0, 8000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału z modulacją fazy')
plt.ylim(bottom=-60)
plt.show()

W stosunku do modulacji częstotliwości, amplituda sygnału modulującego w przypadku modulacji fazy musi być mniejsza.

Interesujący przypadek zachodzi gdy obie częstotliwości, modulująca i nośna, są identyczne lub ich stosunek jest równy stosunkowi małych liczb naturalnych. Wówczas powstaje sygnał harmoniczny. Poniżej przypadek dla obu częstotliwości równych 440 Hz i amplitudy sygnału modulującego równej 10.

In [74]:
modulator = np.sin(2 * np.pi * np.arange(2048) * 440 / fs)
mod_pm = np.sin(2 * np.pi * np.arange(2048) * 440 / fs + 10 * modulator)
widmo_pm = 20 * np.log10(np.abs(np.fft.rfft(mod_pm * np.hamming(2048))) / 1024)
plt.plot(np.fft.rfftfreq(2048, 1 / fs), widmo_pm)
plt.xlim(0, 8000)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Widmo sygnału z modulacją fazy, $f_c$=$f_m$=440 Hz')
plt.ylim(bottom=-60)
plt.show()

Metoda ta jest podstawą syntezy dźwięku metodą modulacji częstotliwości (FM synthesis) – nazwa metody została nieprawidłowo dobrana, gdyż wykorzystywana jest modulacja fazy (w technice cyfrowej jest ona łatwiejsza do implementacji niż modulacja częstotliwości). Metoda syntezy FM była wykorzystywana m.in w instrumentach firmy Yamaha produkowanych w latach 80. ubiegłego stulecia, najbardziej popularnym instrumentem był Yamaha DX7. Stosowano tam wielokrotną modulację fazy w układzie sześciu generatorów sygnału sinusoidalnego, łączonych w różne kombinacje.

Podsumowanie

  • Sygnał sinusoidalny generujemy obliczając wartości pulsacji (kąta fazowego) i stosując funkcję sin dla każdej wartości.
  • Jeżeli nie dysponujemy funkcją obliczającą wartości sinusów, możemy je aproksymować za pomocą szeregu Taylora. Można też wykorzystać filtr IIR na granicy stabilności jako oscylator.
  • Najczęściej stosowane sygnały harmonoczne to sygnał trójkątny, piłokształtny i prostokątny. Cyfrowe generowanie tych sygnałów "z definicji" prowadzi do powstania aliasingu, zniekształcającego sygnał.
  • Sygnały harmoniczne o ograniczonym paśmie można uzyskać poprzez sumowanie składowych sinusodialnych o czestotliwościach nie przekraczających częstotliwości Nyquista. Generując sygnały tą metodą, zniekształcamy ich postać czasową.
  • Można generować sygnały okresowe tworząc ich widmo i wykonując odwrotne przekształcenie Fouriera. Nie możemy w ten sposób uzyskac dowolnej częstotliwości sygnału.
  • Mając wygenerowaną pewną liczbę próbek tworzących okres sygnału, można wygenerować sygnał o dowolnej częstotliwości za pomoca interpolacji. Należy tak ograniczyć widmo sygnału, aby nie powstawal aliasing.
  • Generatory tablicowe składają się zwykle z kilku tablic okresu sygnału o różnych długościach, umożliwiających interpolację w określonym zakresie bez aliasingu.
  • W samplerach stosuje się wytwarzanie sygnału przez odczytanie z pamięci próbek dźwięku. Transjent, pierwsza część próbki, jest odtwarzany raz, interpolacja powoduje zniekształcenia czasowe. Pętla, druga część próbki, jest odtwarzana cyklicznie.
  • Sygnał świergotowy możemy wygenerować zmieniając chwilową częstotliwość sygnału sinusoidalnego, w sposób np. liniowy lub logarytmiczny.
  • Szum biały można wygenerować stosując generator liczb pseudolosowych o rozkładzie równomiernym. Jeżeli nie mamy takiego generatora, możemy samodzielnie stworzyć liniowy generator kongruentny.
  • W sinusoidalnym sygnale nośnym możemy modulować amplitudę, częstotliwość i fazę za pomocą sygnału modulującego. Analogowa modulacja jest stosowana np. do transmisji fal radiowych, w technice cyfrowej modulacja tego typu jest rzadko stosowana.