Základní funkce pro zpracování signálů v knihovně SciPy: FFT a její varianty

Ve čtvrtém článku o knihovně SciPy se budeme věnovat dalším funkcím určeným pro zpracování signálů. Minule jsme si ukázali funkce, které 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í (resp. signály transformovat do této oblasti).

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).

Výpočet rychlé Fourierovy transformace

Rychlá Fourierova transformace (FFT) je založena na následujícím výpočtu, který transformuje signál x[n] z časové oblasti na signál X[k] v oblasti frekvenční:

X[k] = Σn=0N−1 x[n] · e−i·2π·k·n/N,  pro k = 0,...,N−1.

Výpočet e^{výraz s i} proběhne přesně podle známého Eulerova vzorce, který vlastně ve zkrácené podobě vyjadřuje vztah mezi trojicí Taylorových rozvojů (rozvoje exponenciální funkce o základu e a rozvoje sinu a kosinu):

eiθ = cos θ + i sin θ

Podívejme se nyní, jakým způsobem se rychlá Fourierova transformace vypočítá prakticky, a to konkrétně pro vstupní signál se čtyřmi vzorky, tedy pro N=4. Úprava původního vzorce pro N=4 je v tomto případě triviální:

X[k] = Σn=03 x[n] · e-i·2π·k·n/N, pro k = 0..3.

V dalším kroku předpokládejme, že vzorky vstupního signálu x budou mít hodnoty [1, 2, 3, 4]. Dosadíme tyto hodnoty do vzorce pro k = 0..3 a dostaneme následující rovnice:

X[0] = 1 + 2 + 3 + 4
X[1] = 1 + 2·e-iπ/2 + 3·e-iπ + 4·e-i3π/2
X[2] = 1 + 2·e-iπ + 3·e-i2π + 4·e-i3π
X[3] = 1 + 2·e-i3π/2 + 3·e-i3π + 4·e-i9π/2

Nyní musíme použít výše uvedený Eulerův vzorec a vypočítat jednotlivé členy, ve kterých se vyskytuje e^xi:

e-jπ/2 = cos(−π/2) + j·sin(−π/2) = 0 − j = −j
e-iπ = cos(−π) + i sin(−π) = cos π + i·(−sin π) = −1 + 0i = −1
e-i3π/2 = cos(−3π/2) + i sin(−3π/2) = cos(3π/2) − i sin(3π/2) = 0 − i(−1) = i
e-i2π = cos(−2π) + i sin(−2π) = cos(2π) − i sin(2π) = 1 + 0i = 1
e-i9π/2 = cos(−π/2) + i sin(−π/2) = 0 − i = −i
...
...
...

Tedy v našem konkrétním případě:

X[0] = 1 + 2 + 3 + 4 = 10
X[1] = 1 + 2·e-iπ/2 + 3·e-iπ + 4·e-i3π/2 = 1 + 2(−i) + 3(−1) + 4(i) = 1 − 2i − 3 + 4i = −2 + 2i
X[2] = 1 + 2·e-iπ + 3·e-i2π + 4·e-i3π = 1 + 2(−1) + 3(1) + 4(−1) = 1 − 2 + 3 − 4 = −2
X[3] = 1 + 2·e-i3π/2 + 3·e-i3π + 4·e-i9π/2 = 1 + 2(i) + 3(−1) + 4(−i) = 1 + 2i − 3 − 4i = −2 − 2i

Výsledkem tedy je čtveřice hodnot X[0] až X[3] vypsaná do následující tabulky:

Funkce fft

Pro výpočet rychlé Fourierovy transformace je v knihovně SciPy určena funkce nazvaná jednoduše fft, která je uložena v podbalíčku scipy.fft:

fft(x, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, *, plan=None)
    Compute the 1-D discrete Fourier Transform.
 
    This function computes the 1-D *n*-point discrete Fourier
    Transform (DFT) with the efficient Fast Fourier Transform (FFT)
    algorithm [1]_.
 
    Parameters
    ----------
    x : array_like
        Input array, can be complex.
    n : int, optional
        Length of the transformed axis of the output.
        If `n` is smaller than the length of the input, the input is cropped.
        If it is larger, the input is padded with zeros. If `n` is not given,
        the length of the input along the axis specified by `axis` is used.
    axis : int, optional
        Axis over which to compute the FFT. If not given, the last axis is
        used.
    norm : {"backward", "ortho", "forward"}, optional
        Normalization mode. Default is "backward", meaning no normalization on
        the forward transforms and scaling by ``1/n`` on the `ifft`.
        "forward" instead applies the ``1/n`` factor on the forward transform.
        For ``norm="ortho"``, both directions are scaled by ``1/sqrt(n)``.
 
        .. versionadded:: 1.6.0
           ``norm={"forward", "backward"}`` options were added
    overwrite_x : bool, optional
        If True, the contents of `x` can be destroyed; the default is False.
        See the notes below for more details.
    workers : int, optional
        Maximum number of workers to use for parallel computation. If negative,

