Základní funkce pro zpracování signálů v knihovně SciPy

Ve třetí části seriálu o knihovně SciPy se zaměříme na dvě oblasti, které spolu částečně souvisejí. Nejprve se zmíníme o datových sadách, které knihovna SciPy nabízí a které následně použijeme v demonstračních příkladech. Ve druhé části článku si ukážeme některé funkce používané pro zpracování signálů, a to jak signálů jednorozměrných, tak i vícerozměrných. Pro ukázku možností, které SciPy v této oblasti nabízí, bude použita jedna předpřipravená datová sada obsahující delší záznam EKG.

Ještě před spuštěním demonstračních příkladů je nutné do projektového souboru přidat další balíček, jenž je nazvaný pooch.. Jedná se o balíček, který SciPy volá při stahování datových sad (viz též https://pypi.org/project/pooch/):

$ uv add pooch

Projektový soubor pyproject.toml může vypadat takto:

[project]
name = "scipy-lib"
version = "0.1.0"
description = "Scipy examples"
readme = "README.md"
requires-python = "==3.11.*"
dependencies = [
    "matplotlib>=3.10.8",
    "pooch>=1.9.0",
    "scipy>=1.17.1",
]

Balíček scipy.datasets

Datové sady, které je možné použít pro vyzkoušení tříd a funkcí poskytovaných knihovnou SciPy, jsou součástí balíčku nazvaného scipy.datasets. Ve skutečnosti tento balíček neobsahuje přímo všechna potřebná data (což jsou, jak uvidíme dále, n-rozměrná pole knihovny NumPy), ovšem umožňuje je v případě potřeby stáhnout (což ovšem způsobuje problémy například při sestavování aplikací v prostředí odpojeném od Internetu):

import scipy.datasets as datasets
 
help(datasets)

Nápověda poskytovaná k tomuto balíčku:

Help on package scipy.datasets in scipy:
 
NAME
    scipy.datasets
 
DESCRIPTION
    ================================
    Datasets (:mod:`scipy.datasets`)
    ================================
 
    .. currentmodule:: scipy.datasets
 
    Dataset Methods
    ===============
 
    .. autosummary::
       :toctree: generated/

       ascent
       face
       electrocardiogram
 
    Utility Methods
    ===============
 
    .. autosummary::
       :toctree: generated/
 
       download_all    -- Download all the dataset files to specified path.
       clear_cache     -- Clear cached dataset directory.
    ...
    ...
    ...

Testovací obrázek reprezentovaný ve stupních šedi

První datovou sadou, se kterou se v dnešním článku setkáme, je rastrový obrázek s rozlišením 512×512 pixelů, který obsahuje pouze intenzity pixelů. Jinými slovy je obrázek reprezentován ve stupních šedi. Tento obrázek získáme funkcí scipy.datasets.ascent(). Pokud je již obrázek uložen v cache, je vrácen ve formě n-dimenzionálního pole. V opačném případě je nejprve stažen z Internetu:

import numpy as np
import scipy.datasets as datasets
 
 načtení matice
ascent = datasets.ascent()
 
 zobrazení základních informací o matici
np.info(ascent)

Po spuštění tohoto příkladu by se měly vypsat všechny důležité informace o obrázku:

class:  ndarray
shape:  (512, 512)
strides:  (512, 1)
itemsize:  1
aligned:  True
contiguous:  True
fortran:  False
data pointer: 0x562eca70cf00
byteorder:  little
byteswap:  False
type: uint8

Z tohoto výpisu je patrné, že se jedná o dvourozměrné pole o rozměrech 512×512 prvků typu uint8, tj. každý pixel je reprezentován jedním bajtem (bez znaménka).

Vizualizace testovacího obrázku v nepravých barvách i ve stupních šedi

S rastrovými obrázky popř. s dalšími daty, která jsou reprezentována formou dvourozměrné matice, se v praxi můžeme setkat velmi často. Pro vizualizaci takto strukturovaných dat je možné využít knihovnu Matplotlib a její funkci imshow. Jedná se o přímočarou operaci, jejíž výsledek však nemusí být uspokojující (to ovšem do značné míry záleží na povaze dat):

import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
 načtení matice
ascent = datasets.ascent()
 
 zobrazení matice
plt.imshow(ascent)
 
 uložení matice do rastrového obrázku
plt.savefig("datasets_2.png")
 
 zobrazení grafu
plt.show()

Výsledkem bude obrázek v nepravých barvách:

V případě rastrového obrázku, který obsahuje pixely ve stupních šedi, je nutné při vykreslení použít nepovinný parametr cmap (color map) nastavený na hodnotu gray nebo binary nebo gray_r (což je, poněkud nepřesně řečeno, inverzní barvová paleta ke gray):

import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
 načtení matice
ascent = datasets.ascent()
 
 zobrazení matice s volbou barvové palety
plt.imshow(ascent, cmap="gray")
 
 uložení matice do rastrového obrázku
plt.savefig("datasets_3.png")
 
 zobrazení grafu
plt.show()

Nyní bude výsledek vypadat takto:

Můžete ovšem použít libovolnou další barvovou paletu z následujícího seznamu:

Testovací obrázek reprezentovaný v pravých barvách

Další datovou sadou, kterou můžeme přímo použít, je testovací obrázek s rozlišením 1024×768 pixelů, který je ovšem tentokrát reprezentovaný v pravých barvách (true color). Jedná se o obrázek mývala, jehož originál je dostupný na adrese https://pixnio.com/fauna-animals/raccoons/raccoon-procyon-lotor:

import numpy as np
import scipy.datasets as datasets
 
 načtení matice
raccoon = datasets.face()
 
 zobrazení základních informací o matici
np.info(raccoon)

Po spuštění tohoto velmi jednoduchého skriptu se zobrazí základní informace o datové sadě. Vidíme, že se jedná o trojrozměrné pole s 768×1024×3 prvky. Nejdříve je uveden počet řádků, poté počet pixelů na řádku a nakonec počet barvových složek. Typ pixelů je opět uint8, tedy bajt(y) bez znaménka:

class:  ndarray
shape:  (768, 1024, 3)
strides:  (3072, 3, 1)
itemsize:  1
aligned:  True
contiguous:  True
fortran:  False
data pointer: 0x7f44729fe030
byteorder:  little
byteswap:  False
type: uint8

Vizualizace obrázku s převodem na stupně šedi

Rastrový obrázek reprezentovaný v pravých barvách je možné zobrazit buď opět v pravých barvách nebo ve stupních šedi. Podívejme se nejprve na první případ, který lze řešit následovně:

import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
 načtení matice
raccoon = datasets.face()
 
 zobrazení matice s volbou barvové palety
plt.imshow(raccoon, cmap="gray")
 
 uložení matice do rastrového obrázku
plt.savefig("datasets_5.png")
 
 zobrazení grafu
plt.show()

Výsledek by měl vypadat takto:

Můžeme si ovšem vyžádat původní obrázek přímo ve stupních šedi, tj. reprezentovaný dvourozměrnou maticí:

import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
 načtení matice
raccoon = datasets.face(gray=True)
 
 zobrazení matice s volbou barvové palety
plt.imshow(raccoon, cmap="gray")
 
 uložení matice do rastrového obrázku
plt.savefig("datasets_5_gray.png")
 
 zobrazení grafu
plt.show()

Výsledek by měl nyní vypadat následovně:

Datová sada s jednorozměrným signálem

Poslední datovou sadou nabízenou knihovnou SciPy, se kterou se dnes seznámíme, je jednorozměrné pole (tedy jinými slovy vektor) s EKG (resp. jediným signálem EKG). Toto pole získáme snadno:

import numpy as np
import scipy.datasets as datasets
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 zobrazení základních informací o signálu
np.info(ekg)

Z výpisu základních informací o poli je patrné, že obsahuje 108000 prvků typu float64 (což například v jazyce C odpovídá typu double). Při měření EKG byla použita vzorkovací frekvence 360 Hz a záznam je pětiminutový. To znamená, že počet vzorků je roven 360×5×60=108000, což přesně odpovídá počtu prvků vektoru:

class:  ndarray
shape:  (108000,)
strides:  (8,)
itemsize:  8
aligned:  True
contiguous:  True
fortran:  True
data pointer: 0x7f9d22556010
byteorder:  little
byteswap:  False
type: float64

Vizualizace načteného signálu

Zobrazení jednorozměrného signálu je většinou snadné, protože můžeme použít funkci nazvanou plot naimportovanou přímo z balíčku Matplotlib. V případě, že budeme chtít zobrazit všech 108000 vzorků (což ovšem nemusí být nejlepší nápad), lze postupovat přímočaře:

import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 zobrazení signálu
plt.plot(ekg)
 
 uložení grafu s průběhem signálu
plt.savefig("ekg_1.png")
 
 zobrazení grafu
plt.show()

Výsledný průběh by měl vypadat následovně:

Základní úpravy zobrazení: část signálu popř. každý n-tý vzorek

Konkrétně v případě ukázkového EKG může být výhodné si zobrazit pouze část signálu. K tomuto účelu se dá využít mnoha mechanismů. Pravděpodobně nejjednodušší je využití přetíženého operátoru indexování [] (provedení řezu – slice) tak, jak je to ukázáno v dalším příkladu, který zobrazí pouze první část signálu:

import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 zobrazení části signálu
plt.plot(ekg[0:1000])
 
 uložení grafu s průběhem signálu
plt.savefig("ekg_2.png")
 
 zobrazení grafu
plt.show()

Graf s průběhem EKG bude nyní přehlednější, i když stále nebudeme vidět každý vzorek (muselo by se změnit rozlišení výsledného obrázku):

Alternativně si můžeme, například pro účely náhledu, nechat zobrazit každý n-tý vzorek. K tomuto účelu slouží opět operace řezu (slice), tentokrát doplněná o krok (step). Počáteční a koncový index se nemusí zadávat; pokud není zadán, doplní se za první index nula a signál se projde až do svého konce:

import scipy.datasets as datasets
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 zobrazení části signálu
plt.plot(ekg[::200])
 
 uložení grafu s průběhem signálu
plt.savefig("ekg_3.png")
 
 zobrazení grafu
plt.show()

Výsledek:

Balíček scipy.signal: zpracování signálů knihovnou SciPy

Jedním z nejpoužívanějších balíčků, které nalezneme v knihovně SciPy, je balíček nazvaný přímočaře scipy.signal. V tomto balíčku nalezneme funkce určené pro zpracování jednorozměrných i vícerozměrných signálů:

AME
    scipy.signal
 
DESCRIPTION
    =======================================
    Signal processing (:mod:`scipy.signal`)
    =======================================
 
    .. toctree::
       :hidden:
 
       signal.windows
 
    Convolution
    ===========
 
    .. autosummary::
       :toctree: generated/
 
       convolve           -- N-D convolution.
       correlate          -- N-D correlation.
       fftconvolve        -- N-D convolution using the FFT.
       oaconvolve         -- N-D convolution using the overlap-add method.
       convolve2d         -- 2-D convolution (more options).
       correlate2d        -- 2-D correlation (more options).
       sepfir2d           -- Convolve with a 2-D separable FIR filter.
       choose_conv_method -- Chooses faster of FFT and direct convolution methods.
       correlation_lags   -- Determines lag indices for 1D cross-correlation.
       ...
       ...
       ...

V dnešním článku se zaměříme na jednorozměrné signály; příště si ovšem ukážeme i práci se signály vícerozměrnými.

Funkce s realizací oken (window)

Při zpracování signálů se často pracuje s takzvanými okny. Jedná se o funkce, které mají mimo zadaný rozsah nulové hodnoty. Samotná funkce ale může být libovolná, dokonce ani nemusí být hladká. Nejjednodušším oknem je obdélníkové okno popř. okno trojúhelníkové, ovšem existuje celá řada dalších oken s více či méně spojitým průběhem. Jejich realizace (formou funkcí) je dostupná v podbalíčku nazvaném scipy.signal.windows:

NAME
    scipy.signal.windows
 
DESCRIPTION
    Window functions (:mod:`scipy.signal.windows`)
    ==============================================
 
    The suite of window functions for filtering and spectral estimation.
 
    .. currentmodule:: scipy.signal.windows
 
    .. autosummary::
       :toctree: generated/
 
       get_window              -- Convenience function for creating various windows.
 
       barthann                -- Bartlett-Hann window
       bartlett                -- Bartlett window
       blackman                -- Blackman window
       blackmanharris          -- Minimum 4-term Blackman-Harris window
       bohman                  -- Bohman window
       boxcar                  -- Boxcar window
       chebwin                 -- Dolph-Chebyshev window
       cosine                  -- Cosine window
       dpss                    -- Discrete prolate spheroidal sequences
       exponential             -- Exponential window
       flattop                 -- Flat top window
       gaussian                -- Gaussian window
       general_cosine          -- Generalized Cosine window
       general_gaussian        -- Generalized Gaussian window
       general_hamming         -- Generalized Hamming window
       hamming                 -- Hamming window
       hann                    -- Hann window
       kaiser                  -- Kaiser window
       kaiser_bessel_derived   -- Kaiser-Bessel derived window
       lanczos                 -- Lanczos window also known as a sinc window
       nuttall                 -- Nuttall's minimum 4-term Blackman-Harris window
       parzen                  -- Parzen window
       taylor                  -- Taylor window
       triang                  -- Triangular window
       tukey                   -- Tukey window

Průběh funkce, kterou okno reprezentuje, si pochopitelně můžeme zobrazit. Příkladem je Hannovo okno:

import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 filtr
win = signal.windows.hann(50)
 
 vykreslení výsledného signálu
plt.plot(win)
 
 uložení grafu s průběhem signálu
plt.savefig("window_1.png")
 
 zobrazení grafu
plt.show()

Takto by měl vypadat výsledný graf okenní funkce:

Podobně si můžeme nechat vypočítat a obrazit trojúhelníkové okno:

import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 filtr
win = signal.windows.triang(50)
 
 vykreslení výsledného signálu
plt.plot(win)
 
 uložení grafu s průběhem signálu
plt.savefig("window_2.png")
 
 zobrazení grafu
plt.show()

Výsledek:

A do třetice si ukažme výpočet a zobrazení Lanczosova okna:

import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 filtr
win = signal.windows.lanczos(50)
 
 vykreslení výsledného signálu
plt.plot(win)
 
 uložení grafu s průběhem signálu
plt.savefig("window_3.png")
 
 zobrazení grafu
plt.show()

Tentokrát by měl výsledek vypadat takto:

Konvoluce jednorozměrných signálů

Velmi často volanou operací, která se při zpracování signálů používá, je konvoluce (convolution). Ta může být provedena jak pro jednorozměrné signály (zvukové filtry atd.), tak i pro vícerozměrné signály (známé jsou konvoluční filtry v oblasti zpracování obrazu). Konvoluci dvou signálů, typicky zpracovávaného signálu a konvolučního jádra, lze vypočítat s využitím funkce nazvané přímočaře convolve:

Help on function convolve in module scipy.signal._signaltools:
 
convolve(in1, in2, mode='full', method='auto')
    Convolve two N-dimensional arrays.
 
    Convolve `in1` and `in2`, with the output size determined by the
    `mode` argument.
 
    Parameters
    ----------
    in1 : array_like
        First input.
    in2 : array_like
        Second input. Should have the same number of dimensions as `in1`.
    mode : str {'full', 'valid', 'same'}, optional
        A string indicating the size of the output:
 
        ``full``
           The output is the full discrete linear convolution
           of the inputs. (Default)
        ``valid``
           The output consists only of those elements that do not
           rely on the zero-padding. In 'valid' mode, either `in1` or `in2`
           must be at least as large as the other in every dimension.
        ``same``
           The output is the same size as `in1`, centered
           with respect to the 'full' output.
    method : str {'auto', 'direct', 'fft'}, optional
        A string indicating which method to use to calculate the convolution.
 
        ``direct``
           The convolution is determined directly from sums, the definition of
           convolution.
        ``fft``
        ...
        ...
        ...

Příklad konvoluce signálu s okny

V dnešním článku používáme v demonstračních příkladech především jednorozměrný signál se záznamem EKG. Tento signál využijeme i při výpočtu konvoluce, konkrétně konvoluce tohoto signálu s okenními funkcemi (resp. se signály, které jsou pomocí těchto funkcí vypočteny). První příklad vypočte konvoluci EKG s Hannovým oknem:

import scipy.datasets as datasets
import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 jen prvních 1000 vzorků
ekg = ekg[0:1000]
 
 filtr
win = signal.windows.hann(50)
 
 výpočet konvoluce
filtered = signal.convolve(ekg, win, mode="same") / sum(win)
 
 vykreslení výsledného signálu
plt.plot(filtered, "r")
 
 uložení grafu s průběhem signálu
plt.savefig("convolve_1d_1.png")
 
 zobrazení grafu
plt.show()

Výsledek:

Zobrazit si můžeme jak průběh původního signálu, tak i výsledku konvoluce:

import scipy.datasets as datasets
import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 jen prvních 1000 vzorků
ekg = ekg[0:1000]
 
 filtr
win = signal.windows.hann(50)
 
 výpočet konvoluce
filtered = signal.convolve(ekg, win, mode="same") / sum(win)
 
 vykreslení původního i výsledného signálu
plt.plot(ekg)
plt.plot(filtered, "r")
 
 uložení grafu s průběhem signálu
plt.savefig("convolve_1d_2.png")
 
 zobrazení grafu
plt.show()

Nyní získáme graf se dvěma průběhy:

Konvoluce EKG s trojúhelníkovým oknem:

import scipy.datasets as datasets
import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 jen prvních 1000 vzorků
ekg = ekg[0:1000]
 
 filtr
win = signal.windows.triang(50)
 
 výpočet konvoluce
filtered = signal.convolve(ekg, win, mode="same") / sum(win)
 
plt.plot(filtered, "r")
 
 uložení grafu s průběhem signálu
plt.savefig("convolve_1d_3.png")
 
 zobrazení grafu
plt.show()

Výsledek bude v tomto případě vypadat následovně:

Porovnání konvolucí s různými okenními funkcemi

Pochopitelně můžeme porovnat i výsledek konvolucí EKG s různými okenními funkcemi, a to jak numericky, tak i vizuálně:

import scipy.datasets as datasets
import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 jen prvních 1000 vzorků
ekg = ekg[0:1000]
 
 filtry
win1 = signal.windows.triang(50)
win2 = signal.windows.hann(50)
 
 výpočet konvoluce
filtered1 = signal.convolve(ekg, win1, mode="same") / sum(win1)
filtered2 = signal.convolve(ekg, win2, mode="same") / sum(win2)
 
plt.plot(filtered1, "b")
plt.plot(filtered2, "r")
 
 uložení grafu s průběhem signálu
plt.savefig("convolve_1d_4.png")
 
 zobrazení grafu
plt.show()

Výsledek porovnání v numerické podobě (pouze několik prvních vzorků):

-0.095 -0.094
-0.103 -0.102
-0.110 -0.110
-0.117 -0.117
-0.123 -0.125
-0.130 -0.132
-0.136 -0.139
-0.142 -0.146
-0.147 -0.153
-0.153 -0.159
-0.158 -0.164
-0.163 -0.170
-0.168 -0.175
-0.173 -0.179
-0.177 -0.183
-0.180 -0.187
-0.184 -0.190
-0.187 -0.192
-0.190 -0.195
-0.193 -0.197
-0.195 -0.198
-0.197 -0.199
-0.198 -0.201
-0.200 -0.201
-0.201 -0.202
-0.201 -0.203
-0.201 -0.203
-0.202 -0.203
-0.202 -0.204
-0.202 -0.204
-0.201 -0.203
-0.201 -0.203
-0.201 -0.203
-0.200 -0.202
-0.200 -0.202
-0.199 -0.201
-0.198 -0.200
-0.197 -0.199
-0.196 -0.198
-0.194 -0.196
-0.193 -0.195
-0.191 -0.193
-0.188 -0.191
-0.186 -0.189
-0.183 -0.186
-0.179 -0.184
-0.176 -0.180
-0.172 -0.177
-0.168 -0.173
-0.164 -0.169
-0.159 -0.165
-0.155 -0.160
-0.150 -0.155
-0.144 -0.149
-0.139 -0.143
-0.133 -0.137
-0.127 -0.130
-0.120 -0.123
-0.114 -0.116
-0.107 -0.109
-0.101 -0.102
-0.094 -0.094
-0.087 -0.087
-0.081 -0.079
-0.074 -0.072
-0.068 -0.065
-0.062 -0.058
-0.057 -0.052
-0.052 -0.046
-0.047 -0.040
-0.043 -0.035
-0.039 -0.030
-0.035 -0.026
-0.033 -0.023
-0.030 -0.020
-0.029 -0.019
-0.027 -0.018
-0.027 -0.017
-0.026 -0.018
-0.026 -0.018
-0.026 -0.020
-0.027 -0.022
-0.028 -0.024
-0.029 -0.026
-0.030 -0.029
-0.032 -0.032
-0.033 -0.035
-0.034 -0.037
-0.036 -0.040
-0.037 -0.042
-0.039 -0.044
-0.040 -0.046
-0.041 -0.047
-0.042 -0.048
-0.043 -0.048
-0.043 -0.048
-0.043 -0.048
-0.043 -0.047
-0.042 -0.046
-0.041 -0.045

Výsledek porovnání ve vizuální podobě:

Převzorkování (resampling)

Další poměrně často používanou operací je převzorkování, který je v knihovně SciPy realizován funkcí nazvanou scipy.signal.resample. Tuto operaci lze provádět několika způsoby (v časové či frekvenční oblasti) a může se současně aplikovat i vhodné okno:

resample(x, num, t=None, axis=0, window=None, domain='time')
    Resample `x` to `num` samples using the Fourier method along the given `axis`.
 
    The resampling is performed by shortening or zero-padding the FFT of `x`. This has
    the advantages of providing an ideal antialiasing filter and allowing arbitrary
    up- or down-sampling ratios. The main drawback is the requirement of assuming `x`
    to be a periodic signal.
 
    Parameters
    ----------
    x : array_like
        The input signal made up of equidistant samples. If `x` is a multidimensional
        array, the parameter `axis` specifies the time/frequency axis. It is assumed
        here that ``n_x = x.shape[axis]`` specifies the number of samples and ``T`` the
        sampling interval.
    num : int
        The number of samples of the resampled output signal. It may be larger or
        smaller than ``n_x``.
    t : array_like, optional
        If `t` is not ``None``, then the timestamps of the resampled signal are also
        returned. `t` must contain at least the first two timestamps of the input
        signal `x` (all others are ignored). The timestamps of the output signal are
        determined by ``t[0] + T * n_x / num * np.arange(num)`` with
        ``T = t[1] - t[0]``. Default is ``None``.
    axis : int, optional
        The time/frequency axis of `x` along which the resampling take place.
        The Default is 0.
    window : array_like, callable, string, float, or tuple, optional
        If not ``None``, it specifies a filter in the Fourier domain, which is applied
        before resampling. I.e., the FFT ``X`` of `x` is calculated by
        ``X = W * fft(x, axis=axis)``. ``W`` may be interpreted as a spectral windowing
        function ``W(f_X)`` which consumes the frequencies ``f_X = fftfreq(n_x, T)``.
    ...
    ...
    ...

První (naivní) ukázka resamplingu

Pokusme se nyní získat prvních 250 vzorků EKG a převzorkovat tuto část signálu do podoby pouze s padesáti vzorky. Vlastně tak snižujeme vzorkovací frekvenci, ovšem na již existujícím signálu. Výsledné signály (původní i převzorkovaný) následně vykreslíme:

import numpy as np
import scipy.datasets as datasets
import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 jen prvních 250 vzorků
ekg = ekg[0:250]
 
 resampling
resampled = signal.resample(ekg, 50)
 
 vykreslení výsledného signálu
plt.plot(ekg, "b")
plt.plot(resampled, "r")
 
 uložení grafu s průběhem signálu
plt.savefig("resample_1.png")
 
 zobrazení grafu
plt.show()

Výsledek ovšem nebude dokonalý, protože na x-ové ose nepoužíváme stejné jednotky (nový signál má mnohem menší množství vzorků, než signál původní):

Oprava jednotek na x-ové ose

Našim cílem bude zobrazit jak původní úsek EKG, tak i jeho převzorkovanou variantu se stejnými časovými jednotkami na x-ové ose. To se provede relativně snadno – konkrétně předáním vektoru s x-ovými souřadnicemi do funkce matplotlib.pyplot.plot. Vektor musí pro oba signály začínat a končit stejnými hodnotami, ovšem počet hodnot (a krok mezi nimi) bude pochopitelně odlišný:

import numpy as np
import scipy.datasets as datasets
import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 jen prvních 250 vzorků
ekg = ekg[0:250]
 
 resampling
resampled = signal.resample(ekg, 50)
 
 vykreslení výsledného signálu
plt.plot(np.linspace(0, 1000, len(ekg)), ekg, "b")
plt.plot(np.linspace(0, 1000, len(resampled)), resampled, "r")
 
 uložení grafu s průběhem signálu
plt.savefig("resample_2.png")
 
 zobrazení grafu
plt.show()

Nyní bude výsledek téměr korektní, ovšem jen téměř (liší se koncové body a proto jsou signály v částečně odlišném měřítku):

Korektní varianta ve vektorech s x-ovými souřadnicemi vynechává poslední hodnotu, takže signály nebudou končit ve stejném čase (to nám však nevadí, spíš naopak – nedojde k protažení jednoho ze signálů):

import numpy as np
import scipy.datasets as datasets
import scipy.signal as signal
 
import matplotlib.pyplot as plt
 
 načtení signálu
ekg = datasets.electrocardiogram()
 
 jen prvních 250 vzorků
ekg = ekg[0:250]
 
 resampling
resampled = signal.resample(ekg, 50)
 
 vykreslení výsledného signálu
plt.plot(np.linspace(0, 1000, len(ekg), endpoint=False), ekg, "b")
plt.plot(np.linspace(0, 1000, len(resampled), endpoint=False), resampled, "r")
 
 uložení grafu s průběhem signálu
plt.savefig("resample_3.png")
 
 zobrazení grafu
plt.show()

Grafický výsledek:

Můžeme provést i numerické ověření:

import numpy as np
 
x1 = np.linspace(0, 1000, 10)
x2 = np.linspace(0, 1000, 20)
 
print(x1)
print(x2)
 
x1 = np.linspace(0, 1000, 10, endpoint=False)
x2 = np.linspace(0, 1000, 20, endpoint=False)
 
print(x1)
print(x2)

Obě dvojice vektorů se budou lišit. Povšimněte si sice stejné délky vektorů, ovšem odlišného kroku i jiných koncových hodnot:

[   0.          111.11111111  222.22222222  333.33333333  444.44444444
  555.55555556  666.66666667  777.77777778  888.88888889 1000.        ]
[   0.           52.63157895  105.26315789  157.89473684  210.52631579
  263.15789474  315.78947368  368.42105263  421.05263158  473.68421053
  526.31578947  578.94736842  631.57894737  684.21052632  736.84210526
  789.47368421  842.10526316  894.73684211  947.36842105 1000.        ]

vs.:

[  0. 100. 200. 300. 400. 500. 600. 700. 800. 900.]
[  0.  50. 100. 150. 200. 250. 300. 350. 400. 450. 500. 550. 600. 650.
 700. 750. 800. 850. 900. 950.]

Manipulace se signály ve frekvenční oblasti

Funkce, které jsme si až doposud ukázali, jsou určeny pro zpracování signálů v časové oblasti. Ovšem knihovna SciPy (pochopitelně) obsahuje i funkce, které dokážou zpracovat signály v oblasti frekvenční. Do této kategorie patří především funkce provádějící rychlou Fourierovu transformaci (FFT – Fast Fourier Transform), diskrétní sinovou transformací (tu jsem popravdě nikdy neviděl použitou v praxi) a taktéž diskrétní kosinovou transformaci (DCT – Discrete Cosine Transform), s níž jsme se kdysi seznámili v souvislosti s grafickým formátem JPEG (viz též odkazy uvedené na konci článku). S těmito transformacemi, které je možné provádět jak s jednorozměrnými signály, tak i se signály vícerozměrnými, se podrobněji seznámíme v navazujícím článku.

Repositář s demonstračními příklady

Odkazy na Internetu

1. SciPy homepage
   https://scipy.org/
2. SciPy (Wikipedia)
   https://en.wikipedia.org/wiki/SciPy
3. The Hertzsprung–Russell diagram
   https://scipython.com/book2/chapter-9-data-analysis-with-pandas/problems/p92/the-hertzsprung-russell-diagram/
4. Linear algebra (scipy.linalg)
   https://docs.scipy.org/doc/sci­py/reference/linalg.html
5. Frequently Asked Questions – SciPy
   https://scipy.org/faq/
6. SciPy – Introduction
   https://www.tutorialspoin­t.com/scipy/scipy_introduc­tion.htm
7. LAPACK — Linear Algebra PACKage
   https://www.netlib.org/lapack/
8. LAPACK (Wikipedia)
   https://en.wikipedia.org/wiki/LAPACK
9. LAPACK na GitHubu
   https://github.com/Reference-LAPACK/lapack
10. SciPy in Python
    https://pythonguides.com/scipy/
11. scipy.linalg.det
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.linalg.det.html#scipy.li­nalg.det
12. scipy.linalg.inv
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.linalg.inv.html#scipy.li­nalg.inv
13. scipy.linalg.solve
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.linalg.inv.html#scipy.li­nalg.inv
14. Algebra
    https://cs.wikipedia.org/wiki/Algebra
15. Lineární algebra
    https://cs.wikipedia.org/wi­ki/Line%C3%A1rn%C3%AD_alge­bra
16. Lineární rovnice
    https://cs.wikipedia.org/wi­ki/Line%C3%A1rn%C3%AD_rov­nice
17. Soustava lineárních rovnic
    https://cs.wikipedia.org/wi­ki/Soustava_line%C3%A1rn%C3%AD­ch_rovnic
18. Norma matice
    https://cs.wikipedia.org/wi­ki/Norma_matice
19. Matrix norm
    https://en.wikipedia.org/wi­ki/Matrix_norm
20. Norma (vektoru)
    https://cs.wikipedia.org/wi­ki/Norma_(matematika)
21. Frobeniův skalární součin
    https://cs.wikipedia.org/wi­ki/Frobeni%C5%AFv_skal%C3%A1rn%C3%AD_sou­%C4%8Din
22. BLAS (Basic Linear Algebra Subprograms)
    https://www.netlib.org/blas/
23. Basic Linear Algebra Subprograms (Wikipedia)
    https://en.wikipedia.org/wi­ki/Basic_Linear_Algebra_Sub­programs
24. Operace s daty uloženými v binárních souborech v knihovnách NumPy a Pandas
    https://www.root.cz/clanky/operace-s-daty-ulozenymi-v-binarnich-souborech-v-knihovnach-numpy-a-pandas/
25. Operace s daty uloženými v binárních souborech v knihovnách NumPy a Pandas (dokončení)
    https://www.root.cz/clanky/operace-s-daty-ulozenymi-v-binarnich-souborech-v-knihovnach-numpy-a-pandas-dokonceni/
26. NumPy Home Page
    http://www.numpy.org/
27. NumPy v1.10 Manual
    http://docs.scipy.org/doc/num­py/index.html
28. NumPy (Wikipedia)
    https://en.wikipedia.org/wiki/NumPy
29. OpenBLAS: An optimized BLAS library
    https://www.openblas.net/
30. Integrovaná vývojová prostředí ve Fedoře: IPython a IPython Notebook
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-ipython-a-ipython-notebook/
31. Integrovaná vývojová prostředí ve Fedoře: praktické použití IPython Notebooku a knihovny Numpy
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-prakticke-pouziti-ipython-notebooku-a-knihovny-numpy/
32. Integrovaná vývojová prostředí ve Fedoře: praktické použití IPython Notebooku a knihovny Numpy (2.část)
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-prakticke-pouziti-ipython-notebooku-a-knihovny-numpy-2-cast/
33. Integrovaná vývojová prostředí ve Fedoře: vykreslování grafů s využitím knihoven Numpy a matplotlib
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-vykreslovani-grafu-s-vyuzitim-knihoven-numpy-a-matplotlib/
34. Integrovaná vývojová prostředí ve Fedoře: vykreslování grafů s využitím knihoven Numpy a matplotlib (2.část)
    http://mojefedora.cz/integrovana-vyvojova-prostredi-ve-fedore-vykreslovani-grafu-s-vyuzitim-knihoven-numpy-a-matplotlib-2-cast/
35. Piecewise linear interpolation
    https://docs.scipy.org/doc//sci­py/tutorial/interpolate/1D­.html
36. Statistical functions (scipy.stats)
    https://docs.scipy.org/doc/sci­py/reference/stats.html
37. scipy.stats.linregress
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.stats.linregress.html
38. numpy.poly1d
    https://numpy.org/doc/sta­ble/reference/generated/num­py.poly1d.html
39. scipy.optimize.curve_fit
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.optimize.curve_fit.html#sci­py.optimize.curve_fit
40. scipy.interpolate.make_smoothing_spline
    https://docs.scipy.org/doc/sci­py/reference/generated/sci­py.interpolate.make_smoot­hing_spline.html#scipy.in­terpolate.make_smoothing_spli­ne
41. Famous Curves Index
    https://mathshistory.st-andrews.ac.uk/Curves/
42. Curve (Wikipedia)
    https://en.wikipedia.org/wiki/Curve
43. Mathematical curves
    https://www.2dcurves.com/index.html
44. Curves (Wolfram MathWorld)
    https://mathworld.wolfram­.com/topics/Curves.html
45. Smooth Curve (Wolfram MathWorld)
    https://mathworld.wolfram­.com/SmoothCurve.html
46. Spirals (Wolfram MathWorld)
    https://mathworld.wolfram­.com/topics/Spirals.html
47. An Interactive Introduction to Splines
    https://ibiblio.org/e-notes/Splines/Intro.htm
48. Algebraic curve
    https://en.wikipedia.org/wi­ki/Algebraic_curve
49. Transcendental curve
    https://en.wikipedia.org/wi­ki/Transcendental_curve
50. Algebraic Curves
    https://mathworld.wolfram­.com/topics/AlgebraicCurves­.html
51. Elliptic Curves
    https://mathworld.wolfram­.com/topics/EllipticCurves­.html
52. Curves we (mostly) don't learn in high school (and applications)
    https://www.youtube.com/wat­ch?v=3izFMB91K_Q
53. Catenary arch
    https://en.wikipedia.org/wi­ki/Catenary_arch
54. Parabolic arch
    https://en.wikipedia.org/wi­ki/Parabolic_arch
55. Wattova křivka
    https://www.geogebra.org/m/gNh4bW9r
56. Fifty Famous Curves, Lots of Calculus Questions, And a Few Answers
    https://elepa.files.wordpres­s.com/2013/11/fifty-famous-curves.pdf
57. Faux, I.D. a Pratt, M.J.: Computational Geometry for Design and Manufacture,
    Ellis Horwood Ltd., Wiley & Sons, 1979
58. Wallace A.: Differential Topology,
    Benjamin/Cummings Co., Reading, Massachussetts, USA, 1968
59. Glossary of Bridge Terminology
    http://sdrc.lib.uiowa.edu/en­g/bridges/WaddellGlossary/Glos­sC.htm
60. Brachistochrona
    https://cs.wikipedia.org/wi­ki/Brachistochrona
61. Missions: Cassini
    https://solarsystem.nasa.gov/mis­sions/cassini/overview/
62. Giovanni Domenico Cassini
    https://en.wikipedia.org/wi­ki/Giovanni_Domenico_Cassi­ni
63. Cassini Ovals
    https://mathworld.wolfram­.com/CassiniOvals.html
64. Geocentrismus
    https://cs.wikipedia.org/wi­ki/Geocentrismus
65. Who was Giovanni Cassini?
    https://www.universetoday­.com/130823/who-was-giovanni-cassini/
66. Special plane curves
    http://xahlee.info/Special­PlaneCurves_dir/ConicSecti­ons_dir/conicSections.html
67. Interpolace
    https://mathonline.fme.vut­br.cz/pg/Algoritmy/05_APROX_KRIV­KY.htm
68. Lagrange Polynomial Interpolation
    https://pythonnumericalmet­hods.berkeley.edu/notebook­s/chapter17.04-Lagrange-Polynomial-Interpolation.html
69. Python Program for Lagrange Interpolation Method (with Output)
    https://www.codesansar.com/numerical-methods/python-program-lagrange-interpolation-method.htm
70. Smooth Paths Using Catmull-Rom Splines
    https://qroph.github.io/2018/07/30/smo­oth-paths-using-catmull-rom-splines.html
71. Lecture 11: Linear Interpolation Again – Bézier Curves
    http://www.math.kent.edu/~re­ichel/courses/intr.num.com­p.1/fall09/lecture12/bez.pdf
72. Geometrie/Úvod do křivek
    https://cs.wikibooks.org/wi­ki/Geometrie/%C3%9Avod_do_k%C5%99i­vek
73. B-Spline Curves and Surfaces (1)
    http://www.cad.zju.edu.cn/ho­me/zhx/GM/006/00-bscs1.pdf
74. Praktické ukázky možností aplikace Mandelbulber při tvorbě animací
    https://www.root.cz/clanky/prakticke-ukazky-moznosti-aplikace-mandelbulber-pri-tvorbe-animaci/
75. Kochanek–Bartels spline
    https://en.wikipedia.org/wi­ki/Kochanek%E2%80%93Bartel­s_spline
76. class KochanekBartels
    https://splines.readthedoc­s.io/en/latest/_modules/spli­nes.html#KochanekBartels
77. Konvoluce
    https://cs.wikipedia.org/wi­ki/Konvoluce
78. Korelace
    https://cs.wikipedia.org/wi­ki/Korelace
79. Fourierova transformace
    https://cs.wikipedia.org/wi­ki/Fourierova_transformace
80. Lanczos resampling
    https://en.wikipedia.org/wi­ki/Lanczos_resampling
81. Obrázek se schodištěm
    https://pixnio.com/people/accent-to-the-top
82. Obrázek s mývalem
    https://pixnio.com/fauna-animals/raccoons/raccoon-procyon-lotor
83. Sample-rate conversion
    https://en.wikipedia.org/wiki/Sample-rate_conversion
84. Seriál Grafické formáty
    https://www.root.cz/serialy/graficke-formaty/
85. JPEG – král rastrových grafických formátů?
    https://www.root.cz/clanky/jpeg-kral-rastrovych-grafickych-formatu/
86. Ztrátová komprese obrazových dat pomocí JPEG
    https://www.root.cz/clanky/ztratova-komprese-obrazovych-dat-pomoci-jpeg/
87. Programujeme JPEG: transformace a podvzorkování barev
    https://www.root.cz/clanky/pro­gramujeme-jpeg-transformace-a-podvzorkovani-barev/
88. Programujeme JPEG: diskrétní kosinová transformace (DCT)
    https://www.root.cz/clanky/pro­gramujeme-jpeg-diskretni-kosinova-transformace-dct/
89. Programujeme JPEG: Kvantizace DCT koeficientů
    https://www.root.cz/clanky/pro­gramujeme-jpeg-kvantizace-dct-koeficientu/
90. Programujeme JPEG: Huffmanovo kódování kvantovaných DCT složek
    https://www.root.cz/clanky/pro­gramujeme-jpeg-huffmanovo-kodovani-kvantovanych-dct-slozek/
91. Programujeme JPEG: Interní struktura souborů typu JFIF/JPEG
    https://www.root.cz/clanky/pro­gramujeme-jpeg-interni-struktura-souboru-typu-jfifjpeg/
92. Programujeme JPEG: Načtení informací ze souborů typu JFIF/JPEG
    https://www.root.cz/clanky/pro­gramujeme-jpeg-nacteni-informaci-ze-souboru-typu-jfifjpeg/
93. Programujeme JPEG: Progresivní JPEG a informace EXIF
    https://www.root.cz/clanky/pro­gramujeme-jpeg-progresivni-jpeg-a-informace-exif/