Na základě toho, jaký konkrétní backend je při výpočtu zvolen, se výpočet provede buď na CPU nebo na GPU:

Ověření výpočtu zavoláním funkce scipy.fft.fft

Ve druhé kapitole uvedený výpočet, který byl prozatím proveden ručně, si pochopitelně můžeme ověřit na počítači, a to konkrétně právě s využitím knihovny Scipy. FFT vypočteme zavoláním funkce scipy.fft.fft, které předáme vektor obsahující stejné hodnoty, jaké jsme použili výše, tedy konkrétně [1, 2, 3, 4]. Hodnoty reprezentující signál lze předat v různém formátu, my využijeme n-rozměrné pole knihovny NumPy:

from scipy.fft import fft
 
import numpy as np
 
# signál v časové oblasti
t = np.array([1, 2, 3, 4])
 
# rychlá Fourierova transformace
f = fft(t)
 
for x in f:
    print(x)

Výsledky, které tento jednoduchý skript po svém spuštění vypíše, si pochopitelně můžeme porovnat s ručně vypočtenými hodnotami:

(10-0j)
(-2+2j)
(-2-0j)
(-2-2j)

To přesně odpovídá tabulce uvedené ve druhé kapitole.

Inverzní rychlá Fourierova transformace

Kromě transformace signálu z časové oblasti do oblasti frekvenční s využitím (přímé) rychlé Fourierovy transformace se v praxi setkáme i s opačným požadavkem, tedy konkrétně s požadavkem na transformaci z oblasti frekvenční zpět do oblasti časové. Je to vlastně logické, protože například potřebujeme provést jen několik manipulací ve frekvenční oblasti (filtry), ovšem jak vstupem do výpočtu, tak i jeho výsledkem, je diskrétní signál v oblasti časové (tedy jednotlivé vzorky – samply).

Zpětná rychlá Fourierova transformace se vypočte funkcí nazvanou ifft, která patří do stejného podbalíčku jako funkce fft, s níž jsme se již ve stručnosti seznámili:

ifft(x, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, *, plan=None)
    Compute the 1-D inverse discrete Fourier Transform.
 
    This function computes the inverse of the 1-D *n*-point
    discrete Fourier transform computed by `fft`.  In other words,
    ``ifft(fft(x)) == x`` to within numerical accuracy.
 
    The input should be ordered in the same way as is returned by `fft`,
    i.e.,
 
    * ``x[0]`` should contain the zero frequency term,
    * ``x[1:n//2]`` should contain the positive-frequency terms,
    * ``x[n//2 + 1:]`` should contain the negative-frequency terms, in
      increasing order starting from the most negative frequency.
 
    For an even number of input points, ``x[n//2]`` represents the sum of
    the values at the positive and negative Nyquist frequencies, as the two
    are aliased together. See `fft` for details.

Tuto funkci si můžeme poměrně snadno otestovat. Nejprve převedeme nám již známý signál [1, 2, 3, 4] do frekvenční oblasti a získáme čtveřici komplexních hodnot. A následně provedeme zpětnou transformaci, jejímž výsledkem by měl být původní signál:

from scipy.fft import fft, ifft
 
import numpy as np
 
t = np.array([1, 2, 3, 4])
 
# rychlá Fourierova transformace
f = fft(t)
 
for x in f:
    print(x)
 
print("-" * 10)
 
t2 = ifft(f)
 
for x in t2:
    print(x)
 

Tento skript po svém spuštění nejdříve vypíše výsledek dopředné rychlé Fourierovy transformace a následně i výsledek zpětné rychlé Fourierovy transformace:

(10-0j)
(-2+2j)
(-2-0j)
(-2-2j)
----------
(1+0j)
(2+0j)
(3-0j)
(4+0j)

Výsledek opravdu odpovídá vstupnímu signálu, protože sice dostaneme vektor komplexních hodnot, ovšem jejich imaginární složky jsou nulové.

Vizualizace FFT

Vzhledem k tomu, že výsledkem výpočtu rychlé Fourierovy transformace je u jednorozměrných signálů vektor s komplexními hodnotami, je možné si takový vektor zobrazit formou grafu resp. dvou grafů (průběh reálných složek a složek imaginárních, popř. absolutní hodnota a argument). Podívejme se ovšem, jaké hodnoty získáme například po výpočtu FFT kosinusu:

from scipy.fft import fft
 
import numpy as np
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100))
 
# rychlá Fourierova transformace
f = fft(t)
 
for x in f:
    print(x)

Teoreticky by měl být ve výsledném vektoru jeden nenulový prvek a další prvky nulové (nebo téměř nulové), ovšem ve skutečnosti získáme tento vektor:

(1.0000000000000004-0j)
(50.21819887226786+1.578170477977807j)
(-0.3415180830425068-0.021486496555723024j)
(-0.1271061001420848-0.012015063976087002j)
(-0.06739661481092062-0.008514172458434345j)
(-0.04185057087113764-0.006628479244688276j)
(-0.028480087577818183-0.005432867265545906j)
(-0.020583129025292392-0.004600874438040437j)
...
...
...
(-0.028480087577818183+0.005432867265545906j)
(-0.04185057087113764+0.006628479244688276j)
(-0.06739661481092062+0.008514172458434345j)
(-0.1271061001420848+0.012015063976087002j)
(-0.3415180830425068+0.021486496555723024j)
(50.21819887226786-1.578170477977807j)

Přepočet jednotek po FFT vypadá zhruba takto:

N = 100            ; počet samplů
t = 1s             ; doba samplování
R = 100/1 = 100 Hz ; samplovací frekvence
 
k = 0..99
f=R k/N            ; hodnoty (Hz), které by měly být na horizontální ose

Samozřejmě si můžeme průběh FFT kosinu vizualizovat:

from scipy.fft import fft
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100))
 
# rychlá Fourierova transformace
f = fft(t)
 
# vykreslení výsledného signálu
plt.plot(f)
 
# uložení grafu s průběhem signálu
plt.savefig("fft_4.png")

Výsledek bude vypadat následovně:

Zobrazení průběhu pouze s kladnými frekvencemi

Funkce fft ve skutečnosti vrací vektor, který v prvním prvku obsahuje stejnosměrnou složku, prvky A[1:n/2] hodnoty pro kladné frekvence a prvky A[n/2:] hodnoty pro frekvence záporné. Při vizualizaci je vhodné tyto záporné frekvence, tedy druhou polovinu vektoru, zahodit, k čemuž můžeme použít operátor řezu (slice):

from scipy.fft import fft
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100))
 
# rychlá Fourierova transformace
f = fft(t)[1:len(t)//2]
 
# vykreslení výsledného signálu
plt.plot(f)
 
# uložení grafu s průběhem signálu
plt.savefig("fft_5.png")

Výsledek bude lépe odpovídat realitě:

Funkce fftshift

V některých situacích nás ovšem mohou zajímat jak hodnoty odpovídající kladným frekvencím, tak i hodnoty odpovídající frekvencím záporným. Typický požadavek je následující: stejnosměrná složka má být zobrazena uprostřed, nalevo od ní záporné frekvence a napravo od ní naopak frekvence kladné. Samozřejmě by bylo možné ručně přesunout prvky výsledného vektoru, ovšem ve skutečnosti nám knihovna SciPy nabízí hotové řešení představované funkcí s poměrně výmluvným názvem fftshift:

fftshift(x, axes=None)
    Shift the zero-frequency component to the center of the spectrum.
 
    This function swaps half-spaces for all axes listed (defaults to all).
    Note that ``y[0]`` is the Nyquist component only if ``len(x)`` is even.
 
    Parameters
    ----------
    x : array_like
        Input array.
    axes : int or shape tuple, optional
        Axes over which to shift.  Default is None, which shifts all axes.
 
    Returns
    -------
    y : ndarray
        The shifted array.
 
    See Also
    --------
 
    :func:`ifftshift`
        The inverse of `fftshift`.
 
 
    Notes
    -----
 
    **Array API Standard Support**
 
    `fftshift` has experimental support for Python Array API Standard compatible
    backends in addition to NumPy. Please consider testing these features
    by setting an environment variable ``SCIPY_ARRAY_API=1`` and providing
    CuPy, PyTorch, JAX, or Dask arrays as array arguments. The following
    combinations of backend and device (or other capability) are supported.
 
    ====================  ====================  ====================
    Library               CPU                   GPU
    ====================  ====================  ====================
    NumPy                 ✅                     n/a
    CuPy                  n/a                   ✅
    PyTorch               ✅                     ✅
    JAX                   ✅                     ✅
    Dask                  ✅                     n/a
    ====================  ====================  ====================
 
        See :ref:`dev-arrayapi` for more information.
 
    Examples
    --------
    >>> import numpy as np
    >>> freqs = np.fft.fftfreq(10, 0.1)
    >>> freqs
    array([ 0.,  1.,  2., ..., -3., -2., -1.])
    >>> np.fft.fftshift(freqs)
    array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.])
 
    Shift the zero-frequency component only along the second axis:
 
    >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
    >>> freqs
    array([[ 0.,  1.,  2.],
           [ 3.,  4., -4.],
           [-3., -2., -1.]])
    >>> np.fft.fftshift(freqs, axes=(1,))
    array([[ 2.,  0.,  1.],
           [-4.,  3.,  4.],
           [-1., -3., -2.]])

Ukázka použití funkce fftshift

Funkce fftshift se aplikuje na výsledek získaný funkcí fft. V tom nejjednodušším případě může výpočet a vizualizace FFT vypadat takto:

from scipy.fft import fft, fftshift
 
import numpy as np
import matplotlib.pyplot as plt
 
# signál
t = np.cos(np.linspace(0, 2.0*np.pi, 100))
 
# rychlá Fourierova transformace
f = fft(t)
 
# posun signálu ve frekvenční oblasti
f2 = fftshift(f)
 
# vykreslení výsledného signálu
plt.plot(f2)
 
# uložení grafu s průběhem signálu
plt.savefig("fft_6.png")

Nyní je skutečně stejnosměrná složka umístěna doprostřed grafu:

Výpočet FFT pro funkci Hannova okna

V předchozím článku jsme se seznámili  takzvanými okny, což jsou funkce používané pro předzpracování signálů s využitím konvoluce. Ovšem samotné okno je vlastně digitálním signálem, takže si pro něj můžeme vypočítat jeho obraz ve frekvenční oblasti. Celý výpočet i s vizualizací výsledků může vypadat následovně:

import scipy.signal as signal
from scipy.fft import fft, fftshift
 
import numpy as np
import matplotlib.pyplot as plt
 
# filtr
win = signal.windows.hann(50)
 
# rychlá Fourierova transformace
f = fft(win, 1024) / (len(win)/2.0)
 
# rozsah normalizovaných frekvencí na x-ové ose
freq = np.linspace(-0.5, 0.5, len(f))
 
response = np.abs(fftshift(f / abs(f).max()))
 
# vykreslení výsledného signálu
plt.plot(freq, response)
 
# uložení grafu s průběhem signálu
plt.savefig("fft_7.png")
 
# zobrazení grafu
plt.show()

Výsledek:

Převod na decibely (logaritmická jednotka)

Poměrně často se setkáme se zobrazením spektrálního průběhu takovým způsobem, že na vertikální osu jsou vyneseny decibely (popř. například u zvuků jsou i frekvence vyneseny na logaritmickou stupnici). Převod výsledku FFT (po posunu) na decibely je snadný, protože můžeme využít „vektorovou“ operaci numpy.log10:

response = 20 * np.log10(np.maximum(response, 1e-10))

Zakomponování tohoto výpočtu do skriptu je stejně triviální:

import scipy.signal as signal
from scipy.fft import fft, fftshift
 
import numpy as np
import matplotlib.pyplot as plt
 
# filtr
win = signal.windows.hann(50)
 
# rychlá Fourierova transformace
f = fft(win, 1024) / (len(win)/2.0)
 
# rozsah normalizovaných frekvencí na x-ové ose
freq = np.linspace(-0.5, 0.5, len(f))
 
response = np.abs(fftshift(f / abs(f).max()))
 
# dB
response = 20 * np.log10(np.maximum(response, 1e-10))
 
# vykreslení výsledného signálu
plt.plot(freq, response)
 
# uložení grafu s průběhem signálu
plt.savefig("fft_8.png")
 
# zobrazení grafu
plt.show()

Výsledkem bude nyní tento graf:

Úprava grafu odstraněním první hodnoty z vypočteného vektoru

Výsledný graf ve skutečnosti není zcela symetrický, i když by tak měl vypadat. K „dokonalosti“ (což je v oblasti zpracování signálů dosti nepatřičný termín) postačuje odstranit první prvek z posunutého signálu. Ovšem stejně tak musíme odstranit tento prvek z x-ové osy, tedy z vektoru freq:

import scipy.signal as signal
from scipy.fft import fft, fftshift
 
import numpy as np
import matplotlib.pyplot as plt
 
# filtr
win = signal.windows.hann(50)
 
# rychlá Fourierova transformace
f = fft(win, 1024) / (len(win)/2.0)
 
# rozsah normalizovaných frekvencí na x-ové ose
freq = np.linspace(-0.5, 0.5, len(f))
 
response = np.abs(fftshift(f / abs(f).max()))
 
# dB
response = 20 * np.log10(np.maximum(response, 1e-10))
 
# vykreslení výsledného signálu
plt.plot(freq[1:], response[1:])
 
# uložení grafu s průběhem signálu
plt.savefig("fft_9.png")
 
# zobrazení grafu
plt.show()

Opět se pochopitelně podíváme na výsledek této operace:

Fourierova transformace obdélníkového signálu

Zajímavé bude zjistit, jak bude vypočtena a vizualizována rychlá Fourierova transformace u obdélníkového signálu. U takového signálu totiž můžeme předpokládat rozvinutí do nekonečného množství frekvencí, takže nám nezbývá, než si tento předpoklad ověřit (ovšem u signálu s konečným množstvím samplů):

import scipy.signal as signal
from scipy.fft import fft, fftshift
 
import numpy as np
import matplotlib.pyplot as plt
 
# signal
t = np.zeros(50) + np.ones(50) + np.zeros(50)
 
# rychlá Fourierova transformace
f = fft(t, 1024) / (len(t)/2.0)
 
# rozsah normalizovaných frekvencí na x-ové ose
freq = np.linspace(-0.5, 0.5, len(f))
 
# vykreslení výsledného signálu
plt.plot(f)
 
# uložení grafu s průběhem signálu
plt.savefig("fft_A.png")
 
# zobrazení grafu
plt.show()

Zobrazení výsledku výpočtu FFT, ovšem bez posunu pomocí funkce fftshift:

Posun FFT funkcí fftshift

Na x-ové ose je ovšem většinou vhodné zobrazit korektní rozsah frekvencí (absolutně v Hz nebo relativně v rozsahu –1/2 až 1/2), takže se nevyhneme výpočtu těchto frekvencí a použití funkce fftshift. Taktéž provedeme převod na decibely:

# rozsah normalizovaných frekvencí na x-ové ose
freq = np.linspace(-0.5, 0.5, len(f))
 
response = np.abs(fftshift(f / abs(f).max()))

Upravený skript bude nyní vypadat následovně:

import scipy.signal as signal
from scipy.fft import fft, fftshift
 
import numpy as np
import matplotlib.pyplot as plt
 
# signal
t = np.zeros(50) + np.ones(50) + np.zeros(50)
 
# rychlá Fourierova transformace
f = fft(t, 1024) / (len(t)/2.0)
 
# rozsah normalizovaných frekvencí na x-ové ose
freq = np.linspace(-0.5, 0.5, len(f))
 
response = np.abs(fftshift(f / abs(f).max()))
 
# vykreslení výsledného signálu
plt.plot(freq, response)
 
# uložení grafu s průběhem signálu
plt.savefig("fft_B.png")
 
# zobrazení grafu
plt.show()

Nyní je výsledný průběh názornější:

Diskrétní kosinová transformace (DCT)

V předchozím textu jsme se zabývali především dopřednou i zpětnou rychlou Fourierovou transformací, ovšem knihovna SciPy nabízí i výpočty dalších typů transformací. Při zpracování statických snímků i videa se velmi často setkáme s diskrétní kosinovou transformací, kterou si přiblížíme v navazujících kapitolách.

Diskrétní kosinová transformace (DCT – Discrete Cosine Transform) patří do skupiny metod provádějících takzvané transformační kódování nad diskrétním (vzorkovaným) jednorozměrným či vícerozměrným signálem. Do stejné skupiny patří i optimální KLT (Karhunen-Loéve Transform – Karhunen-Loéveho transformace) či známá FFT (Fast Fourier Transform – rychlá diskrétní Fourierova transformace), kterou jsme se (ve vztahu ke knihovně SciPy) zabývali v předchozích kapitolách. Základem prakticky všech v praxi používaných transformačních kódování je v případě zpracování obrazu, který můžeme pro tyto účely považovat za dvourozměrný vzorkovaný signál, nalezení korelace mezi sousedními, popř. i vzdálenějšími pixely – jedná se o takzvanou mezipixelovou redundanci. Při práci s videem se kromě mezipixelové redundance uplatní i redundance mezisnímková.

Většinou se při transformacích provádí převod (či mapování, resp. zobrazení) zpracovávaného signálu z časové (prostorové) oblasti do oblasti frekvenční, protože existuje předpoklad, že například obrazy reálných předmětů neobsahují mnoho energie ve vyšších frekvencích a je tedy vhodné shromáždit co největší množství relevantních dat do relativně malého množství koeficientů. Kódování by mělo v důsledku vést k celkovému snížení počtu bitů nesoucích vizuální informaci (psychovizuální redundance). V dalším textu se již soustředíme na popis diskrétní kosinové transformace.

Existuje více typů diskrétních kosinových transformací, které jsou v literatuře označovány římskými číslicemi jako DCT I až DCT IV. My si v dalším textu představíme „pouze“ DCT II. Ta je zdaleka nejpoužívanější a tvoří základ pro většinu praktických úloh, které DCT nějakým způsobem využívají. Dvourozměrná diskrétní kosinová transformace typu II je v současnosti použita například v mnoha souborových formátech i formátech určených pro přenos videa. Například se jedná o standardy JPEG, MPEG 1 (do značné míry zastaralé, podobně jako další dva zmíněné standardy), MPEG 2, MPEG 4, JVT (H.263), H.264, H.261, HEVC(H.265) VP9 atd. Sice existují i jiné transformace vhodné pro některé úlohy, například velmi populární vlnková transformace (Wavelet Transform), KLT atd., ty však z několika důvodů (implementačních, obchodních i patentových) prozatím DCT v praxi nepřekonaly.

V knihovně SciPy je možné výpočet DCT provést s využitím funkce, která se jmenuje jednoduše – dct a nalezneme ji ve stejném podbalíčku, jako funkce určené pro výpočet Fourierovy transformace:

dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, orthogonalize=None)
    Return the Discrete Cosine Transform of arbitrary type sequence x.
 
    Parameters
    ----------
    x : array_like
        The input array.
    type : {1, 2, 3, 4}, optional
        Type of the DCT (see Notes). Default type is 2.
    n : int, optional
        Length of the transform.  If ``n < x.shape[axis]``, `x` is
        truncated.  If ``n > x.shape[axis]``, `x` is zero-padded. The
        default results in ``n = x.shape[axis]``.
    axis : int, optional
        Axis along which the dct is computed; the default is over the
        last axis (i.e., ``axis=-1``).
    norm : {"backward", "ortho", "forward"}, optional
        Normalization mode (see Notes). Default is "backward".
    overwrite_x : bool, optional
        If True, the contents of `x` can be destroyed; the default is False.
    workers : int, optional
        Maximum number of workers to use for parallel computation. If negative,
        the value wraps around from ``os.cpu_count()``.
        See :func:`~scipy.fft.fft` for more details.
    orthogonalize : bool, optional
        Whether to use the orthogonalized DCT variant (see Notes).
        Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.

Základní vlastnosti DCT

Diskrétní kosinová transformace má některé vlastnosti, které jsou důležité a mnohé velmi užitečné v praxi. Patří mezi ně:

1. Energie jednorozměrného, dvourozměrného či n-rozměrného signálu je po provedení diskrétní kosinové transformace většinou soustředěna v několika málo koeficientech. Toho se využívá při následném kvantování vypočtených koeficientů (využito například ve známém rastrovém formátu JPEG).
2. Výsledek DCT transformace je blízký optimální KLT transformaci, především z důvodu minimalizace takzvané zbytkové korelace. DCT je však implementačně mnohem jednodušší než KLT (ostatně základy výpočtu DCT si uvedeme hned v navazující kapitole).
3. Existují rychlé implementace přímé i inverzní DCT. Při optimalizacích se většinou využívá symetrie výpočtů a takzvaná separabilita (například dvourozměrná DCT lze vypočítat s využitím série jednorozměrných DCT).
4. Při implementaci je možné použít jak aritmetiku v plovoucí řádové čárce (FP – Floating Point), tak i aritmetiku prováděnou v pevné řádové čárce (FX – Fixed Point).
5. Rekurzivní struktura – vhodné při implementaci vícerozměrných DCT.
6. DCT má ovšem i některé nectnosti, především blokovou strukturu, která může být viditelná například při nastavení velkých komprimačních faktorů (malého koeficientu Q) při komprimaci pomocí JPEGu. Totéž platí i při zápisu videa do některého z formátů MPEG (ukázky si uvedeme v navazujícím článku).

Jednorozměrná diskrétní kosinová transformace (1D DCT)

Jednorozměrná diskrétní kosinová transformace typu II (zkráceně ji v dalším textu budeme označovat jako 1D DCT) slouží ke zpracování jednorozměrného diskrétního (tj. vzorkovaného) signálu, získaného například ze zvukové stopy CDčka. Nemusí se však jednat pouze o celočíselné vzorky, protože pomocí DCT lze zpracovat i reálný či dokonce komplexní signál. Každou hodnotu tohoto signálu označme symbolem s(n), kde n je pozice (index) vzorku v signálu. Po aplikaci jednorozměrné diskrétní kosinové transformace získáme sadu transformovaných koeficientů, které budeme značit t(k).

Vztah pro výpočet 1D DCT má následující tvar:

kde N značí řád/velikost DCT (v praktických aplikacích typicky 8 či 16) a hodnota indexu k leží v rozsahu 0..N-1. Konstanta c(k) je rovna (1/N)1/2 v případě, že k=0, v ostatních případech platí c(k)=(2/N)1/2:

S využitím jednorozměrné diskrétní kosinové transformace tedy můžeme zpracovat vždy maximálně N vzorků, v případě jejich většího počtu se vzorec periodicky opakuje pro dalších N vzorků vstupního signálu. Teoreticky je možné za N zvolit délku celého vstupního signálu, v praxi se však tato činnost neprovádí, protože přináší některé implementační problémy (relativně pomalá operace, velká spotřeba paměti, nevýhodnost z hlediska ztrátové komprimace atd.).

Ukázka výpočtu bázových funkcí DCT

Samotný výpočet bázových funkcí DCT je možné relativně snadno naprogramovat. Následující program (psaný pro změnu v céčku) vytvoří trojici obrázků s bázovými funkcemi použitými v jednorozměrné diskrétní kosinové transformaci:

/*----------------------------------------------------------------------------- */
/* Tento program vytvori obrazky s prubehem bazovych funkci pouzitych           */
/* v jednorozmerne diskretni kosinove transformaci (1D DCT).                    */
/*----------------------------------------------------------------------------- */
 
#include <assert.h>
#include <math.h>
#include <png.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
/* datova struktura popisujici pixmapu */
typedef struct {
    unsigned int   width;
    unsigned int   height;
    unsigned char *pixels;
} pixmap;
 
/*----------------------------------------------------------------------------- */
/* vytvoreni pixmapy ve formatu truecolor                                       */
/*----------------------------------------------------------------------------- */
pixmap *createPixmap(unsigned int width, unsigned int height) {
    pixmap      *px;
    unsigned int size;
 
    assert(width > 0);
    assert(height > 0);
    /* alokace pameti pro hlavicku pixmapy */
    px = (pixmap *)malloc(sizeof(pixmap));
    /* kontrola alokace pameti */
    assert(px);
    px->width = width;
    px->height = height;
    size = 3 * width * height;
    /* alokace pameti pro pixely */
    px->pixels = (unsigned char *)malloc(sizeof(unsigned char) * size);
    /* kontrola alokace pameti */
    assert(px->pixels);
    return px;
}
 
/*----------------------------------------------------------------------------- */
/* ulozeni pixmapy do externiho souboru                                         */
/*----------------------------------------------------------------------------- */
int savePixmap(const pixmap *pixmap, const char *file_name) {
    FILE                *fout;
    const unsigned char *p = pixmap->pixels;
    int                  code = 0;
    int                  scanline;
    png_structp          png_ptr = NULL;
    png_infop            info_ptr = NULL;
 
    char *title = NULL;
 
    if (p == NULL || file_name == NULL) {
        return -1;
    }
 
    /* open file for writing in binary mode */
    fout = fopen(file_name, "wb");
    if (fout == NULL) {
        fprintf(stderr, "Could not open file %s for writing\n", file_name);
        code = 1;
        goto FINALISE;
    }
 
    /* initialize write structure */
    png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
    if (png_ptr == NULL) {
        fprintf(stderr, "Could not allocate write struct\n");
        code = 1;
        goto FINALISE;
    }
 
    /* initialize info structure */
    info_ptr = png_create_info_struct(png_ptr);
    if (info_ptr == NULL) {
        fprintf(stderr, "Could not allocate info struct\n");
        code = 1;
        goto FINALISE;
    }
 
    /* setup Exception handling */
    if (setjmp(png_jmpbuf(png_ptr))) {
        fprintf(stderr, "Error during png creation\n");
        code = 1;
        goto FINALISE;
    }
 
    png_init_io(png_ptr, fout);
 
    /* write header (8 bit colour depth) */
    png_set_IHDR(png_ptr, info_ptr, pixmap->width, pixmap->height, 8,
                 PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
                 PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
 
    /* set the title */
    if (title != NULL) {
        png_text title_text;
        title_text.compression = PNG_TEXT_COMPRESSION_NONE;
        title_text.key = "Title";
        title_text.text = title;
        png_set_text(png_ptr, info_ptr, &title_text, 1);
    }
 
    png_write_info(png_ptr, info_ptr);
 
    /* write image data */
    for (scanline = 0; scanline < pixmap->height; scanline++) {
        png_write_row(png_ptr, p);
        p += pixmap->width * 3;
    }
 
    /* end write */
    png_write_end(png_ptr, NULL);
 
FINALISE:
    if (fout != NULL) {
        fclose(fout);
    }
    if (info_ptr != NULL) {
        png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
    }
    if (png_ptr != NULL) {
        png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
    }
 
    return code;
}
 
/*----------------------------------------------------------------------------- */
/* ulozeni barvy pixelu na souradnice x, y                                      */
/*----------------------------------------------------------------------------- */
void putpixel(const pixmap *pix, int x, int y, int r, int g, int b) {
    unsigned char *p;
    assert(pix);
    assert(pix->pixels);
    assert(x >= 0 && x < pix->width);
    assert(y >= 0);
    assert(y < pix->height);
    p = pix->pixels + 3 * (x + y * pix->width);
    *p++ = r;
    *p++ = g;
    *p = b;
}
 
/*----------------------------------------------------------------------------- */
/* vykresleni horizontalni usecky                                               */
/*----------------------------------------------------------------------------- */
void hline(const pixmap *pix, int x1, int x2, int y, int r, int g, int b) {
    int i;
    for (i = x1; i <= x2; i++) {
        putpixel(pix, i, y, r, g, b);
    }
}
 
/*----------------------------------------------------------------------------- */
/* vykresleni vertikalni usecky (povoleno prohozeni souradnic y1 <=> y2)        */
/*----------------------------------------------------------------------------- */
void vline(const pixmap *pix, int x, int y1, int y2, int r, int g, int b) {
    int i;
    if (y1 < y2) {
        for (i = y1; i <= y2; i++) {
            putpixel(pix, x, i, r, g, b);
        }
    } else {
        for (i = y2; i <= y1; i++) {
            putpixel(pix, x, i, r, g, b);
        }
    }
}
 
/*----------------------------------------------------------------------------- */
/* vykresleni bazovych funkci jednorozmerne DCT                                 */
/*----------------------------------------------------------------------------- */
void draw1d_dct(const pixmap *pix, int n) {
#define PI 3.1415927
    int u;
    /* pozicovani prubehu bazovych funkci v pixmape */
    int   xoffset = 25 * 8 / n, yoffset = 50;
    int   xdelta = 50 * 8 / n, ydelta = 100;
    float amplitude = 40.0;
 
    /* vykreslit vsech "n" bazovych funkci */
    for (u = 0; u < n; u++) {
        int x;
        /* pozice horizontalni osy pro danou bazovou funkci */
        int   y0 = yoffset + ydelta * u;
        float y;
        /* horizontalni osa */
        hline(pix, 0, pix->width - 1, y0, 255, 128, 0);
        /* vypocet vsech bodu bazove funkce */
        for (x = 0; x < n; x++) {
            float alfa = PI * (2.0 * x + 1.0) * u / (2.0 * n);
            /* hodnota bazove funkce pro dane "x" */
            y = amplitude * cos(alfa);
            /* printf("%d %d -> %f\n", u, x, alfa); */
            /*  vykresleni vertikalni hodnoty */
            vline(pix, xoffset + xdelta * x, y0, y0 - y, 0, 255, 128);
            /* zvyrazneni koncoveho bodu pomoci horizontalni carky */
            hline(pix, xoffset + xdelta * x - 3, xoffset + xdelta * x + 3, y0 - y, 0,
                  255, 128);
        }
    }
}
 
void process_dct(int n) {
#define FNAME_LENGTH 50
    char    full_name[FNAME_LENGTH];
    pixmap *pix = createPixmap(400, n * 100);
 
    draw1d_dct(pix, n);
    snprintf(full_name, FNAME_LENGTH, "dct_%02d.png", n);
    savePixmap(pix, full_name);
}
 
/*----------------------------------------------------------------------------- */
/* hlavni funkce konzolove aplikace                                             */
/*----------------------------------------------------------------------------- */
int main(int argc, char **argv) {
    printf("\nprocessing...\n");
 
    process_dct(8);
    process_dct(16);
    process_dct(32);
 
    printf("\ndone!\n");
    return 0;
}
 
/*----------------------------------------------------------------------------- */
/* finito                                                                       */
/*----------------------------------------------------------------------------- */

Vypočtené a vizualizované výsledky budou vypadat následovně:

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/
94. The Fourier Analysis –The Fast Fourier Transform (FFT) Method
    https://www.electronics-lab.com/article/the-fourier-analysis-the-fast-fourier-transform-fft-method/
95. Understanding the output of FFT
    https://howthefouriertran­sformworks.com/understanding-the-output-of-an-fft/
96. What factors go into choosing DCT type 1 over type 2 or type 3 etc?
    https://dsp.stackexchange­.com/questions/31611/what-factors-go-into-choosing-dct-type-1-over-type-2-or-type-3-etc
97. The Discrete Cosine Transform (DCT)
    https://ccrma.stanford.edu/~jos/mdft/Dis­crete_Cosine_Transform_DCT­.html
98. Eulerův vzorec
    https://cs.wikipedia.org/wi­ki/Euler%C5%AFv_vzorec
99. Komplexní analýza
    https://cs.wikipedia.org/wi­ki/Komplexn%C3%AD_anal%C3%BDza
100. What is the difference between a Fourier transform and a cosine transform?
     https://dsp.stackexchange­.com/questions/13/what-is-the-difference-between-a-fourier-transform-and-a-cosine-transform